TOI-837 b is a Young Saturn-sized Exoplanet with a Massive 70 𝑀 ⊕ Core

We present an exhaustive photometric and spectroscopic analysis of TOI-837, a F9/G0 35 Myr young star, hosting a transiting exoplanet, TOI-837 b, with an orbital period of ∼ 8 . 32 d. Utilising data from TESS and ground-based observations, we determine a planetary radius of 0 . 818 + 0 . 034 − 0 . 024 𝑅 J for TOI-837 b. Through detailed HARPS spectroscopic time series analysis, we derive a Doppler semi-amplitude of 34 . 7 + 5 . 3 − 5 . 6 m s − 1 , corresponding to a planetary mass of 0 . 379 + 0 . 058 − 0 . 061 𝑀 J . The derived planetary properties suggest a substantial core of approximately 70 𝑀 ⊕ , constituting about 60% of the planet’s total mass. This finding poses a significant challenge to existing theoretical models of core formation. We propose that future atmospheric observations with JWST could provide insights into resolving ambiguities of TOI-837 b, offering new perspectives on its composition, formation, and evolution.


INTRODUCTION
Characterising the properties of young exoplanets (< 1 Gyr) at different stages is crucial to understanding the evolution and populations of exoplanets.These planets are often elusive to detection using indirect techniques like transit and radial velocity (RV) methods, as strong stellar activity generates stellar signals that often overshadows their signals in photometric and spectroscopic data.Recently, tens of young transiting exoplanets have been discovered (e.g., Bouma et al. 2020;David et al. 2018;Hobson et al. 2021;Martioli et al. 2021;Mann et al. 2021;Rizzuto et al. 2020;Barragán et al. 2022b).This success is largely attributed to missions like K2 (Howell et al. 2014) and the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) that provide thousands of light curves of young stars.
After identifying a young transiting exoplanet, the next step typically involves spectroscopic observations to observe the star's RV variations caused by the planet to measure its mass.However, separating signals from the planet and star in radial velocity data remains challenging.Gaussian Processes (GPs) have emerged as a favoured method for modelling activity-induced radial velocities, offering a flexible way to represent stochastic variations, such as the quasiperiodic stellar signals (see e.g., Aigrain & Foreman-Mackey 2023).Multiple authors have used GPs, and variations of them such as multidimensional GPs (see e.g., Rajpaul et al. 2015;Barragán et al. ★ oscar.barrragan@physics.ox.ac.uk 2022a), to detect RV planetary signals on spectroscopic time series of active young stars (e.g., Barragán et al. 2019bBarragán et al. , 2023;;Mantovan et al. 2023;Mallorquín et al. 2023;Zicher et al. 2022).Nevertheless, although widely embraced by the community, there exists a notable concern that the GP activity model may inadvertently mask or alter potential planetary signals (see e.g., discussions in Ahrer et al. 2021;Rajpaul et al. 2021).A recent example of this comes from Suárez Mascareño et al. (2021).The authors, utilising GPs, discovered that the planets orbiting the young V 1298 Tau exhibited significantly higher density than anticipated.However, these detections have been challenged by the exoplanet community.Blunt et al. (2023) suggested that Suárez Mascareño et al. (2021)'s detections may have been affected by over-fitting, thereby questioning the reliability of their RV Doppler detections.This is crucial, especially in cases where the measurements are in tension with models.Bouma et al. (2020, hereafter B20) reported the discovery and validation of the young transiting exoplanet around TOI-837.The star is a young (∼ 35Myr) F9/G0 dwarf star in the southern open cluster IC 2602.TOI-837's main identifiers and parameters are given in Table 1.They used TESS cycle 1, the GAIA mission, multi-band photometry and spectroscopic observations to validate the planetary nature of a transit signal.The transit object corresponds to a planet with a radius of ∼ 0.77  J and a period of 8.32 d.In this manuscript, we present a spectroscopic follow-up and exhaustive analysis of TOI-837, to characterise further the nature of TOI-837 b.This paper is structured as follows: The photometric and spectroscopic data of  (K) 5995 ± 79 to this data are outlined in Section 3. A discussion of the findings is provided in Section 4, and the paper concludes with a summary of the key outcomes in Section 5.This manuscript is part of a series of papers under the project GPRV: Overcoming stellar activity in radial velocity planet searches funded by the European Research Council (ERC, P.I. S. Aigrain).We acknowledge that while this paper was being reviewed, we learned that Damasso et al. (2024) were performing an independent analysis of the public HARPS data of TOI-837.Should the reader wish to compare their results with ours, it would serve to further validate/test the findings presented in this manuscript.
2 TOI-837 DATA 2.1 TESS data TOI-837 (TIC 460205581) was observed by the TESS mission on cycles 1, 3 and 5.During cycle 1, TOI-837 was observed in Sectors 10 and 11 from 2019 March 26 to 2019 May 20 with a 2 min cadence.The TESS Science Processing Operations Center (SPOC; Jenkins et al. 2016) transit search (Jenkins 2002;Jenkins et al. 2010Jenkins et al. , 2020) ) discovered a transiting signal with a period of 8.3 d in TOI-837's light curve.This was announced in the TESS SPOC Data Validation Report (DVR; Twicken et al. 2018;Li et al. 2019)  It should be noted that TOI-837 is situated in a relatively dense stellar region.Utilising Gaia DR2 data, B20 identified two sources brighter than  = 16 within the same TESS pixel as TOI-837: TIC 847769574 ( = 14.6), which was argued a possible co-moving companion within the IC 2602 moving group, and TIC 460205587 ( = 13.1), a background giant star.B20 employed speckle imaging and photometric ground-base observations (refer to Section 2.2) to rule out potential false-positive scenarios and confirm that the observed transits occur on TOI-837, thereby validating the planetary nature of the transiting signal.More recently, Behmard et al. (2022), using Gaia DR3 data, suggested that TIC 847769574 is likely gravitational bound to TOI-837.This companion is consistent with being an M dwarf, with an estimated mass of 0.393 ± 0.034, ⊙ and a projected separation of approximately 330 AU, suggesting that the TOI-837 system may be binary, with TIC 847769574 being TOI-837 (B).However, the findings from B20 continue to support the scenario of the planet transiting TOI-837 (A).For simplicity, this manuscript will continue to refer to the planet as TOI-837 b instead of TOI-837 (A) b.We note that any extra flux coming from contaminating stars within TESS pixels (named crowding) is corrected for the PDCSAP light curves by the SPOC (Stumpe et al. 2012).Therefore, no further processing is needed for the analyses in this manuscript where we use PDCSAP light curves.
During cycles 3 and 5, TOI-837 was reobserved in Sectors 37 and 38 (from 2021 April 02 to 2021 May 26, with 2 min cadence) and Sectors 63 and 64 (2023 Mar 10 to 2023 May 04, with 20 s cadence), respectively.We downloaded the Presearch Data Conditioning Simple Aperture Photometry (PDCSAP; Smith et al. 2012;Stumpe et al. 2012Stumpe et al. , 2014) ) light curves for TOI-837 from the Mikulski Archive for Space Telescopes (MAST).Figure 1 shows the TOI-837's light curves for all the three TESS cycles.TOI-837 will be re-observed by TESS in Sector 90 during the first semester of 2025.
To perform the transit analysis we flattened the TESS light curves (see Sect. 3.2).The TESS light curves for TOI-837 exhibit variations outside of transits, possibly due to stellar or instrumental variability.The TESS light curves were detrended using the open-source code citlalicue (Barragán et al. 2022a).In essence, citlalicue employs GPs, as implemented in george (Ambikasaran et al. 2015), to model light curve fluctuations outside of transits.We provided citlalicue with normalised light curves and the ephemeris for the transit signal.To focus on removing low-frequency signals, we used data binned at 3-hour intervals and masked all transits when fitting the GP with a Quasi-Periodic kernel (as in Ambikasaran et al. 2015).We applied an iterative maximum Likelihood optimisation coupled with a 5-sigma clipping technique to identify the best model for light curve variations outside of transits.Subsequently, we divided the entire light curves based on this model, resulting in a flattened curve showing only transit signals.It is worth noting that each TESS cycle was detrended individually.The detrended light curves for all TESS cycles are presented in Figure 1.

Ground-based photometric data
B20 performed multi-band ground based follow-up of TOI-837.For full details on the acquisition and reduction we refer the reader to that paper.Four transit events were monitored using the 0.36-meter telescope at El Sauce Observatory, situated in Chile's Río Hurtado Valley.The data collection occurred during specific nights: the Cousins-R band was used on April 1 and 26 of 2020, the Cousins-I band on May 21, 2020, and the Johnson-B band on June 14, 2020.TOI-837 was also observed by the Antarctic Search for Transiting Exoplanets (ASTEP) telescope at the Concordia base on the Antarctic Plateau, on the dates of May 12, May 29, June 14, and June 23 of 2020.We downloaded all these public data to include them in our analyses.

HARPS spectroscopic observations
We acquired 75 high-resolution (R ≈ 115 000, wavelength range of 380 to 690 nm) spectra of TOI-837 with the High Accuracy Radial velocity Planet Searcher (HARPS; Mayor et al. 2003) spectrograph.The instrument is mounted at the 3.6 m ESO telescope at La Silla Observatory, Chile.
The observations were carried out between October 03 2022 and January 16 2023, as part of the observing program 0110.C-4341(A) (PI: Yu).The typical exposure time per observation was 1800 s, except for the observations before October 25 where the exposure time was set to 900 s.This produced spectra with a typical signal-tonoise (S/N) of 50 (30 for the 900s exposures) at 550 nm.
We post-processed the 1-D spectra produced by the official data reduction software (DRS).All the spectra (  , ) were first continuum normalised by RASSINE (Cretignier et al. 2020b), an upper envelope method fitting the spectra continua  (  , ).The normalised spectra time series  (  , ) = (  , )/ (  , ) was then post-processed by YARARA (Cretignier et al. 2021).Given the limited S/N of the observations, we restricted the cleaning by YARARA to cosmic rays and tellurics which are the only notable features at such a level of flux precision.From this step onwards, four observations were rejected by YARARA based on an anomalous number of detected outliers.Those spectra were all very low S/N observations (< 25).Once spectra were corrected of systematics, the normalised flux spectra  (  ) were scaled back to absolute flux units spectra (  ) by using a reference continuum 1  ref () where the bolometric flux of the spectra was preserved during the scaling: This step is equivalent in practice to the colour correction usually performed on the order-by-order CCFs to remove any dependency with airmass or weather condition (Lovis 2007;Cretignier 2022).We extracted the RVs and activity indicators by a Gaussian fit on CCFs obtained using the G2 mask of the HARPS DRS since for fast-rotating stars, the empirical line list (see e.g., Cretignier et al. 2020a) obtained from the observations could be imperfect due to a large number of blends.The full width at half maximum (FWHM) of the CCFs shows a consequent broadening with a value around FWHM ∼ 26 km s −1 in agreement with the fast rotation rate of the star (see Sect. 3.1).Since the CCFs are considerably broadened, we also extracted the 1 Sometimes referred as a "reference colour".
RVs with a Generalised normal distribution (GND) as introduced in Heitzmann et al. (2021).GND has for advantage to contains an extra shape parameter  that affects the kurtosis of the profile in order to produce a more "top-hat flatten" behaviour.When  = 2, GND is fully equivalent to a Gaussian.We extracted the RVs using both a completely free GND model (() time dependent) and a more rigid one with the shape parameter  fixed to the mean value (() = 2.8).The RVs obtained with these alternative parametrizations are fully consistent with the RVs obtained with the Gaussian fit.Our HARPS RV measurements have a typical error of 10 m s −1 and a RMS of 100 m s −1 .Table 2 lists the extracted HARPS RV, RV GND, RV GND_2.8,FWHM, Bisector span (BIS), calcium lines S-index ( HK ), and Hydrogen-alpha (H  ) time series.

Periodograms
As a first check to test the information contained in our HARPS spectroscopic time series, we ran a General Lomb-Scargle (GLS; Zechmeister & Kürster 2009) periodogram on them.Figure 2 shows the periodogram of all the time series.We can see that all of the spectroscopic time series peak around 3 and 1.5 d, which correspond to the rotational period of the star and its first harmonic.The peaks at the stellar rotation period and its first harmonic suggest the presence of active regions on the stellar surface (e.g., Aigrain et al. 2012).No substantial peak corresponding to the orbital period of TOI-837 b is evident in the raw RVs.We also check for signals at longer periods (up to 100 d, not shown in Fig. 2) and did not find any significant peak in any time series.

Atmospheric parameters
We used the HARPS spectra to derive the atmospheric parameters of TOI-837.As this star is a fast rotator, we analysed the spectra using spectral synthesis.Two different methods were employed, FASMA-synthesis2 and PAWS3 , as well as two different sets of spectra, the DRS spectra and the YARARA-processed spectra (see Sect. 2.3).Each set of spectra was shifted in the lab frame and coadded and the analyses were performed on this co-added spectrum.
FASMA-synthesis made use of the MARCS atmospheric models (Gustafsson et al. 2008) and the radiative transfer code MOOG4 .Microand macroturbulent velocities were calculated through the calibration relations mentioned in Tsantaki et al. (2013) and Doyle et al. (2014), respectively.More details on this method can be found in Tsantaki et al. (2020).
PAWS uses the functionalities from iSPEC (Blanco-Cuaresma 2019) and is described in more detail in Freckelton et al. (2024).We used the Kurucz atmospheric models in our analysis (Kurucz 1993) and skipped the initial step of the equivalent width analysis as this is less reliable for fast rotators where the spectral lines are broadened.
These four analyses provided values for the effective temperature, the surface gravity, the metallicity, and the projected rotational velocity.All values are in agreement with their 1-sigma errors.For our adopted parameters, listed in Table 1, we choose to take a weighted average of all the results.The high value of the projected rotational velocity (16.8 km/s) confirms the fast-rotating nature of TOI-837.
All other parameters are close to Solar values.

Mass, radius, and age
To calculate the fundamental stellar parameters, we made use of the isochrones (Morton 2015) package.Our stellar models came from the MESA Isochrones and Stellar Tracks (MIST -Dotter 2016) and we used the nested sampler MultiNest (Feroz et al. 2019) for our likelihood analysis.As input data, we used the photometric  magnitudes B, V, J, H, Ks, W1, W2, W3, the Gaia DR3 parallax, and the spectroscopically derived effective temperature and metallicity.
In a similar way as described in Mortier et al. (2020), we ran the code four times, each time changing the spectroscopic input to the individually obtained results and keeping the other data the same.Narrow bounds were set on the stellar age, between 30 and 46 Myr, following B20, which helps constrain the stellar mass given strong correlations between mass and age for young stars.
After checking that all parameters from the individual runs agreed within errors, the final parameters and errors, listed in Table 1, were obtained from the median and 16th/84th percentile of the combined posterior distributions of the four runs.The mass and radius values are, unsurprisingly, close to Solar values.From these, we can recalculate a much more precise value of the surface gravity.It is lower than, but consistent with, the spectroscopically derived surface gravity.Using the stellar radius and the projected rotational velocity, we can furthermore also place an upper limit on the stellar rotation period.We find that the rotation period has a maximum value of 3.17 ± 0.04 d.This is fully compatible with the values obtained with TESS and spectroscopic time series (see Sect. 3.3).

Transit analysis
For all the subsequent planet analyses we used the code pyaneti (Barragán et al. 2019a;Barragán et al. 2022a).In all our runs we sample the parameter space with 250 walkers using the Markov chain Monte Carlo (MCMC) ensemble sampler algorithm implemented in pyaneti (Barragán et al. 2019a) which is based on Foreman-Mackey et al. (2013).The posterior distributions are created with the last 5000 iterations of converged chains.We thinned our chains by a factor of 10 giving a distribution of 125 000 points for each sampled parameter.
To speed up the transit modelling, we only model data chunks spanning a maximum of three hours on either side of each transit mid-time.TESS and ground-base photometry were acquired with short cadence (< 2 min).Therefore we can assume that this data can be described by instantaneous evaluations of the transit models and we do not need model re-sampling (see e.g., Gandolfi et al. 2018).In total we have six different photometric data sets named: TESS 2 min data (cycles 1 and 3), TESS 25s data (cycle 5), ASTEP, and El Sauce R, I, and B bands.
To model the transits of TOI-837 b, we need to set priors for the following set of parameters: time of mid-transit,  0 ; orbital period,  orb ; orbital eccentricity,  and angle of periastron,  using the polar parametrisation

√
cos  ★ and √  sin  ★ ; scaled planetary radius  p / ★ ; the stellar density,  ★ ; and the the limb darkening parameters  1 and 𝑞 2 for each band (following Mandel & Agol 2002; Kipping  2013, models and parametrisations).We assume a black object transit and we sample for a single  p / ★ for all bands (B20 showed that there is no variation of  p / ★ in different bands).The scaled semimajor axis, / ★ , is recovered from  ★ and Kepler's third law (see e.g., Winn 2010).The model also includes photometric jitter term per data set to penalise the likelihood.We start by assuming that the orbit is circular, so we fix For the rest of parameters we set wide informative uniform priors as listed in Table 3.
To start our transit analysis, it is worth mentioning that TOI-837 b's transit looks V-shaped, indicating that the transit could be grazing (where the grazing condition is given when  > 1 −  p / ★ ).B20 reported  = 0.94 ± 0.013 and  p / ★ = 0.08 ± 0.01, suggesting that the transit is either grazing or nearly grazing.We note that B20 performed this analysis using the available data at the time: TESS cycle 1 and the ground-based data.We performed an analysis using exactly this data set, and we obtained fully consistent results of  = 0.923 +0.019 −0.011 and  p / ★ = 0.079 +0.006 −0.003 .The small differences with B20 can be explained by the different ways to deal with the detrending and/or different parametrisations.
We then repeated the analysis including the new data of the TESS cycles 3 and 5.We obtained more precise values of  = 0.921 +0.014 −0.009 and  p / ★ = 0.080 +0.004 −0.002 .These values are still consistent with a grazing transit.Figure 3 shows a correlation plot between ,  ★ , and  p / ★ .We can see that the grazing condition covers only the low probability tail of the posterior distribution.Therefore, we are able to put strong constraints on the inferred planetary radius.
It is also worth mentioning that the stellar density recovered from the transit analysis,  ★ = 1.58 ± 0.17 g cm −3 , is fully consistent with the stellar parameters derived in Section 3.1.This suggests that the planetary orbit is nearly circular.To test this further, we performed another transit analysis where we used a Gaussian prior on the stellar density using the value in Table 1.We also sample for the orbit eccentricity using the polar parametrisation

√
sin  ★ .We inferred √  cos  ★ = 0.08 +0.12 −0.15 and √  sin  ★ = 0.10 +0.39 −0.35 that relates to an eccentricity of  = 0.10 +0.17 −0.06 .These results are consistent with a nearly circular orbit for TOI-837 b.It is also worth mentioning that the model with a circular orbit is strongly preferred over the model with eccentricity with a difference of Akaike Information Criterion (Δ AIC) of 52.We then conclude that the photometric data strongly prefers a model assuming a circular orbit.Figure 4 shows TOI-837 b's transits for the different data sets with the respective inferred models.

Characterising the stellar signal
Because of its youth, TOI-837 has strong stellar signals in its light curve and RVs time series.We performed several analyses of the light curves and spectroscopic time series to analyse the time scales over which the stellar signal evolves.

Stellar signal in the TESS light curves
We first performed a onedimensional GP regression of the TESS light curves to study the behaviour of the stellar signal at the different TESS cycles.Since we are interested in modelling the long-term evolution of the light curve, we binned the data into 3 hour chunks and masked the transits out of the light curve.
We model the covariance between two data points at times   and   for each time-series as where  is an amplitude term, and  ,  is the Quasi-Periodic (QP) kernel given by whose hyperparameters are:  GP , the GP characteristic period;  p , the inverse of the harmonic complexity; and  e , the long-term evolution timescale.
We sample 5 parameters in each run: four GP hyperparameters (,  GP ,  e ,  p ), and a jitter term included in the Gaussian Likelihood.We set wide uniform priors for all the parameters: for  e we set a uniform prior between 0.1 and 100 d, for  p between 0.1 and 5, and for  GP between 2.5 and 3.5 d.We sample the parameter space using the built-in MCMC sampler in pyaneti using the same criteria as in Sect.3.2.
Table 4 shows the inferred  GP ,  e , and  p hyperparameters for all the light curves.Figure 5 shows the inferred posterior distributions for the  GP ,  e , and  p parameters for all the modelled time series.Figure 6 shows the TESS light curves with the inferred GP model.The first thing we observe is that the recovered  GP is consistent within the maximum rotational period derived in Section 3.1.However, the rotational period obtained from TESS Cycle 1 is consistent with Cycles 3 and 5 just at 3 sigma.This could be evidence of differential rotation.This can also be caused by the relatively short evolution time-scale recovered that is smaller than the inferred period (see e.g., Rajpaul et al. 2015).We note that all the hyperparameters seem consistent between the different light curves.This is unusual for young stars given that each TESS cycle happens two years after another and the evolution of stellar activity can manifest as different signals that can be explained by different GP hyperparameters (e.g., Barragán et al. 2021).

Stellar signal in the activity indicators
We performed an individual analysis of each activity indicator as the one presented in Barragán et al. (2023).We analyse the FWHM, BIS,  HK , and H  spectroscopic time series.Note that, in comparison with Section 3.3.1,we are now facing a more difficult signal Equilibrium temperature (c)  eq () F [] refers to a fixed value , U [, ] to an uniform prior between  and , N [, ] to a Gaussian prior with mean  and standard deviation , and J [, ] to the modified Jeffrey's prior as defined by Gregory (2005, eq. 16).
Inferred parameters and errors are defined as the median and 68.3% credible interval of the posterior distribution.
Assuming an albedo of zero.
Derived using  p =   p  −2 p . Derived using sampled parameters following Southworth et al. (2007).characterisation given that the spectroscopic time series are not as well sampled as the TESS light curves.
In Section 2.3.1 we show how all the spectroscopic time-series periodograms show significant peaks around 3 and 1.5 d, which are related to the stellar rotational period and its first harmonic, respec-tively.To better characterise the periodic signals and analyse further signal complexity, we performed a GP regression on each time series using the 1D modelling setup as the one presented in Sect.3.3.1.We sample for 6 parameters in each run, the kernel hyperparameters (,  GP ,  e ,  p ), an offset, and a jitter term.The priors used were uni- form with the ranges: for  e between 0.1 and 100 d (that corresponds to the HARPS observing window),  p between 0.1 and 2, and for  GP between 2.5 and 3.5 d.
Table 4 and Figure 5 summarise the inferred GP hyperparameters.Figure 6 shows the spectroscopic time-series data together with the inferred GP model.The first thing we note is that all activity indica-  tors can recover a period that is consistent with the stellar rotation period.This is consistent with the results found in the periodogram analysis (even the double peak observed in the H  periodogram is presented in the recovered posterior, see Sect.2.3.1).We were expecting a better agreement between the hyperparameters recovered from the spectroscopic time series (see Barragán et al. 2023).From Figure 6 we can see by eye that the recovered process in each time series looks significantly different.The processes recovered to explain FWHM and BIS Span present variations that are shorter than the time sampling of the HARPS data, this behaviour can be a result of overfitting (see Blunt et al. 2023).We can also see that the Solid-coloured lines show the corresponding inferred signal coming from our GP regression for TESS Cycles 1 (brown), 3 (pink), and 5 (grey).Bottom panel: TOI-837 spectroscopic time-series for FWHM, BIS Span,  HK , and H  .The corresponding measurements are shown with black circles with error bars with a semi-transparent error bar extension accounting for the inferred jitter.Solid-coloured lines show the corresponding inferred signal coming from our GP regression, while light-coloured shaded areas show the one and two-sigma credible intervals of the corresponding GP model.can be caused by relatively high white noise that does not allow to characterise the complexity of the underlying signal.
These results suggest that the data cadence is not good enough to characterise the stellar signal in the activity indicators.The rotation period of the star of 3 days is short compared with our sampling strategy of taking one point per night.Three points may not be enough to sample the complexity of the stellar signal, especially if this signal is expected to vary significantly in short time scales, as suggested by the TESS photometry.Therefore, our GP regression may not be able to characterise in full the stellar signal.Furthermore, the rotation period of the star is close to the integer of 3 d.This has as a consequence that observations do not sample properly the phase of the stellar rotation given that observations happen at similar times during the whole observing run.This is worsened in the case in which the stellar signal evolves drastically between each period.We conclude that for our spectroscopic time series, we can constrain the rotation period of the star.However, any further complexity pattern that may exist cannot be characterised by our dataset.Because of this lack of information on the activity indicators, we conclude that they do not bring enough information of the stellar signal into a multidimensional GP regression analysis.In Appendix A we show the results and a discussion on the results of applying a multidimensional GP framework to this data set.

RV analysis
We also performed onedimensional GP modelling to the RVs time series.We first ran a model identical to the GP analyses in Sect.3.3.2.This means that we sample for the QP kernel hyperparameters (,  GP ,  e ,  p ), and a jitter term.As a mean function, we included just an offset (with no Keplerian model yet).For the GP hyperparameters, we use the same priors as the ones described in Sect.3.3.2.We run this analysis for the Gaussian, GND and GND_2.8RV time series.For the Gaussian RVs we recover the hyperparameters  GP = 3.00 ± 0.02 d,  e = 42 +25 −16 d, and  p = 0.39 +0.22 −0.13 ; for GND  GP = 3.00±0.01d,  e = 69 +22 −24 d, and  p = 0.42 +0.26 −0.14 ; and GND_2.8  GP = 3.00±0.01d,  e = 67 ± 23 d, and  p = 0.42 +0.25  −0.14 .We can see that the hyperparameters are consistent between the different RV flavours.We can also see that besides the period, the hyperparameters do not compare with any of the hyperparameters obtained from the other spectroscopic time series (see Table 4).
We then repeated the analysis but this time we added a mean function that includes an offset and a Keplerian model to account for the Doppler effect of TOI-837 b.For the planet model we assume a circular orbit, we set Gaussian priors for the planet ephemeris based on the analysis in Sect.3.2.For the Doppler semi-amplitude we set a uniform prior between 0 and 100 m s −1 .We inferred hyperparameters fully consistent with the no-planet models.We also recovered a Doppler semi-amplitude of  = 36.6± 6.5 m s −1 for the Gaussian RVs,  = 34.7 ± 5.6 m s −1 for the GND RVs, and  = 34.6 ± 5.5 m s −1 for the GND_2.8time series.All these values are fully consistent between them.
When comparing all models, the ones with the planetary signal are strongly favoured with a Δ AIC > 20.It is also worth mentioning that the recovered jitter term for the model with the planet signal is significantly smaller than the one without a Keplerian signal.The jitter values reduce from 38 m s −1 to 30 m s −1 , from 37 m s −1 to 27 m s −1 , and 36 m s −1 to 26 m s −1 , for the Gaussian, GND and GND_2.8 time series, respectively.This is also evidence that the model that better adjusts the data is the one where the planetary signal is included.The GND_2.8 RV time series yields the most precise Doppler semi-amplitude and exhibits the lowest jitter among the RV series evaluated.Based on these findings, we determine that the GND_2.8series represents the optimal dataset for our analyses.Consequently, we will focus exclusively on this time series in our subsequent work.

Final joint model
Based on the analyses presented in this section, our final model for the photometric data of TOI-837 is the transit model described in Section 3.2, together with the RV model to the GND_2.8data described in Sect.3.4.For completeness, we ran a final joint model of photometry and RVs to characterise TOI-837 b.The whole set of sampled parameters and priors are shown in Table 3.
We also explored a solution where we allowed for an eccentricity.We inferred an orbital eccentricity of 0.10 +0.06 −0.05 .However, the model including a circular orbit is slightly preferred with Δ AIC = 4.This suggests that the orbit of TOI-837 b is (close to) circular and the data we have is not enough to measure any small deviation from it.
Figure 7 shows the RV time series, while Fig. 8 shows the phasefolded Doppler signal.Table 3 shows the inferred sampled parame-ters, defined as the median and 68.3% credible interval of the posterior distribution.Table 3 also shows the derived planetary and orbital parameters.

RV detection tests
We know a priori that TOI-837 b exists and its nature is consistent with being planetary (B20).Therefore, we expect that there is a Doppler signal larger than zero that is consistent with the ephemeris of the transiting signal.If we assumed that TOI-837 b has the same properties as similar planets with its radius, we estimated a mass of ∼ 0.11  J that would generate a Doppler signal of ∼ 10 m s −1 (Chen & Kipping 2017).Due to its youth (∼ 35 Myr), we would expect a rather lower density object with a significantly smaller mass (e.g., Baraffe et al. 2008;Fortney et al. 2007).Our detected Doppler semi-amplitude is 34.7 +5.3  −5.6 m s −1 which translates to a mass of 0.379 +0.058 −0.061  J for TOI-837 b.This mass together with the inferred radius of 0.818 +0.034 −0.024  J results in a planetary density of 0.89 +0.20  −0.18 g cm −3 .This unexpected high density on young giant planets has been found previously (Suárez Mascareño et al. 2021), but the RV detection of such planets are challenged by the community (e.g., Blunt et al. 2023).In this section, we test the reliability of our detection of the induced Doppler signal on TOI-837.

Planet signal in periodograms
As shown in Section 2.3.1 and Figure 2 there is no evidence of the planetary signal in RV GLS periodograms of TOI-837.This is expected given that any planetary signal would be dwarfed by the relatively high stellar activity signal.
We then took the residuals that we obtained from the no-planet onedimensional GP analyses presented in Sect.3.4 for the tree different RV time series and we applied a GLS periodograms that are shown in the bottom panel of Figure 2. We can see that the peaks that correspond to the stellar signal disappeared from the periodograms, indicating that the GP model was able to remove the intrinsic stellar signal in the RV time series.Furthermore, a significant peak appeared close to the orbital period of the planet.This result is unexpected.GPs have been proven to engulf planetary signals, especially when the latter ones are not included in the model (e.g., Ahrer et al. 2021;Rajpaul et al. 2021).We speculate that if a relatively strong coherent signals is present in our data set, it may survive a GP regression if it evolves with significantly different time scales than the stellar signal.This assumption is consistent with the recovered Doppler signal of TOI-837 b.
To explore further the existence of the planetary signal in the periodogram, we ran a ℓ1 periodogram (Hara et al. 2017) on the raw GND_2.8RV time series.Briefly, the ℓ1-periodogram is a variation of the Lomb-Scargle periodogram methodology that has reduced sensitivity towards outliers and non-Gaussian noise.Figure 9 shows the ℓ1-periodogram applied to the raw RV time-series of TOI-837.The three most significant peak happens at 1.5, 8.31 and 3.09 d that correspond to the first harmonic of the rotation period of the star, the planet orbit, and the rotation period of the star, respectively.This suggests that there is a strong deep-seated signal that coincides with the orbital period of TOI-837 b.This result supports our previous speculation that there is evidence of strong coherent signal in the RVs that is consistent with our TOI-837 b's recovered Doppler signal.

Cross-validation
The next test we perform is cross-validation to check for overfitting as suggested by Blunt et al. (2023).To do so, we employ a data partition- ing strategy similar to the principles of k-fold cross-validation (see e.g., Fazekas & Kovacs 2024) to systematically divide our dataset into distinct subsets.Specifically, we first randomised the dataset to ensure that each partition is representative of the overall data distribution.Subsequently, we segmented the dataset into five partitions.In each iteration, we use four out of these five partitions (constituting approximately 80% of the data) for the modelling, while the remaining partition (approximately 20% of the data) is masked out.This approach guarantees that each unique subset serves as the test set exactly once, eliminating any potential overlap between masked-out data across all iterations.Figure 10 shows the five different subsamples created with this approach.
We then perform a RV analysis (accounting for the planet signal) as the one presented in Sect.3.4 for each sub-sample, we call each one of these runs as sub-test , where  runs between 1 and 5. Figure 10 shows the data and recovered model for all cases.The recovered Doppler semi-amplitude for each case is fully consistent with the value reported in Sect.3.4, but with slightly larger error bars, as expected due to the smaller number of data modelled.We can see that, by eye, the masked-out points fall well within the limits of the inferred models in each case, suggesting that the predictive model in each case can explain the masked-out observations (see also discussion in Luque et al. 2023).We can also see that the predictive model for each sub-set (black lines) falls well within the confidence intervals of the predictive distribution of the model of the full data set (blue lines and shaded region).These results suggest that our model does not suffer overfitting.
To assess further the reliability of our cross-validation, we perform a Two-sample Kolmogorov-Smirnov (KS;Kolmogorov 1933;Smirnov 1948) test.Briefly, the KS test is a non-parametric statistical test that is utilised to validate the hypothesis that two subsets of data originate from the same continuous distribution.This test compares the empirical cumulative distribution functions (ECDFs) of two samples without making any assumptions about the underlying distribution of data.The maximum absolute difference between the ECDFs of the two samples, known as the KS statistic, serves as the basis for evaluating the null hypothesis that the two samples are drawn from the same distribution.A low KS statistic, accompanied by a high p-value, indicates a failure to reject the null hypothesis, suggesting that the differences between the two samples could be attributed to random variation, and hence, they can be considered as coming from the same distribution.
We apply the KS test to each one of our sub-tests in which we compare the residuals of the modelled data with the residuals of the masked-out data (i.e., masked-out points minus predictive distribution at their respective times).If the KS test proves that both samples are explained from the same distribution, then we can conclude that the predictive model explains well the masked-out data.This implies that the inferred model has the ability to explain unseen data, thus, it does not overfit.Figure 10 shows the ECDF for the modelled and masked-out data for each one of the sub-tests.The p-values of the five sub-tests are significantly larger than 0.05, therefore, we can conclude that both samples are consistent with being described by the same underlying distribution.These results give us confidence that our RV model for TOI-837 does not suffer from overfitting.

Injection tests
Previous research indicates that complex models can generate false planet-like RV time series, particularly in stars with high activity levels (e.g., Rajpaul et al. 2016).To check if this is our case, we performed injection tests similar to the ones presented in Barragán et al. (2019b); Zicher et al. (2022).We used citlalatonac (Barragán et al. 2022a) to simulate synthetic RV time series at the same time-stamps as our HARPS data.We first took the predictive distribution for the stellar signal that we obtained for TOI-837 in Sect.3.4.We then added correlated noise using a squared exponential kernel with a length scale of one day, and the same amplitude as the jitter term obtained from the real data to simulate the red noise in our data.To finish, we also added white noise for each synthetic observation according to the nominal measurement uncertainty of each HARPS datum.We created 100 different realisations of stellar-like signals.
For each stellar-like signal, we created 3 different types of RV time series injecting signals with Doppler semi-amplitudes of 0 (assuming there is no planet), 10 (expected signal of this planet if it were inflated), and 35 m s −1 (the recovered signal).This leads to a total of 300 synthetic RV time series with similar noise properties and the same time sampling as the real data.
We model each synthetic dataset using a one-planet model and onedimensional GP configuration as described in Sect.3.4.We plot the recovered posterior distribution of the Doppler semi-amplitude of the planet for all the runs in Figure 11.
We first check the percentage of the simulations in which we can claim a detection.For this purpose, we consider that a detection has occurred if the median of the posterior is larger than 3-, where  is half the interval between the 16.5 th and 83.5 th percentiles.We found that we can claim a detection of 3%, 28%, and 100% for the synthetic time series with injected signals of 0, 10, and 35 m s −1 , respectively.We then revise which percentage of these detections is consistent with our detection in the real data set.We check what fraction of the time the recovered semi-amplitude in the synthetic time series is within 1- of the recovered posterior distribution from the real data set.We found that this condition is filled in 0%, 1%, and 62% of the cases for the 0, 10, and 35 m s −1 time series, respectively.These results suggest that if the real Doppler signal presented in our RV time-series was 10 m s −1 we would detect the 35 m s −1 signal with a probability of 1%.The 62% recovered in the 35 m s −1 case is consistent with the expected value taking into account Poisson counting errors.
For the remainder of the paper, we will assume that our mass and radius measurements of TOI-837 b are reliable.Below we describe the physical implications of this for the planet's properties.

TOI-837b's properties
With a mass of 0.379 +0.058  −0.061  J and radius of 0.818 +0.034 −0.024  J , TOI-837 b's position in the mass-radius diagram is highlighted in Fig. 12. Figure 12 also shows a mass-radius diagram for giant exoplanets (0.5 <  J < 2  J and 0.1 <  p < 13  J ) detected with a precision better than 30% in radius and mass.Considering that the evolution of planets is influenced by their distance from their host star, we only present planets that orbit within a semi-major axis range of 0.08 to 1 AU.This selection is informed by the semi-major axis of TOI-837 b, ∼ 0.09 AU, and is supported by models that forecast a comparable evolution for planets situated within this range of semimajor axes (see Fortney et al. 2007).We also over plot the Fortney et al. (2007)'s models with similar parameters to those of TOI-837 b.This corresponds to models for planets around a solar-like star, with an age of 32 Myr, at a distance of 0.1 AU with cores of 25 (pink), 50 (orange) and 100 (green)  ⊕ .We performed a radial basis function interpolation on the aforementioned Fortney et al. (2007)'s models to find the core mass that best describes the mass and radius of TOI-837 b.We found that the position and age of TOI-837 b would be consistent with a young planet with a core mass of 68  ⊕ , assuming a 50-50% ice-silicate core (black line in Fig. 12).This core mass corresponds to ∼ 60% of TOI-837 b total mass of ∼ 125  ⊕ .
From figure 12, we can see that if well there are planets that fall close to TOI-837 b in the mass vs radius diagram, they all are older than 1 Gyr and they should be compared with older planet's models.This can be done from figure 12 where we also show the Fortney et al. (2007)'s models for planets aged 3 Gyr with orbital semi-major axis of 0.1 AU.The first thing we see is that the older counterparts of TOI-837 b are consistent with planets with cores ≲ 50  ⊕ .If TOI-837 b itself was older, we would estimate a significantly smaller core of ∼ 50  ⊕ , our estimation of ∼ 70  ⊕ is a pure consequence of the expected inflated state, hence larger radius, given its youth.In this line of thought, TOI-837 b stands out as the only planet with a semi-major axis of ∼ 0.1 AU consistent with a ∼ 70  ⊕ core.
We then relaxed our assumptions and proceed to compare TOI-837 b with planets with similar mass and radius, but without constrain on their semi-major axes.We found that K2-60 b (Eigmüller et al. 2017), HD 149026 b (Sato et al. 2005), and TOI-1194b (Wang et al. 2023) are denser than TOI-837 b and are consistent with core masses between 50−100 ⊕ despite being significantly older (> 1 Gyr).This suggests that planets with relatively similar masses to TOI-837 b and large cores can exist.The question now is what physical mecha-nisms can create such massive cores.Johnson & Li (2012) theorised that a star's high metallicity might reflect a primordial circumstellar disc with elevated dust content that could potentially enhance core formation.Both, HD 149026 b and TOI-1194 b have super-solar atmospheric metallicities that are consistent with this model.In fact, Bean et al. (2023) observed a significant atmospheric metal enrichment on the atmosphere of HD 149026 b that is explained by a bulk heavy element abundance of 66% of the planetary mass.However, the solar metallicity of TOI-837 and K2-60 does not support this hypothesis, suggesting that other factors may also play a critical role in the formation of massive cores.Boley et al. (2016)  result from the merging of tightly packed orbiting inner planets that formed during the initial phases of the circumstellar disc's evolution.
Another possibility is that TOI-837 b has a less massive core, but the models are biased towards larger radii for young planets.Coreaccretion models of planetary evolution predict planets of < 100 Myr to be at the early stages of their contraction phase, showing very large radii and low densities (Baraffe et al. 2008;Fortney et al. 2007, see also 1 Myr models in Fig. 12).However, Fortney et al. (2007)'s models do not include a formation mechanism and can be arbitrarily large and hot at very young ages (see Marley et al. 2007).This can lead to giant planet's models with radii several tenths of a Jupiter radius larger than one computed taking into account formation mechanisms.This would lead as a consequence to estimate more massive cores for a given planet's mass and radius.For the remainder of the paper we will assume that TOI-837 b can be described by Fortney et al. (2007)'s models assuming a core of 70  ⊕ .
It is also worth to compare TOI-837 b with the properties of other exoplanets that are younger than 100 Myr.According to the NASA exoplanet archive (Akeson et al. 2013), there are five such exoplanets documented with precise measurements of age, mass, and radius.These planets include AU Mic b and c (22 Myr, Plavchan et al. 2020;Martioli et al. 2021;Zicher et al. 2022;Donati et al. 2023), V1298 Tau b and e (23 Myr, Suárez Mascareño et al. 2021;Finociety et al. 2023), andHD 114082 b (15 Myr, Zakhozhay et al. 2022).The planets orbiting AU Mic have properties that suggest they are inflated when compared to exoplanets orbiting older stars.This aligns with existing theoretical models that predict younger planets might exhibit such inflated characteristics due to active and intense processes in their early developmental stages.Although these planets are smaller than Neptune, with radii smaller than 4  ⊕ , they showcase a distinct set of physical processes which might differ significantly from those observed in larger exoplanets like TOI-837 b.Moreover, the exoplanets V1298 Tau b and e, along with HD 114082 b, are consistent with giant planets, as illustrated in Fig. 12.Interestingly, these planets are denser than what current theoretical models predict for their age and type (see discussions in Suárez Mascareño et al. 2021;Finociety et al. 2023;Zakhozhay et al. 2022).This tendency suggests that existing theoretical models and/or mass/radius young exoplanet measurements require further dedicated studies.In case our measurement methods are correct, this implies a reconsideration of planetary formation and evolution theories.

Atmospheric characterisation perspectives
Atmospheric observations could play a crucial role in differentiating among various potential formation mechanisms for TOI-837 b.According to the mass-metallicity relationship empirically established by Thorngren et al. (2016), our mass measurement of TOI-837 b suggests a predicted metal fraction ratio between the planet and its host star of  planet / star = 14.7 ± 2.5.With TOI-837's stellar metallicity measured at [Fe/H] = 0.01 ± 0.04, this translates into an estimated bulk metallicity fraction for TOI-837 b of 0.21 ± 0.04.Bean et al. (2023) presents a relationship between atmospheric and bulk metallicity fractions.If we assume that TOI-837 b follows this trend, we could expect its atmospheric metallicity to surpass solar values by a factor of 5 − 8.This significant enhancement in atmospheric metallicity would not only support the substantial core mass inferred from our observations but also provide insight into the planet's formation conditions and subsequent evolutionary history.
To predict the observability of the atmosphere of TOI-837 b with JWST, we leverage packages picaso5 and pandexo6 .We forward model the emission spectrum of TOI-837 b using picaso (Batalha et al. 2019), based on the parameters derived in this work.We make use of the integrated 1D climate modelling for the temperaturepressure (TP) profile and chemical abundance solution (Mukherjee et al. 2023), including all opacities available in picaso.The TP profile includes a convective zone only in the deepest atmospheric layers.We then use pandexo to estimate the signal-to-noise of the dayside emission signal of TOI-837 b with JWST/NIRSpec G395H; the predicted observations using three eclipses are shown in Fig. 13, assuming solar metallicity and C/O.At least two eclipse observations are needed to distinguish between a 1× and 10× solar metallicity scenario.Using the chi-square statistic for hypothesis testing, we would expect to distinguish between these cases to a significance of 11.1 with two eclipse observations (totalling 10.01 hrs), rising to 17.1 with three eclipses (Fig. 13).Constraining the atmospheric metallicity will help to break the degeneracy with interior composition and the implications for formation.

CONCLUSIONS
Our investigation into the TOI-837 system and its intriguing companion, TOI-837 b, unveils a young Saturn-sized exoplanet that defies conventional expectations with its unexpected massive core.Our exhaustive analysis of data from TESS ground-based observations, and HARPS spectroscopic enabled us to determine TOI-837 b's radius at 0.818 +0.034 −0.024  J and mass at 0.379 +0.058 −0.061  J , translating to a density of 0.89 +0.20  −0.18 g cm −3 .This density together with its age and distance to the star suggest a core mass of approximately 70  ⊕ , accounting for 60% of the planet's total mass.Such a substantial core within a relatively young planetary body presents a challenging scenario for  (2007).Dotted and dashed lines extend these predictions to planets aged 1 Myr and 3 Gyr, respectively.Black line shows the interpolated model for a 32 Myr old and 68  ⊕ core.We also highlighted other giant exoplanets younger than 100 Myr with red squares.The methodology for generating this diagram follows the approach outlined in Barragán et al. (2018).
current models of planet formation and core accretion, especially due to the relatively low stellar metallicity.
The unique characteristics of TOI-837 b underscore the urgency for advanced atmospheric characterisation.Eclispe observations with JWST could to offer unparalleled insights into the composition of TOI-837 b.A measurement of the planetary atmospheric bulk metal fraction will potentially elucidate the true nature of its significant core.Such future studies are crucial for breaking the current degeneracies in planet composition models and could revolutionise our understanding of planetary formation.
We also leave open the possibility that our RV detection could not be accurate.Despite the tailored campaign with a high cadence of observations, the apparent fast evolution of the strong stellar signal could generate biased measurements of the Doppler semi-amplitude.We showed with several statistical tests that our RV planetary detection is robust, nevertheless, further observations could help us to test this further.For example, tailored near Infrared RV campaigns, where stellar activity is less significant, with instruments such as the Near Infra Red Planet Searcher (NIRPS) could help to assert our detection in the optical.information in the activity indicators.Therefore, we have chosen to adopt the more conservative onedimensional regression approach for the main paper.This decision is driven by our judgement regarding the use of an activity indicator that alone does not constrain the stellar signal.
and designated by the TESS Science Office as TESS Object of Interest (TOI; Guerrero et al. 2021) TOI-837.01 (hereafter TOI-837 b).

Figure 1 .
Figure 1.TESS light curves for TOI-837 in cycles 1, 3 and 5 (from top to bottom).TESS data are shown with grey points with the out-of-transit variability model over-plotted in red.The resulting flattened light curve is presented with blue points.TOI-837's transit positions are marked with a blue triangle.All plots have the same -range to ease comparison of data and signals.

Figure 2 .
Figure 2. GLS periodograms of the spectroscopic time series.GLS periodograms of the spectroscopic time series are presented in panels a) to f).Panel a) displays periodograms for Gaussian (solid line), GND (dotted line), and GND_2.8 (dashed line) extracted RVs, which appear overlapped.Panel f) depicts a GLS periodogram for these RVs after removing the stellar model (see Sect 3.4).The horizontal dashed line indicates the 1% False Alarm Probability (FAP).The dot-dashed and dotted black lines represent the stellar rotation period and its first harmonic, respectively, while the vertical red line marks the orbital period of TOI-837 b.

Figure 3 .
Figure 3. Correlation plot for the main parameters of the transit analysis.

Figure 4 .
Figure 4. Phase-folded light curves of TOI-837 b for different data sets.Nominal observations are shown in light grey.Coloured circles represent 10min binned data.Transit models are shown with a solid black line.

Figure 6 .
Figure 6.Top panel: TOI-837's TESS light curves and GP inferred models.The corresponding measurements (3 h binned data) are shown with black circles.Solid-coloured lines show the corresponding inferred signal coming from our GP regression for TESS Cycles 1 (brown), 3 (pink), and 5 (grey).Bottom panel: TOI-837 spectroscopic time-series for FWHM, BIS Span,  HK , and H  .The corresponding measurements are shown with black circles with error bars with a semi-transparent error bar extension accounting for the inferred jitter.Solid-coloured lines show the corresponding inferred signal coming from our GP regression, while light-coloured shaded areas show the one and two-sigma credible intervals of the corresponding GP model.

Figure 7 .Figure 8 .
Figure 7. TOI-837's RV time-series after being corrected by inferred offsets.Top panel: RV data together with full, stellar and planetary signal inferred models; RV data with stellar signal model subtracted; and RV residuals.Measurements are shown with black circles, error bars, and a semi-transparent error bar extension accounting for the inferred jitter.The solid lines show the inferred full model coming from our multi-GP, light-shaded areas showing the corresponding GP model's one and two-sigma credible intervals.We also show the inferred stellar (dark blue line) and planetary (light green line) recovered signals with an offset for better clarity.Middle panel: RV data with stellar signal model subtracted and planetary model.Bottom panel: RV residuals.

Figure 10 .
Figure 10.Top panel: Cross-validation results with RV data sub-sets.Dark circles represent modelled data, red circles indicate masked-out data.Black solid line shows the recovered model for each sub-test.We also show the full model and 3 −  confidence interval of the full model recovered from the whole data set (same inferred model as Figure7) with a blue line and light shaded areas, respectively.Bottom panel: Empirical Cumulative Distribution Functions (ECDF) for the residuals of the modelled (black) and masked-out (red) data for all sub-tests.Each sub-plot shows the corresponding KS statistics.

Figure
Figure Doppler semi-amplitude distributions from injection tests.The blue line shows results from original data; other coloured lines from synthetic data.Labels indicate injected values.

Figure 12 .
Figure12.Mass vs Radius Diagram.Grey data points with error bars show exoplanets with mass and radius measurements with a precision of 30% or better and orbital semi-major axis between 0.08 and 1 AU.Data taken from the NASA Exoplanet Archive as of March 2024 (https://exoplanetarchive.ipac.caltech.edu/;Akeson et al. 2013).The location of TOI-837 b is highlighted by a blue circle.Solid lines depict theoretical models predicting the mass-radius relationship for planets of varying core masses, assuming an age of 32 Myr and a distance of 0.1 AU from a Sun-like star, based onFortney et al. (2007).Dotted and dashed lines extend these predictions to planets aged 1 Myr and 3 Gyr, respectively.Black line shows the interpolated model for a 32 Myr old and 68  ⊕ core.We also highlighted other giant exoplanets younger than 100 Myr with red squares.The methodology for generating this diagram follows the approach outlined inBarragán et al. (2018).

Figure 13 .
Figure13.The predicted emission spectrum of TOI-837 b.The red data points show the spectrum attainable with three eclipse observations, using JWST/NIRSpec G395H, assuming a solar composition in both metallicity and C/O ratio.For comparison, the forward-modelled spectra assuming 10× and 100× metallicity are also shown (with solar C/O).

Table 1 .
Main identifiers and parameters for TOI-837.

Table 2 .
HARPS spectroscopic measurements.The full version of this table is available in a machine-readable format as part of the supplementary material.

Table 3 .
Modelled and derived parameters for TOI-837 b.

Table 4 .
Recovered hyperparameters for onedimensional GP regression for different stellar time series.