TESS and ESPRESSO discover a super-Earth and a mini-Neptune orbiting the K-dwarf TOI-238

The number of super-Earth and mini-Neptune planet discoveries has increased significantly in the last two decades thanks to transit and radial velocity surveys. When it is possible to apply both techniques, we can characterise the internal composition of exoplanets, which in turn provides unique insights on their architecture, formation and evolution. We performed a combined photometric and radial velocity analysis of TOI-238 (TYC 6398-132-1), which has one short-orbit super-Earth planet candidate announced by NASA's TESS team. We aim to confirm its planetary nature using radial velocities taken with the ESPRESSO and HARPS spectrographs, to measure its mass and to detect the presence of other possible planetary companions. We carried out a joint analysis by including Gaussian processes and Keplerian orbits to account for the stellar activity and planetary signals simultaneously. We detected the signal induced by TOI-238 b in the radial velocity time-series, and the presence of a second transiting planet, TOI-238 c, whose signal appears in RV and TESS data. TOI-238 b is a planet with a radius of 1.402$^{+0.084}_{-0.086}$ R$_{\oplus}$ and a mass of 3.40$^{+0.46}_{-0.45}$ M$_{\oplus}$. It orbits at a separation of 0.02118 $\pm$ 0.00038 AU of its host star, with an orbital period of 1.2730988 $\pm$ 0.0000029 days, and has an equilibrium temperature of 1311 $\pm$ 28 K. TOI-238 c has a radius of 2.18$\pm$ 0.18 R$_{\oplus}$ and a mass of 6.7 $\pm$ 1.1 M$_{\oplus}$. It orbits at a separation of 0.0749 $\pm$ 0.0013 AU of its host star, with an orbital period of 8.465652 $\pm$ 0.000031 days, and has an equilibrium temperature of 696 $\pm$ 15 K. The mass and radius of planet b are fully consistent with an Earth-like composition, making it likely a rocky super-Earth. Planet c could be a water-rich planet or a rocky planet with a small H-He atmosphere.


Introduction
The discovery and characterisation of exoplanets is one of the most exciting fields in astronomy today.The diversity and complexity of these worlds challenge our understanding of planetary formation and evolution, as well as offer potential clues for the origin and distribution of life in the universe.To fully explore these aspects, it is essential to obtain accurate measurements of the physical properties of exoplanets, such as their masses, radii, densities, atmospheres, orbits, and interactions with their host stars.
The total number of known exoplanets keeps increasing.With main contributions from transit satellites such as the Kepler Space Telescope (Borucki et al. 2010;Howell et al. 2014) and the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015), as well as radial velocity surveys, the number of confirmed detections has surpassed 5500 1 (Christiansen 2022).Despite this ⋆ Based [in part] on Guaranteed Time Observations collected at the European Southern Observatory under ESO programme 1102.C-0744.by the ESPRESSO Consortium. 1 NASA Exoplanet Archive (Akeson et al. 2013) promising evolution in terms of detection, the characterisation of these objects is a matter that has not been addressed with the same level of success.From the large sample of confirmed exoplanets, only ∼20% have their true dynamical mass constrained, and only ∼10% of the complete sample have a measurement of their density with uncertainties lower than 33%, with one of the main limitations being the uncertainty in the planetary mass.This limits studies of planetary formation or atmospheric characterisation, for which a precise mass measurement is required (Batalha et al. 2019).
There is an observed bimodal distribution of the size of small exoplanets (Fulton et al. 2017).Atmospheric mass loss mechanisms such as photo evaporation (Owen & Wu 2017) or corepowered mass loss (Ginzburg et al. 2018) are both able to adequately reproduce it, assuming the planets have rocky cores.This has led to the interpretation that super-Earths and sub-Neptunes have formed accreting only dry condensates within the water ice line, and their difference in sizes is attributed to the absence or presence of extended volatile-rich atmospheres.However, given their bulk densities, another possibility is that sub-Neptune planets are water-rich worlds (Léger et al. 2004).The most promising way of discriminating between both scenarios is through atmospheric characterisation of exoplanets around the valley of the distribution, which first requires producing a significant sample of transiting low-mass exoplanets with precise masses.
The most widely used method to measure exoplanet masses is radial velocity (RV), which relies on detecting the Doppler shift induced by the gravitational tug of an orbiting planet on its host star.The amplitude of this shift depends on several factors, such as the orbital period, eccentricity, inclination, and the planet-to-star mass ratio.In the past years, RV instruments such as HARPS (Mayor et al. 2003), HARPS-N (Cosentino et al. 2012), and more recently, CARMENES (Quirrenbach et al. 2014) and ESPRESSO (Pepe et al. 2021), have demonstrated that it is possible to detect planets with masses similar to the mass of the Earth using RVs (e.g., Anglada-Escudé et al. 2016;Astudillo-Defru et al. 2017;Suárez Mascareño et al. 2020, 2023).However, this is only possible in the case of low-mass stars.For planets orbiting solar-type and K-type stars, we continue to be restricted to the detection and mass measurement of super-Earths (rocky planets with masses between 1 and 10 M ⊕ Mayor et al. 2011).
TESS is a NASA mission that was launched in April 2018 with the primary goal of finding transiting planets around bright, nearby stars Ricker et al. (2015).TESS observes a large fraction of the sky using four cameras that cover 24x96 degree fields of view.Each field is observed for about 27 days, with a cadence of 2 minutes for selected targets and 200 seconds for full-frame images.TESS image data are reduced and analyzed by the Science Processing Operations Center (SPOC) at NASA Ames Research Center.The SPOC pipeline performs photometry on the images to produce light curves (LC), which then removes instrumental and systematic effects and performs transit searches on them (Jenkins et al. 2016).This pipeline also performs various validation tests to assess the reliability of transit signals and to rule out false positives caused by instrumental or astrophysical effects.The most promising candidates are then labelled as TESS Objects of Interest (TOIs) and are made available to the community for further follow-up observations.TOI-238 (TYC 6398-132-1, TIC 9006668) is a bright (V = 10.75 mag) K-dwarf.It was observed by the TESS satellite in sectors 2, 29 and 69 (August-September of 2018, 2020and 2023, respectively).It is suspected to host a super-Earth (Radius ∼ 1.6 R ⊕ ) in a very short orbit (orbital period ∼ 1.27 days).This result was recently validated by Mistry et al. (2023).
In this paper we present a combined analysis of the TESS data with RV measurements from the ultra-stable spectrographs ESPRESSO and HARPS that allow us to obtain a high-precision mass determination for the transiting candidate TOI-238.01.Our combined analysis reveals the presence of a transiting sub-Neptune with a period of ∼8.46 days.

TESS observations
TESS is an all-sky photometric transit survey with the main objective of finding planets smaller than Neptune orbiting bright stars.Such stars are amenable to follow-up observations that aim to determine planetary masses and atmospheric compositions.In its primary mission, TESS has conducted high-precision photometry of more than 200,000 stars over two years of observations until 4 July 2020.TESS is currently in its second Extended Mission.All targets are made available to the community as target pixel files (TPFs) and calibrated LCs.TESS light curve files include the timestamps, simple aperture photometry (SAP) fluxes (Twicken et al. 2010;Morris et al. 2020), and pre-search data conditioned simple aperture photometry (PDCSAP) fluxes (Smith et al. 2012;Stumpe et al. 2012Stumpe et al. , 2014)).The SAP flux comes from the sum of the calibrated pixels within the TESS optimal photometric aperture, while the PDCSAP flux corresponds to the SAP flux values corrected for instrumental variations and for crowding.In the process, it usually also removes activity variations with periods similar to the length of a sector.The optimal photometric aperture is defined such that the stellar signal has a high signal-to-noise ratio, with minimal contamination from the background.The TESS detector bandpass spans from ∼ 530 to 1060 nm. Figure 1 shows the TPF files of the field of TOI-238, obtained using the publicly available tpfplotter code (Aller et al. 2020).This code overplots all sources from the Gaia Data Release 3 (DR3) catalogue (Gaia Collaboration et al. 2021) with a magnitude contrast up to ∆m=8 mag on top of the TESS TPFs.There is one faint Gaia sources at the edge of the photometric aperture around TOI-238 automatically selected by the pipeline (∆m=4 mag).Considering the flux-ratio between both, this will have a negligible effect on the planetary radius, assuming the planet transits the brighter star (Ciardi et al. 2015).Mistry et al. (2023), and our own efforts (see Sect.2.2), showed TOI-238.01 transits the central brightest star within aperture, which leads us to consider the extracted TESS light curve to be free of significant contamination from nearby stars.TOI-238 was observed by TESS at the 2-min cadence integrations in sectors 2, 29 and 69.These sectors were processed by the SPOC pipeline (Jenkins et al. 2016) and searched for transiting planet signatures with an adaptive, wavelet-based transit detection algorithm (Jenkins et al. 2010).The SPOC pipeline identified a planet candidate at an orbital period of 1.27241 days2 , and TOI-238 was announced as a TESS object of interest (TOI) via the dedicated MIT TESS data alerts public website3 .The SPOC conducted a transit search of Sector 2 on 04 October 2018 with an adaptive, noise-compensating matched filter (Jenkins 2002;Jenkins et al. 2010Jenkins et al. , 2020)), producing a threshold crossing event for which an initial limb-darkened transit model was fitted (Li et al. 2019) and a suite of diagnostic tests were conducted to help make or break the planetary nature of the signal (Twicken et al. 2018).The TESS Science Office reviewed the vetting information and issued an alert on 29 November 2018 (Guerrero et al. 2021).The signal was repeatedly recovered as additional observations were made in sector 29.The host star is located within 18.07 ± 2.69 arcsec of the source of the transit signal.
We used SAP flux data in our study.The TESS light curve of TOI-238 consists of three sectors of 27.4, 26.2 and 25.8 days, respectively, separated by 706 and 1067 days, respectively.The data has a median precision of 430 parts per million (ppm) with a root mean square (RMS) of 3007 ppm, 3295, and 2881 ppm for sectors 2, 29 and 69, respectively.Figure 2 shows the three TESS light curves.All LCs and TPF files were downloaded from the Mikulski Archive for Space Telescopes4 , which is a NASA funded project.

LCOGT
The TESS pixel scale is ∼ 21 ′′ pixel −1 and photometric apertures typically extend out to roughly 1 ′ , generally causing multiple stars to blend in the TESS photometric aperture.To determine the true source of the TESS detection, we acquired ground-based time-series follow-up photometry of the field around TOI-238 as part of the TESS Follow-up Observing Program (TFOP; Collins 2019) 5 .
We observed a full transit window of TOI-238.01 continuously for 128 minutes in Pan-STARRS z-short band on UTC 2022 November 26 at the Las Cumbres Observatory Global Telescope (LCOGT) (Brown et al. 2013) 1 m network node of the Teide Observatory, Spain.The 4096 × 4096 LCOGT SINISTRO cameras have an image scale of 0 ′′ .389 per pixel, resulting in a 26 ′ × 26 ′ field of view.The images were calibrated by the standard LCOGT BANZAI pipeline (McCully et al. 2018) and differential photometric data were extracted using AstroImageJ (Collins et al. 2017).The TOI-238.01 SPOC pipeline transit depth of 379 ppm is generally too shallow to reliably detect with ground-based observations, so we instead checked for possible nearby eclipsing binaries (NEBs) that could be contaminating the TESS photometric aperture and causing the TESS detection.To account for possible contamination from the wings of neighboring star point spread functions (PSFs), we searched for NEBs out to 2 ′ .5 from TOI-238.If fully blended in the SPOC aperture, a neighboring star that is fainter than the target star by 8.6 magnitudes in TESS-band could produce the SPOC-reported flux deficit at mid-transit (assuming a 100% eclipse).To account for possible TESS magnitude uncertainties and possible deltamagnitude differences between TESS-band and Pan-STARRS zshort band, we included an extra 0.5 magnitudes fainter (down to TESS-band magnitude 19.0).We calculated the RMS of each of the 13 nearby star light curves (binned in 10 minute bins) that meet our search criteria and found that the values are smaller by at least a factor of 3 compared to the required NEB depth in each respective star.Our analysis ruled out an NEB blend as the cause of the detection of the planet candidate TOI-238.01 by the SPOC pipeline.All light curve data are available on the EXOFOP-TESS website 6 5 https://tess.mit.edu/followup 6https://exofop.ipac.caltech.edu/tess/target.php?id= 9006668

ESPRESSO observations
The Échelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO) is a fibre-fed high resolution echelle spectrograph installed at the Very Large Telescope (VLT) telescope array in the ESO Paranal Observatory, Chile (Pepe et al. 2013;González Hernández et al. 2018;Pepe et al. 2021).ESPRESSO has a resolving power of approximately R ∼ 140 000 over a spectral range from ∼380 to ∼788 nm and has been designed to attain a long-term radial velocity precision of 10 cm•s −1 .It is contained in a vacuum vessel to avoid spectral drifts due to temperature and air pressure variations, thus ensuring its stability.Observations were carried out using a Fabry-Perot (FP) etalon as simultaneous calibration.The FP offers the possibility of monitoring the instrumental drift with a precision better than 10 cm•s −1 without the risk of contamination of the stellar spectra by the ThAr saturated lines (Wildi et al. 2010).ESPRESSO can be used on any VLT unit telescope (UT), allowing for an efficient observation and a high-cadence observations.More information can be found on the ESPRESSO user manual. 7SPRESSO is equipped with its own Data Reduction Software (DRS) and provides extracted and wavelength-calibrated spectra, as well as RV measurements.The DRS automatically measures radial-velocity through a Gaussian fit of the cross correlation function (CCF) of the spectrum with a binary mask, computed from a stellar template (Fellgett 1955;Baranne et al. 1996;Pepe et al. 2000).The DRS, version 3.0.0, is available to download from the ESO pipeline website8 .
We obtained 77 individual spectra as part of the ESPRESSO Guaranteed Time Observations (GTO), within programme ID 1102.C-744 (PI: F.Pepe), between May 2019 and December 2021.Measurements were taken in ESPRESSO's 1UT high resolution mode, with fast readout and 1×1 detector binning (HR11).We used 15 minutes of integration time, for a signal to noise ratio (S/N) of 68 in each slice at 550 nm.More information on the different observing modes can be found on the ESO instrument page 9 .
In June 2019, ESPRESSO underwent an intervention to update the fibre link, improving the instrument's efficiency by up to 50% (Pepe et al. 2021) (2022).We used data reduced with DRS 3.0.0and decided not to split the time-series before and after the lamp change.Therefore, we consider a single dataset of 73 ESPRESSO spectra.

HARPS observations
In combination with the ESPRESSO data, we include 50 spectra obtained between May 2019 and July 2019 with the High Accuracy Radial velocity Planet Searcher (HARPS) (Mayor et al. 2003) under the programme ID 1102.C-0249 (PI: D. Arm-strong).HARPS is a fibre-fed high resolution echelle spectrograph installed at the 3.6 m ESO telescope in La Silla Observatory, Chile.It has a resolving power of R ∼ 115 000 over a spectral range from ∼380 to ∼690 nm and has been designed to attain very high long-term radial-velocity precision.It is contained in temperature-and pressure-controlled vacuum vessels to avoid spectral drifts due to temperature and air pressure variations, thus ensuring its stability.HARPS is equipped with its own pipeline providing extracted and wavelength-calibrated spectra, as well as RV measurements and other data products such as cross-correlation functions and their bisector profiles.All observations have been carried out with simultaneous calibration, using the FP.TOI-238 was typically observed once per night using an exposure time of 1800 s, achieving a median S/N of 63 at 550 nm.As in the case of ESPRESSO, the DRS provides RV measurements determined by a Gaussian fit of the CCF of the spectrum with a binary mask computed from a stellar template.

Radial velocities
Using the ESPRESSO and HARPS spectra, we extracted RVs using the S-BART algorithm (Silva et al. 2022a).S-BART builds a high signal-to-noise template by co-adding all the existing observations.Then, unlike usual template-matching algorithms, S-BART uses a single RV shift to describe simultaneously the RV differences between all orders of a given spectrum and the template.The algorithm estimates the posterior distribution of RV shifts after marginalising with respect to a linear model for the continuum levels of the spectra and template, using a Laplace approximation of the posterior probability distribution, an analytical expression of the posterior probability derived from a Gaussian fitting with a mean equal to the maximum a posteriori probability.From the Laplace approximation to the posterior distribution, the mean is used as the estimated RV and the standard deviation as the estimated RV uncertainty for each epoch.With S-BART, and after a sigma-clip based on the measured uncertainty, we measured an RV RMS of 6.2 m•s −1 and 5.55 m•s −1 for HARPS and ESPRESSO data, respectively.We measure median RV internal uncertainties of 1.1 m•s −1 and 0.58 m•s −1 , respectively.The data is spread along a baseline of 952 days.The velocities obtained with S-BART are consistent with those extracted using other template matching software, such as SER-VAL (Zechmeister et al. 2018) or using the CCF technique.For more details on the motivation behind the choice of RV extraction, see Appendix A. Figure 2 shows the RV time-series.We investigated the effect of telluric contamination and the different corrections by comparing RVs coming from uncorrected spectra, RVs from spectra using the automatic correction by S-BART, based on TELFIT (Gullikson et al. 2014), and RVs from spectra corrected following Allart et al. (2022).We opted to use S-BART's automatic correction, which performed very similar to the correction of Allart et al. (2022), and could be applied to all our data in a homogenous way.Appendix B provides with more details.

CCF derived products
The presence of active regions on the stellar disc affects the flux emitted by the star and its velocity field, distorting the shape of the lines and the PSF measured by the instrument.This effect manifests itself as changes in the width, depth, and symmetry of the CCF.These changes are usually related to the changes in flux, or their gradient, caused by the active regions; and their scale is related to the coverage of active regions, their contrast, and v sin i of the star.
To measure these variations, we use the full width at half maximum (FWHM), the bisector span (Queloz et al. 2001) and the contrast of the CCF.All those quantities are automatically provided by the ESPRESSO and HARPS DRS.Over our complete time-series, the CCF FWHM shows a dispersion of 17.0 m•s −1 , with a median uncertainty of 1.9 m•s −1 .The dispersion of the bisector span time-series is 6.2 m•s −1 with a median uncertainty of 1.9 m•s −1 .The contrast of the CCF has a dispersion of 0.249 with a median uncertainty of 0.011.Figure 2 shows the FWHM and bisector span time-series.Figure D.1 shows the CCF contrast time-series.

Telemetry data
Modern RV spectrographs are designed to minimise instrumental effects caused by the changes in their environments.However, small effects either linked to the stability of both instruments or to the extraction of the velocities can be present in the data.This potentially biases some of the results (Suárez Mascareño et al. 2023).We collected the time-series of the temperature of the Échelle gratings of ESPRESSO and HARPS, and the Barycentric Earth Radial Velocity (BERV) as proxies for variations potentially induced by instrumental changes or imperfect data extraction.We measure an RMS of the variations of the temperature of the Échelle gratings of ESPRESSO and HARPS of 11 mK and 5 mK, respectively.See Appendix C for more details.

Spectral activity indicators
In addition to the CCF indicators, we studied the time-series of several standard chromospheric indicators.

Ca II H&K
The intensity of the emission of the cores of the Ca II H&K lines is linked to the strength of the magnetic field of the star, which in turn is well correlated with the stellar rotation period of FGKM stars.The measured emission intensity of the line cores also changes when active regions move across the stellar disc, helping us trace the rotation of the star.We calculate the Mount Wilson S -index for the ESPRESSO data following Vaughan et al. (1978).We define two triangular-shaped passbands with a FWHM of 1.09 Å centred at 3968.470 Å and at 3933.664 Å for the Ca II H&K line cores.For the continuum we use two bands 20 Å in width centred at 3901.070 Å (V) and 4001.070Å (R).Then the S-index is defined as: where ÑH , ÑK , ÑR , and ÑV are the mean fluxes per wavelength unit in each passband, while α and β are calibration constants fixed at α = 1.111 and β = 0.0153 .The S index (S MW ) serves as a measurement of the Ca II H&K core flux normalised to the neighbour continuum.We measure a median S-index of 0.320 ± 0.064, which using the calibration of Suárez Mascareño et al. (2015) translates to a log 10 (R ′ HK ) of -4.74 ± 0.27.Following the activity-rotation calibration for K-dwarfs of Suárez Mascareño et al. (2016), we estimate a rotation period of 28 ± 16 days, with the uncertain-ties being a result of the large uncertainties in the photometric magnitudes shown in Table 1.

Hα
Similar to Ca II H&K, the emission in the core of the Hα line (or filling, in the case of low activity stars) is related to the strength of the magnetic field, and the presence of active regions on the stellar disc.We can also use it to track the motion of said regions across the stellar disc, and thus the stellar rotation.We compute the Hα values following the definition of Gomes da Silva et al. (2011). Figure D.3 shows the Hα index time-series.We measure a median Hα index value of 0.02845 ± 0.00012.

Na I
The sodium D lines have been shown to be good chromospheric indicators for low-mass stars.We compute the time-series of the Na I D fluxes following the definition of Díaz et al. (2007).Figure D.4 shows the Na I index time-series.We measure a median Na I index value of 0.00534 ± 0.00045.

Stellar parameters of TOI-238
We obtained stellar atmospheric parameters from the highresolution ESPRESSO spectra via the ARES+MOOG method (Sousa 2014).The ARES+MOOG method is a curve-of-growth technique based on neutral and ionised iron lines.Equivalent widths of these spectral lines were measured automatically from the stacked ESPRESSO spectrum using ARESv2 (Sousa et al. 2015).We determine effective temperature (T eff ), surface gravity (log g), iron abundance ([Fe/H]), and microturbulent velocity by imposing excitation and ionisation equilibria.For this purpose, we used the radiative transfer code MOOG (Sneden 1973), that assumes local thermodynamic equilibrium (LTE) employing a grid of ATLAS plane-parallel model atmospheres (Kurucz 1993).Subsequently, we corrected the surface gravity for accuracy (Mortier et al. 2014) and added systematic errors in quadrature to our precision errors for the effective temperature, surface gravity, and iron abundance (Sousa et al. 2011).To compute the stellar mass and radius, we used PARAM (da Silva et al. 2006;Rodrigues et al. 2017).This code matches the stellar parameters T eff , log g, and [Fe/H] obtained in the previous analysis, along with the Gaia DR3 parallax and the V magnitude published in the literature, to a grid of stellar evolutionary tracks and isochrones from PARSEC (Bressan et al. 2012).The iron-to-silicate mass fraction ( f star iron ) of the planetbuilding blocks, as estimated from their host star abundances, can be used as a proxy for the composition of small planets in the planetary system (Adibekyan et al. 2021).We derived the abundances of Mg and Si using the same tools and models as for stellar parameter determination as well as using the classical curve-of-growth analysis method assuming local thermodynamic equilibrium.Although the EWs of the spectral lines were automatically measured with ARES, for Mg which has only three lines available we performed careful visual inspection of the EWs measurements.For the derivation of abundances we closely followed the methods described in Adibekyan et al. (2012Adibekyan et al. ( , 2015)).Table 1 shows the full list of parameters.
We crosschecked our results using the STEPAR code (Tabernero et al. 2019), and by computing the parameters from the spectral energy distribution (SED) using the available photometry.STEPAR provided us with consistent stellar parameters: T eff = 5054 ± 78 K, log g = 4.52 ± 0.20, [Fe/H] = -0.06± 0.04, -4.74 ± 0.27 0 P rot GP [days] 26.0 ± 1.  1.It is important to note that the uncertainties of derived stellar masses and radii are typically smaller than the uncertainties of direct measurements.As these parameters can affect the computation of planetary parameters, we enlarged the error bars of the stellar parameters in all our calculations following the reccomendations of Tayar et al. (2022).Using Gaia DR3 positions, proper movements, RV and parallax, we computed the galactic velocities of TOI-238 following Johnson & Soderblom (1987).We obtain galactic velocities U = 14.13 ± 0.09 km•s −1 , V = 2.83 ± 0.10 km•s −1 , W = 19.50 ± 0.26 km•s −1 .These values are fully compatible with the thin disk of the galaxy and do not match any known young moving group (Gagné et al. 2018).This result indicates the star to be older than 0.8 Gyr, which is consistent with the ages derived from the spectroscopic modelling.

Telemetry and stellar activity
We started by analysing the possible effect of the changes in environmental conditions and imperfect data extraction in the science data.We studied the correlations between the BERV and temperature of the Échelle gratings of ESPRESSO and HARPS, and the different science data.Similar to Suárez Mascareño et al.
(2023), we did not detect any significant correlation between any of them and the radial velocities.In some of the CCF and spectroscopic indicators, we found correlation between the temperature of the optical elements and the BERV.We interpret this as an effect of a subtle change in focus (correlation with Éch.Temp.) and leftover telluric contamination (correlation with BERV).Appendices B and C provide a more detailed analysis.
We then analysed all stellar indicators described in Sect. 2 to search for periodicities related to stellar activity and that could create false-positive detections in RV or bias the amplitude measurements in RV.We modelled each time-series independently, using Gaussian Processes regression (GP; see Rasmussen &Williams 2006 andRoberts et al. 2012).These results are later used to decide the characteristics of the global model.
The GP framework has become one of the most successful methods in the analysis of stellar activity in RV time-series (e.g.Haywood et al. 2014).The stellar noise is described by a covariance with a prescribed functional form and the parameters attempt to describe the physical phenomena to be modelled.The GP framework can be used to characterise the activity signal without requiring a detailed knowledge of the distribution of active regions on the stellar surface, their lifetime, or their temperature contrast.One of the biggest advantage of GPs is that they are flexible enough to effortlessly model quasi-periodic signals, and account for changes in the amplitude, phase, or even small changes in the period of the signal.This flexibility is also one of their biggest drawbacks, as they can easily overfit the data, suppressing potential planetary signals.
We used the newly presented S+LEAF code (Delisle et al. 2022)10 , which extends the formalism of semi-separable matrices introduced with celerite (Foreman-Mackey et al. 2017) to allow for fast evaluation of GP models even in the case of large datasets.The S+LEAF code supports a wide variety of GP kernels with fairly different properties.We opted for the ESP Kernel, which is a very close approximation of the widely used Quasi-Periodic Kernel: where A is the variance of the data, P rot is the periodicity of the data, L is the evolutionary timescale and ω is the length scale of the periodic component.The covariance matrix also includes a term of uncorrelated noise (σ), independent for every instrument.This term is added in quadrature to the diagonal of the covariance matrix to account for all unmodelled noise components, such as uncorrected activity or instrumental instabilities.
The ESP Kernel in the S+LEAF package allows for a explicit configuration of the number of harmonics of the main variability to include.We used the default configuration, which includes the first two harmonics.
Considering the correlations with the telemetry mentioned before, we included low-order polynomials with respect to the temperature of the Échelle grating and with respect to the BERV, when we deemed appropriate.We found that none of them improved on the GP-only models.All the polynomial parameters where consistent with zero within 1σ.Contrary to what was seen in Suárez Mascareño et al. (2023), it seems that any small variation induced by instrumental changes is too subtle to be confidently separated from the activity model.TOI-238 is more active than GJ 1002, which means the ratio between activity and instrumental variations will be much larger.The analysis of some indicators also included a low-order polynomial or a sinusoidal to account for long-term variations.None of these models improved on the GP-only model.
We opted to restrict the analysis to a GP-only model, with floating offsets and jitter terms.We modelled the data using Bayesian inference through the nested sampling (Skilling 2004(Skilling , 2006) ) code Dynesty (Speagle 2020;Koposov et al. 2023).All models had seven free parameters (N f ree ).We went for a number of live points of 10 × N f ree .We used the default bounding option, multi-ellipsoidal decomposition (Feroz et al. 2009), and chose to sample our parameter space using random walk (Skilling 2006), which is very efficient at exploring low-to-mid dimension parameter spaces.
Additionally, for every indicator we study the relationship between the indicator and the derivative of its best fit GP model, with the RV measurements.We measured the slope between both quantities by fitting a linear function, using the measurement uncertainties as weights.We then used the diagonal of the covariance matrix to evaluate the significance of the slope.We consider that there might be a relationship when the slope is ≥2σ different from zero.This information can later be leveraged when building a multi-dimensional GP model.
For all models, we use log-normal priors for the amplitude of the GP, with a mean corresponding to ln(RMS) of the data, and a standard deviation of ln(RMS) of the data.We use the same priors for the jitter values.We explore rotation periods between 5-50 days, and evolutionary timescales of 10-200 days, both in log-uniform scale.Main-sequence K-dwarfs are not expected to have rotation periods longer than 50 days, and the timescale of coherence of the rotation signal is typically of 2-3 rotations (Giles et al. 2017).For the zero-points of the time-series we use a normal prior with a mean of zero and a standard deviation of the RMS of the data.
We obtained good fits for the FWHM and BIS of the CCF, with significant correlations between the activity proxies and their derivatives, and the RV data.We find that correlations between FWHM and BIS, and the RV, show opposite slopes.We also obtained a good fit for the Hα time-series, very similar to that of the FWHM and showing a very similar correlation with the RV as the FWHM.We opted for using the FWHM and the BIS of the CCF in our global model, as both provided different information about the stellar activity variations.
Figure 3 shows the data, model, periodogram and correlation plots of the FWHM of the CCF.The time-series of the FWHM of the CCF show an RMS of 17.0 m•s −1 , with the data having a median uncertainty of 2.1 m•s −1 .For this time-series, we obtain a period of 24.3 +1.4  −1.1 days, with a damping timescale of 22.7 +5.7 −4.7 days.The period measured with the GP is very close to the period detected by the GLS periodogram (Zechmeister & Kürster 2009).We obtain a good fit to the data, with reasonably clean variations during all observing campaigns.After subtracting the model, we measure an RMS of the residuals of 6.7 m•s −1 , and find no significant signals in their periodogram.We report moderate positive relationships between the FWHM and the RV measurements, as well as the gradient of the best-fit model and the RV measurements.
Figure 4 shows the data, model, periodogram and correlation plots of the bisector span of the CCF (Queloz et al. 2001).We measure an RMS of the data of 6.4 m•s −1 , and a median uncertainty of 1.9 m•s −1 .We obtain a period of 25.28 +0.82  −0.65 days, with a damping timescale of 52 +24 −18 days.The period measured with the GP is close to the period detected by the GLS periodogram and to the period measured in the FWHM time-series.The GLS of the data shows a structure with a periodicity of around 1 year.However, the inclusion of additional parameters to account for it is disfavoured in the fitting.After subtracting the model, we measure an RMS of the residuals of 3.7 m•s −1 , and find no significant signals in their periodogram.We find moderate anti-correlations between the bisector span and the RV measurements, and a moderate positive correlation between the gradient of the best-fit model and the RV measurements.This indicates that the RV data is significantly affected by stellar activity variations.
The models for the additional activity proxies can be found in Appendix D.

Global model -Activity only
We built the global model combining TESS photometry and RV, and FWHM and BIS as activity proxies.We modelled the stellar variations using a GP with S+LEAF, again using the ESP Kernel.The S+LEAF code allows to fit simultaneously a GP to several time-series, based on a linear combination of the GP and its derivative, with different amplitudes for each time-series.
The usual assumption is that there exists an underlying function governing the behaviour of the stellar activity, which we denote G(t).G(t) manifests in each time series as a linear combination of itself and its gradient, G ′ (t), with a set of amplitudes for each time series, following the idea of the FF ′ formalism (Aigrain et al. 2012;Rajpaul et al. 2015), as described in equation 3.
with ∆ T S i representing the variations of each time-series, and A i and B i the scaling coefficients of the underlying funciton G(t) and its gradient, G ′ (t), respectively.We modelled the data using Bayesian inference via nested sampling with Dynesty, with the same configuration described before, but with one exception.Now, we sampled the parameter space using random slice sampling, which is well suited for the high-dimensional spaces (Handley et al. 2015a,b) resultant of modelling several time-series at once.We used a number of slices equal to 2 × N f ree .
We fitted the RV combined with the FWHM and/or the BIS following this idea, with moderate success.Depending on the specific model, we derived different timescales, and the residuals after the fit contained periodic signals that we could still attribute to stellar activity.For example, the combined fit of the RV and FWHM left some periodic signals present in the BIS data and vice versa.The combined model of the three time-series resulted in moderate under-fit that left a few periodic signals in the RV which also appeared in the FWHM or BIS data.
Following the conclusions of Tran et al. (2023), we attempted a model with two latent GPs.Tran et al. (2023) demonstrated that in some circumstances a single multi-dimensional GP is not capable of following all activity timescales simultaneously.As the variations in FWHM and BIS appear to show different timescales (see Sect 4.1), and their effect in the RV seems to be different, we built a model that linked the FWHM to the RV, the BIS to the RV, but not the BIS to the FWHM.
We included the TESS data in the global model with its own independent activity model, that is not linked to the RV or to any of the activity proxies.TESS SAP photometry is often contaminated by instrumental trends, and the rotation period of this star is too long to be preserved in the PDCSAP photometry.The variations seen in the SAP photometry are likely a combination of stellar and instrumental variability.We modelled TESS variability with a GP with a simple harmonic oscillator (SHO) to take advantage of its scalability to datasets with thousands of points.We resampled the TESS data into 10 minutes bins to avoid an excessive increase of computing time.The cost of the evaluation of the GP model scales linearly with the number of points.Modelling the data in 2-minute bins takes five times longer than in 10-minute bins, while it does not offer any advantage given the expected characteristics of the transit.
where η = (1 − (2L/P i ) −2 ) 1/2 controls the damping of the oscillator.This Kernel has a power spectrum density: where ω is the angular frequency, ω i is the undamped angular frequency for each component (ω i = 2 π / P i ), S i is the power at ω = ω i for each component, and Q i is the quality factor.S i , P i and Q i are the parameters sampled in the covariance matrix, which are related to the amplitude (C i ), rotation period (P rot ) and timescale of evolution (L) in the following way: Our global activity model, including four different timeseries and three GP kernels, is then defined as follows: with G TESS being a SHO Kernel (following eq.4), and G FWHM and G BIS being ESP Kernels (following eq. 2), with variance equal to unity.The model also includes a floating zero-point for every time-series and instrument, and jitter terms for the RV, FWHM and BIS data.In the case of TESS data, we fix the jitter by estimating the RMS of the data in the flat region between the start of the observations and BJD 2458356 (see panel a of Fig. 2).
For the model of the TESS data we use the following parameters and priors: A TESS is the amplitude of the TESS variability and uses a log-uniform (LU) prior with a range of [0 , 20].P TESS is the period of the variability of the TESS data and uses a normal (N) prior around 28 days (highest peak in the periodogram) with a standard deviation of 10%.L TESS is the evolutionary timescale of the TESS variability.Based on Giles et al. (2017) we use a N prior 4.0 ± 0.3 (54 +20 −14 days).The zero-point of the TESS data uses a N prior centred at 0, with a standard deviation equal to the RMS of the data.
The combined FWHM+BIS+RV model uses the following parameters and priors: A FWHM and A BIS are the amplitudes of the variability in the FWHM and BIS time-series and use uniform (U) priors between 0 and 5× the RMS of their respective time-series.P FWHM and P BIS are the periods of variability of the FWHM and BIS, and use N priors of 25 ± 2.5 days (∼ the weighted mean of all periods measured in the activity indicators).L FWHM and L BIS are the evolutionary timescales of the FWHM and BIS variability and, once again based on Giles et al. (2017), use N priors 4.0 ± 0.3 (54 +20 −14 days, consistent although slightly longer than the values measured in the activity indicators).These priors for the L parameters are 1σ compatibles with the values obtained for the L parameters in the independent FWHM and BIS models.ω FWHM and ω FWHM are the length scale of the periodic components of the variability, and use LU priors with a range of [-5 , 5].A11 RV and A21 RV are the amplitudes of the variability in RV related to the FWHM and BIS variations, respectively, and use U priors [-5× RMS , 5× RMS].A12 RV and A22 RV are the amplitudes of the variability in RV related to the gradient of the FWHM and BIS variations, respectively, and use U priors [-10× RMS , 10× RMS].The priors in the RV amplitude allow for negative values to take into account possible negative correlations between RV and FWHM.The zero points of all datasets use N priors centred at 0, with sigmas equal to the RMS of their respective dataset.The white noise component of all datasets uses LN priors centred around the log(RMS) of their respective datasets, with a sigma of log(RMS).
Figures 5 and 6 show the best-fit model, using the parameters obtained from the posterior distribution of the nested sampling.In the case of the TESS data, we do not gain a lot of information about the stellar activity of the star.However, the GPmodel manages to perform a smooth fit that we use to detrend the data.We apply a Transit Least Squares (TLS, Hippke & Heller 2019) periodogram to the residuals under the fit.The periodogram shows its highest peak at 1.27 days, the period of the candidate TOI 238.01.The periodogram also shows peaks at ∼ 2.54 days (2× 1.27 days) and 8.47 days.
The model on the FWHM+BIS+RV provides significantly different timescales of variability for the FWHM+RV and BIS+RV components.The FWHM+RV component shows a periodicity of 24.49 ± 0.86 days, with a timescale of evolution of 35.5 +6.7 −5.8 days, while the BIS+RV shows a periodicity of 26.53 ± 0.46 days, with a timescale of evolution of 47.8 +10.1 −9.6 days.The omega parameters are 0.76 +0.36  −0.25 , and 0.74 (3σ), respectively.The variability shows a significant RV amplitude related to the FWHM variations (A11 RV = 4.7 +2.1 −1.5 m•s −1 ), a non significant RV amplitude correlated with the gradient of the FWHM variations (A12 RV = -1.6 +4.6  −5.0 m•s −1 ), a significant RV amplitude anti-correlated with the BIS variations (A21 RV = -5.1 +1.2 −1.4 m•s −1 ) and an almost-significant RV amplitude correlated with the gradient of the BIS variations (A22 RV = 8.6 +3.9 −4.0 m•s −1 ).These results show once again that both indicators are affected by different types of variability, and they are both present in the RV data.We apply the best model to the data and subtract it to study the residuals.
The RMS of the residuals after the fit of the RV data goes down to 3.0 m•s −1 from the original 5.86 m•s −1 (49% reduction), with the ESPRESSO and HARPS data retaining a very similar RMS.The RMS of residuals after the fit to the FWHM, and BIS, data is of 6.9 m•s −1 and 4.2 m•s −1 , respectively.These values correspond to a reduction of the variability of 59 % for the FWHM data and 32% for the BIS data.
The GLS periodogram of the residuals has its most dominant peak at 1.27 days (see panel f of Fig. 6) again the period of the candidate TOI 238.01.Some other peaks appear in the periodogram.Most of them are related to the aliases of the 1.27 days periodicity and the gaps inbetween the observations.It is worth noting that the periodogram of the raw velocities shows several significant peaks.A forest of peaks around 25 days (stellar rotation), another at 14 days, and a third one at 8.47 days (see panels c and f of Fig. 6).The latter has a period that matches one of the peaks in the TLS periodogram of the TESS detrended data (Hippke & Heller 2019).It is possible that all of them are artefacts of the rotation signal.Yet, it may be that the GP is incorrectly suppressing some of them due to their closeness to the rotation and its harmonics.

Global model -One planet
Following the results obtained after detrending the data with an activity-only model, we included the possibility of one planet in the TESS and RV data.
To model the transit, we used the pytransit package (Parviainen 2015) with quadratic limb darkening (Mandel & Agol 2002).For the orbital period, we used a N prior of 1.27 ± 0.1 days.To avoid a multi-modal posterior for the time of transit (T0), we parameterised it as a phase centred around the expected last transit of the TESS light curve.The time of transit is defined as: where x[-1] T ES S is the last TESS observation, P b is the orbital period of the planet, and Ph b is the phase, using a U prior [0,1].We sampled the planet radius directly, with LU prior [-5,5] in R ⊕ (corresponding to ∼ 0 to 148 R ⊕ ).The LU prior allows for better sensitivity when sampling very small values (e.g. in the case of having no transit signal).For the impact parameter, we used a U [0,1] prior.We parameterised the eccentricity as e = ( √ e cos(ω)) 2 + ( √ e sin(ω)) 2 and ω = arctan 2( √ e sin(ω), √ e cos(ω)).We then sample √ e cos(ω) and √ e sin(ω) with normal priors of 0 ± 0.3.Eccentricity is typically overestimated in noisy data and datasets with unmodelled sources of variability (Hara et al. 2019).The parameterisation described above favours low eccentricities, which are expected for close-in planets, while still allowing for eccentricities up to 1.To estimate the limb darkening priors we used the LDTK package (Parviainen & Aigrain 2015), which provides limb darkening coefficients, and uncertainties, for a given set of the stellar parameters and a given observing passband.We used N priors centred in the results given by LDTK, with sigmas of 3× the provided uncertainty.Even with simple 2 parameter models, such as the quadratic limb darkening, there can be degeneracies between them.Several combinations can yield indistinguishable results, in particular at low signal to noise.LDTK provides a likelihood term that estimates the log-likelihood of the combination of both parameters.We included it as an additional term in our likelihood evaluation.
For the planetary component in the RV data, we used a Keplerian model: where the true anomaly η is related to the solution of the Kepler equation that depends on the orbital period of the planet P orb and the orbital phase ϕ.This phase corresponds to the periastron time, which depends on the mid-point transit time T 0 , the argument of periastron ω, and the eccentricity of the orbit e.
The P orb , mid-point transit time T 0 argument of periastron ω, and eccentricity, are shared with the transit model and their parameterisation is described above.Instead of the amplitude (K), we sampled the planet mass directly, with LU prior [-5,5] in M ⊕ (corresponding to ∼ 0 to 148 M ⊕ ).
We included the stellar radius and mass in the model as they are needed for the computation of several transit parameters, and the transformations between planet radius and mass, and their respective amplitudes.We sampled them using N priors, centred around the values shown in Table 1.The standard deviation of the distributions is the uncertainty derived in the stellar parameters, and enlarged following the results of Tayar et al. (2022), to account for the systematic uncertainties of the model, which result in R s = 0.733 ± 0.034 R ⊙ and M s = 0.790 ± 0.043 M ⊙ .This way, we ensure proper propagation of uncertainties to the rest of the parameters sampled in the model.
We obtain a significant detection of the candidate TOI 238.01, both in the TESS light curve and the RV data.We measure an orbital period P b of 1.2730991 ± 0.0000029 days and a time of mid-transit at BJD 2460207.0240± 0.0033 days.We measure a radius R b of 1.387 ± 0.086 R ⊕ , and a mass m b of 3.90 ± 0.51 M ⊕ .These values correspond to a planet-star contrast (R planet /R star ) of 0.01704 ± 0.00084 and a RV semi-amplitude of 2.72 ± 0.32 m•s −1 .We measure an impact parameter < 0.61 (σ).We measure a ∆ lnZ of +495.7 with respect to the activityonly model, which extremely favours the model with a planet with respect to the activity-only model.The RMS of the residuals after the fit of the RV data goes down to 1.53 m•s −1 , from the original 5.86 m•s −1 (74% reduction), with the ESPRESSO data showing an RMS of the residuals of 0.93 m•s −1 and the HARPS data showing 2.11 m•s −1 .
Figures 7 and 8 show the best-fit model, using the parameters obtained from the posterior distribution of the nested sampling.Table E.1 shows the priors and measured values for all parameters.The TLS periodogram of the residuals of the TESS light curve shows a dominant peak at 8.47 days.The same peak appears in the GLS periodogram of the activity-induced signal in RV.This suggests that the GP is potentially absorbing an additional planetary signal in the RV data.

Global model -Two planets
We investigated the 8.47 days peak in the TLS periodogram by including a second Keplerian component in the model.We sample the period using a N prior of 8.5 ± 0.5 days.The rest of the components of the planet model are the same described in section 4.3.
We obtain a significant detection of a planetary signal in both the TESS light curve and the RV data.We measure an orbital period P c of 8.465651 ± 0.000035 days and a time of mid-transit at BJD 2460204.9675± 0.0050 days, a radius R c of 2.18 ± 0.18 R ⊕ , and a mass M c of 6.7 ± 1.1 M ⊕ .These values correspond to a R planet /R star of 0.0267 ± 0.0021 and a RV semi-amplitude of 2.46 ± 0.41 m•s −1 .We measure an impact parameter of 0.774 ± 0.057.This model yields a ∆ lnZ of +30.4 with respect to the 1-planet model, and +528.1 with respect to the activity-only model.These values paint the 2-planet model as much more likely than either the 1-planet model or the activity-only model.The RMS of the residuals after the fit of the RV data goes down to 1.61 m•s −1 , from the original 5.86 m•s −1 (72% reduction), with the ESPRESSO data showing an RMS of the residuals of 1.50 m•s −1 and the HARPS data showing 1.75 m•s −1 .The RMS of the residuals is slightly higher than with the 1-planet model, which might be explained by the GP overfitting the data in the 1-planet model.Additionally, we repeated the modelling process in the RV and TESS data independently.We found consistent results for the planetary parameters in our combined analysis, a RV-only analysis and a TESS-only analysis.Figure 11 shows the comparison between the periods and times of mid transits obtained in these three models.We tested the effect of different GP Kernels (such as S+LEAF's MEP Kernel or different combinations of SHO Kernels).We found no significant differences in the planetary parameters, however the activity model would be sometimes very different.We found the ESP Kernel to be the most reliable at producing a smooth model that did not show obvious signs of overfitting.
After fitting a model with a GP and two planetary signals, the periodograms of the residuals of the TESS light curve and the RV don't show any clear peak that could be attributed to a planetary signal.There is, however, a peak at 14 days in the GLS of the activity-induced component of the RV that does not fully fit as one of the harmonics of the rotation period.

Additional planets
Although we find no more prominent periodicities in the TLS periodogram of the residuals of the TESS light curve, or the GLS periodogram of the residuals of the RV, we test for the presence of additional planets that might have been missed.After fitting the model that includes a GP and two planetary signals, the RMS of the residuals of the RV data is 1.61 m•s −1 , with the ESPRESSO data showing an RMS of 1.50 m•s −1 and the HARPS data 1.75 m•s −1 .These values are significantly larger than the internal errors of the RV of 0.58 m•s −1 and 1.10 m•s −1 for ESPRESSO and HARPS, respectively.This large difference can be due to a combination of systematic issues, unmodelled activity variations and/or additional planets in the system.Considering the impact parameter of planet c, we do not expect any additional detectable transiting planet in the data.However, it is possible that the signal of a non-transiting planet has been absorbed by the GP model.The GLS periodogram of the suspected activity-induced RV data showed a peak at 14 days that does not fully fit as one of the harmonics of the measured rotation period (see figure 10).To investigate this signal, and others that might have gone unnoticed, we add an additional sinusoidal component in the RV data.We write the RV variation as: which ensures that T 0 coincides with the time of mid transit.We test two different models: a blind search, in which we define the prior for the orbital period as LU between 10 and 100 days, and a guided fit in which we define the prior for the orbital period as N around 14 ± 2 days.The T 0 and mass are defined as described in section 4.3.
In the case of the blind search, we obtain a power excess at ∼ 30 days orbital period, with a non-significant mass measurement and a wide range of possible values for the phase.While this periodicity is not consistent with the variability measured by the combined GPs, it is consistent with a secondary peak in the periodogram of the FWHM (see Fig. 3).Althought this might not be enough information to fully establish the origin of the signal, it is sufficient to not consider a planetary origin.
In the guided search, the model converges to a period of 13.992 ± 0.033 days and a possible mass of 8.4 +2.2 −3.0 M ⊕ (2.8 σ).We measure a ∆ LogZ of ∼ 3.These results hint at the presence of an additional signal, independent from the activity model and the other two planets, but are not decisive enough to claim a detection.
We assesed the significance of the detection of those two signals using the False inclusion probability (FIP) framework (Hara et al. 2022).This framework uses the posterior distribution of the nested sampling run, and compute the probability of having a planet within a certain orbital frequency interval based on the fraction of samples within that frequency interval.Figure 12 shows the FIP periodogram of the two 3-planet models tested.We identify clear significant peaks for the signals at 1.27 days and 8.47 days.However, the peaks at 14 days and 30 days appear well below the 1% threshold.
With the data at hand, and within the framework of our analysis, there is no clear evidence of the presence of additional planets in the system.There might be an additional non-transiting planet with a period of 14 days, and a minimum mass of ∼ 8.4 M ⊕ , but with the current dataset and analysis technique, we cannot confidently disentangle the signal from the activity-induced signal.

The planetary system of TOI-238
We detect the presence of two planetary signals both in the TESS and RV data.TOI-238 b is a planet with an orbital period of 1.2730991 ± 0.0000029 days and a time of mid-transit at BJD 2460207.0240± 0.0033 days.It has a radius of 1.402 ± 0.086 R ⊕ (R planet /R star = 0.01721 ± 0.00083), and a mass of 3.40 ± 0.46 M ⊕ (K = 2.36 ± 0.32 m•s −1 ).It orbits at a distance of 0.02118 ± 0.00038 au from its parent star.The planet's orbit is consistent with circular, with an upper limit to the eccentricity of e = 0.18 to within 3σ.Assuming a bond albedo of 0.3, we estimate an equilibrium temperature of equilibrium (T eq ) of 1311 ± 28 K. TOI-238 c is a planet with an orbital period of 8.465651 ± 0.000035 days and a time of mid-transit at BJD 2460204.9675± 0.0050 days.It has a radius of 2.18 ± 0.18 R ⊕ (R planet /R star = 0.0267 ± 0.0021), and a mass of 6.7 ± 1.1 M ⊕ (K = 2.46 ± 0.41 m•s −1 ).The planet's orbit is consistent with circular, with an upper limit to the eccentricity of e = 0.34 to within 3σ.It orbits at a distance of 0.0749 ± 0.0013 au.Assuming bond albedos of 0.3 we estimate an equilibrium temperature of 696 ± 15 K. Table 2 shows all the parameters of the planets in the system.Figures 13 and 14 show the phase-folded plots of the TESS light curve and the RV time-series, at the periods of both planets b and c, after subtracting the activity model and, in each case, the other signal.
Figure 15 shows the mass-radius diagram of known exoplanets, with mass and radius measured with uncertainties smaller than 33%, in the range between 1 and 20 M ⊕ and 0.7 and 4.2 R ⊕ .TOI-238 b is a super-Earth whose mass-radius configuration is consistent with an Earth-like composition.TOI-238 c on the other hand is consistent with a water world (Léger et al. 2004) with ∼ 50% water or a rocky core with a very small hydrogen envelope.The planet falls exactly at the intersection of these two tracks, making it impossible to distinguish the correct composition based on the mass and radius alone.The large uncertainty in the radius (∼ 9%) prevents us from fully rejecting both Earth-like compositions or a water-rich planet with small Hydrogen envelope (both within 3σ).A better determination of the radius, e.g. using high-precision photometry of missions such as the CHaracterising ExOPlanet Satellite (CHEOPS, Benz et al. 2021), could constrain better the nature of planet c, however it is unlikely that it would be sufficient to uniquely establish it.Using the parameters described above, we estimate densities for planets b, and c, of 6.8 ± 1.6 g•cm −3 , and 1.83 ± 0.66 g•cm −3 , respectively.Planet b has a density larger than Earth, expected for super-Earths with Earth-like composition.Planet c has a density similar to Neptune.
One of the goals of the ESPRESSO-GTO is probing the thresholds leading to significant evaporation of planetary envelopes.For this purpose, the project investigated systems in which small rocky planets and sub-Neptunes coexisted under intense irradiation.TOI-238 was originally not part of this sample, 728 ± 78 T eq A=0.3 [K] 1311 10204.9675 ± 0.0050 P orb [d] 8.465651 ± 0.000035 as there was a single identified planet candidate, but the discov- ery of the second planet grants it a place inside.To study the dependence of the density of these planets with the respective flux received from their host star, we calculated the incident flux from the parent star (S) as a function of the effective temperature, the radius of the star, and the semi-major axis.We estimate a S = 728 ± 78 S ⊕ for TOI-238 b, and S = 58.3± 6.2 S ⊕ for Figure 16 shows the histogram of the distribution of radii of well characterized transiting planets.The distribution is bimodal, showing the known radius valley (Fulton et al. 2017).TOI-238 b falls right in the peak of low-radius planets, while TOI-238 c is at the inner edge of the large-radius mode.The figure also shows the histograms of the densities, normalized to the density of the Earth-like composition track, similar to what Luque & Pallé (2022) showed.Once again we find a bimodal distribution, separating high-density and low-density planets, with a somewhat deeper valley than what is seen in radius.The planets of TOI-238 fall one on each side.If we restrict the histogram to the densities of planets orbiting K-dwarfs, there is a gap at ρ/ρ ⊕,s ∼ 0.5, separating rocky planets from puffier planets.However, the amount of data is low enough for it to be a statistical anomaly.In both samples, it is difficult to identify a frontier between different types of low-density planets.Figure 17 shows density of planets included in Figure 15, normalized to the density of the mass-radius track of Earth-like composition.It is possible to visually identify the two populations hinted in the density histogram (highlighted as shaded regions).
There are currently several hypotheses behind the origin of these distinct populations of low-mass exoplanets, such as photo-evaporation (Owen & Wu 2017), core-powered mass loss (Ginzburg et al. 2018) or migration after formation beyond the ice-line (Mordasini et al. 2009).Following Owen & Campos Estrada (2020) as described in Cloutier et al. (2020) we estimate that the system is consistent with photevaporation as long as planet c is larger than 1.04 ± 0.11 M ⊕ and planet b is smaller than 12.2 ± 1.7 M ⊕ .In accordance with Cloutier et al. (2020) we estimate that the system is consistent with core-powered mass loss as long as planet c is heavier than ∼ 2.78 M ⊕ and planet b is lighter than ∼ 8.05 M ⊕ .With their current mass measurements, the planetary system of TOI-238 is comfortably consistent with both photo-evaporation and core-powered mass loss scenarios.Additionally, the position of both planets within the mass-radius diagram is consistent with the observation of Luque & Pallé (2022) that in the case of multi-planetary systems including water-worlds, the inner planet is usually a less massive rocky planet, which would be in agreement with the predictions from type I migration models (Mordasini et al. 2009).
We modelled the interior structure of the planets using the machine learning inference model ExoMDN (Baumeister & Tosi 2023).ExoMDN provides full inference for the interior structure of low-mass exoplanets, based on synthetic data created with the TATOOINE code (Baumeister et al. 2020;MacKenzie et al. 2023), using the mass, radius and equilibrium temperature of the planets.The determination of interior structures, based only on those parameters, is a degenerate problem.A potential path to overcome this issue is the use of the elemental abundances of the host star, which may be representative of the bulk abundances of the planet and its atmosphere (Dorn et al. 2015).This approach requires some assumptions about the history of formation and evolution of the system.Another possible approach is through the measurement of the fluid Love numbers.Fluid Love numbers describe the shape of a rotating planet in hydrostatic equilibrium.The second-degree Love number k2 depends on the interior density distribution (Kellermann et al. 2018).In a body with k2 = 0, the entire mass is concentrated in the center, while k2 = 1.5 corresponds to a fully homogeneous body.For a number of exoplanets, the fluid numbers are potentially measurable through second-order effects on the shape of the transit light curve (Akinsanmi et al. 2023).ExoMDN uses the k2 number to break some of the degeneracy between interior configurations.
For the case of the planets of TOI-238, it is not possible to measure the k2 number directly.Based on the position of the planets in the mass-radius diagram (see Fig. 15), we assume the k2 number of Earth for planet b, with an uncertainty of 10% (0.933 ± 0.093).This is roughly equivalent to assuming Earth composition.Planet c does not lie on the track of any planet with a known k2.We opt to not include this constraint in the analysis, which widens the set of potential interior configurations.To account for the uncertainties in the planetary parameters, we draw normal distributions with a σ equal to the uncertainties of each parameter.We run 1000 predictions with 1000 samples each.
For TOI-238 b we estimate a core-mass fraction of 0.43 ± 0.24, a mantle-mass fraction of 0.57 ± 0.24, and negligible mass fractions for water and the atmosphere.We estimate a coreradius fraction of 0.65 ± 0.18 and a mantle-radius fraction of 0.34 ± 0.17.We estimate radii fractions of water and volatiles smaller than 0.01.This interior structure is very similar to the interior structure of the Earth.This is not surprising, since the planet falls in the mass-radius track of the Earth and since we have used Earth's k2 number as a proxy of its composition.
For TOI-238 c we compute a core-mass fraction 0.30 ± 0.18, a mantle-mass fraction of 0.39 ± 0.25, a water-mass fraction of 0.27 ± 0.15 and a mass fraction of volatile elements < 0.018 (99.7%).We estimate a core-radius fraction 0.38 ± 0.10,   a mantle-radius fraction of 0.25 ± 0.17, a water-radius fraction of 0.26 ± 0.12, and atmosphere-radius fraction of 0.10 ± 0.06.The large uncertainty in the parameter allow for very different configurations.For a more sophisticated analysis, it would be important to obtain a more precise determination of the radius, which would require additional high-precision photometric observations.Figure 18 shows the ternary plots of the core-mantlewater mass and radii fraction configurations for both planets.
Recently, Adibekyan et al. (2021) demonstrated that the scaled density of small-sized (R p < 2 R ⊕ ) rocky planets correlates with the iron-to-silicate mass fraction ( f star iron ) of the planet building blocks, as estimated from their host star abundances.
Using the abundances of Mg, Si, and Fe as proxies for the composition of the protoplanetary disk, we calculated f star iron (Santos et al. 2015(Santos et al. , 2017) ) to be 32.7 ± 1.8%, which is similar to thevalue determined from the composition of the Sun ( f sun iron = 33.2± 1.7%).The scaled density of TOI-238 b, obtained from its radius and mass and normalized to the density of a planet with the same mass but Earth-like composition (Dorn et al. 2017), was calculated to be 0.99 ± 0.22.With these parameters, TOI-238 b aligns well with the regression line established in Adibekyan et al. (2021) (see Fig. 19).
The true structure of planets similar to TOI-238 c is very difficult to establish by their mass and radius alone.These planets lie in a region of the parameter space in which several composition tracks are close to each other.The mass-radius configuration TOI-238 c is equally consistent with water-rich world with no extended volatile atmosphere and a rocky core with a small H 2 layer.To fully disentangle these scenarios requires the characterisation of the atmosphere of the planet.
Following Kempton et al. (2018), we compute the transmission spectroscopy metric (TSM) and the emission spetroscopy metric (ESM).The TSM is a metric proportional to the expected transmission spectroscopy S/N ratio and is based on the strength of spectral features and the brightness of the host star.With our results, we estimate a TSM of 5.2 ± 1.1 for planet b, and 35 ± 11 for planet c.With these numbers, and considering the cutoffs by Kempton et al. (2018), we expect the analysis of the atmospheres of either of the planets through transmission spectroscopy to be a challenging effort.The ESM is a metric proportional to the expected S/N ratio of a JWST secondary eclipse at mid-IR wavelengths.We estimate a ESM of 3.03 ± 0.32 for planet b, and 1.82   ± 0.30 for planet c.Once again, none of them are truly favourable candidates.
We compared the results obtained for the TOI-238 system with the values obtained for all planets shown in Fig. 15, using the planetary and stellar parameters from the Planetary Systems Composite Data table from the NASA Exoplanet archive.Figure F.1 in Appendix F shows the TSM and ESM of TOI-238 compared to known low-mass transiting planets with radii and masses measured with a precision of 33%, or better.TOI-238 c ranks average in TSM compared to most planets of its class.Both planets rank average in ESM compared to other similar exoplanets.

Stellar activity of TOI-238
TOI-238 is a K2-dwarf for which we measured a log 10 (R ′ HK ) of -4.74 ± 0.27.This corresponds to an expected rotation period of 28 ± 16 days.To model the activity-induced RV variations without leaving structures that could be related to stellar activity, we needed to use a combination of two multi-dimensional GPterms: one for FWHM+RV and one for BIS+RV.We obtained different timescales in both cases, although not very different.
In our global analysis, we measured GP periods of 24.42 ± 0.86 days and 26.64 ± 0.42 days for the combinations of FWHM+RV, and BIS+RV, respectively.Both these periods are consistent with the expectations.The two measurements are within 2σ of each other, which could indicate that the difference is related to the data and/or modelling procedure, and not to a true astrophysical difference.We measured different timescales of evolution (L).For the FWHM+RV, we measure a L = 35.2+6.9 −5.8 days, while for the BIS+RV we measure L = 48.9+10.9 −8.9 days.
The signal in the bisector seems more stable than the signal in FWHM, but the differences in evolutionary timescale are within 1σ of each other.The posterior distributions are fairly similar to the imposed prior (see Fig. E.1), which suggests that the data does not carry enough information to constrain the evolutionary timescale.We tested the behaviour of modelling with wider priors.The overall results of the model were similar.However, in some cases the GP would heavily overfit the data.We also measured small differences in the length scale of the periodic component (ω).For the FWHM+RV we measure a ω of 0.71 +0.29 −0.20 , while for the BIS+RV we measure a ω of 0.12 +0.21 −0.10 .While the differences are not large, all of them point to the variations measured in the BIS being more stable than those measured in the FWHM, the same as their counterpart signals in the RV data.This picture fits with the data seen in Figs. 3 and 4, in which the FWHM shows larger changes of amplitude and shape over the course of the observing campaigns.This is particularly noticeable in the last year of observations.The origin of this difference is, however, not so clear.One possible explanation is that we are measuring the effect of different phenomena on the surface of the star, with different timescales of evolution.However, given that the data sampling is not always optimal to follow an irregularlyshaped ∼ 25 days variations, it might also be that the sampling itself is creating apparent variations by sometimes measuring the data at times of maximum variation and sometimes at times of minimum variation.The signals in FWHM and BIS having a lag with respect to one another could be the culprit behind the observed differences in shape.We estimate the stellar rotation period by averaging the two measurements we obtained in the global model.We obtain a measurement of the P rot of 25.51 ± We measure a stellar activity standard deviation (stdev) of 22.7 +6.9 −4.5 m•s −1 in FWHM and a stdev of 5.5 +1.3 −1.2 m•s −1 in the BIS time-series.In the RV data, activity-induced RV has stdev of 4.5 m•s −1 .The value is consistent with what would be expected for a star with a log 10 (R ′ HK ) of -4.74 (Suárez Mascareño et al. 2017), however, the large uncertainty in the log 10 (R ′ HK ) and the differences in methodology make the numbers difficult to compare (a standard deviation, as measured with a GP is not strictly the same as an amplitude, measured with a sinusoidal model).We find the activity-induced RV stdev related to the variations correlated to the BIS to be the dominant contribution, with 3.3 m•s −1 , followed by the variations correlated with the gradient of the BIS, with 2.9 m•s −1 .The variations linked to the BIS are anti-correlated, which is expected in the case of spot-dominated variations.The variations correlated with the FWHM account for 2.1 m•s −1 .There are no significant variations correlated with the gradient of the FWHM. Figure 20 shows the scaled variations of the FWHM, BIS, their gradients and the associated RV variations.FWHM, d/dt FWHM, BIS and d/dt BIS are scaled to their variance.All individual RV components are scaled to the variance of the global RV, to show the comparison between them.The different behaviour of the RV variations linked with FWHM and BIS, coupled with the slight difference in timescales, hint at the two indicators tracking different activity regions.The anti-correlation between the BIS and RV, combined with a correlation between its gradient and the RV, are typical signs of spotdominated variations (Queloz et al. 2001;Aigrain et al. 2012).However, the relation between the FWHM and RV variations is more reminescent of plague-dominated variations (Dumusque et al. 2014).
We do not find evidence of the presence of a magnetic cycle.However, the current dataset might not span long enough to allow such a detection.Sairam & Triaud (2022) showed that the P cyc vs P rot relation showed in Suárez Mascareño et al. ( 2016) could be used to provide a prediction on the cycle period of the star.Using the slope and zero point identified for FGK stars, we obtain the relationship log 10 (P cyc /P rot ) = 3.22 + 0.89 • log 10 (1/P rot ).We obtain a prediction for the stellar cycle period of 2340 ± 410 days (6.41 ± 1.12 years).This period is more than three times the baseline of observations, which would explain the lack of a detection.The baseline is, however, long enough to hint towards the lack of a large-amplitude cycle, given that we do not find a significant slope in any of the activity indicators.

Conclusions
We performed a high-precision radial velocity campaign on the system TOI-238, which hosts the planet candidate TOI-238.01, using the high-resolution spectrographs ESPRESSO and HARPS.We confirmed the presence of the planet, both in TESS photometry and radial velocity.   .It orbits at a distance of 0.02118 ± 0.00038 au of its host star, with an orbital period of 1.2730988 ± 0.0000029 days.It has an equilibrium temperature of 1311 K ± 28 and receives a flux 728 ± 78 times Earth's irradiation.We estimate the planet's internal structure to be consistent with Earth's internal structure.We estimate a core-mass fraction of 0.43 ± 0.24, a mantle-mass fraction of 0.57 ± 0.24, and negligible mass fractions for water and the atmosphere.
We discovered a second planet, both in TESS photometry and radial velocity.TOI-238 c is most likely a water-world or a rocky core with a small volatile envelope.It has a radius of 2.18 +0.18  −0.18 R ⊕ and a mass of 6.7 +1.1 −1.1 M ⊕ .It orbits at a distance of 0.0749 ± 0.0013 au of its host star, with an orbital period of 8.465652 ± 0.000031 days.It has an equilibrium temperature of 696 ± 15 K and receives a flux of 58.3 ± 6.2 S ⊕ .With the current data available, the interior structure of this planet is difficult to constrain.A more precise radius and, most likely, a study of its atmosphere are needed to improve the characterisation of the planet.
The RV data gives hints of a potential non-transiting planet.However, the current dataset does not allow us to properly disentangle the signal from the activity-induced RV signal.
The planets lie at both sides of the radius and density valleys.The current status of both planets seems to be compatible with both photo-evaporation and core-powered mass loss, and with migration of the outer planet after having formed beyond the ice line.
We estimate the transit spectroscopy metric (TSM) of both planets.While it should be possible to study their atmospheres via transmission spectroscopy with JWST, their comparatively low TSM indicates it would be challenging to obtain any significant detection of the planetary features.
We study the stellar activity variations of TOI-238.We measured a rotation period of 25.51 ± 0.49 days, and the activityinduced RV signal to have a standard deviation of 4.5 m•s −1 , which is consistent with the expectations of a star of its spectral type and its mean level of chromospheric activity.We find that the largest components of the activity-induced RV signal are correlated with the bisector span of the CCF and its gradient.We do not detect the presence of a magnetic cycle.However, we estimate that if there is a magnetic cycle it would have a period of ∼ 6.4 years.through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products.
This paper made use of data collected by the TESS mission and are publicly available from the Mikulski Archive for Space Telescopes (MAST) operated by the Space Telescope Science Institute (STScI).
This research has made use of the The bulk of the analysis was performed on desktop PC with an AMD Ryzen TM 9 3950X (16 cores, 2 threads per core, 3.5-4.7 GHz) and a server hosting 2x Intel (R) Xeon (R) Gold 5218 (2x16 cores, 2 threads per core, 2.3-3.9GHz).Additionally, this article used flash storage and GPU computing resources as Indefeasible Computer Rights (ICRs) being commissioned at the ASTRO POC project that Light Bridges will operate in the Island of Tenerife, Canary Islands (Spain).The ICRs used for this research were provided by Light Bridges in cooperation with Hewlett Packard Enterprise (HPE) and VAST DATA.
The ESPRESSO CCF velocities were obtained using the default K2 mask available in the ESPRESSO DRS.This mask is based on a high signal-to-noise spectrum of the low-activity K2V star HD 144628 and uses 3470 individual lines to compute the RV.For the HARPS CCF velocities we used the K5 mask available in the HARPS DRS, built using a synthetic spectrum following Queloz (1995) and Pepe et al. (2002), which uses 5129 individual lines.The CCF is computed order by order and summed to obtain the final CCF.In both cases the velocities were compute by fitting a Gaussian function to the CCF and computing its centre.We measured an RV RMS of 6.3 m•s −1 and 5.22 m•s −1 for HARPS and ESPRESSO data, respectively.We measure RV internal uncertainties of 1.6 m•s −1 and 0.73 m•s −1 , respectively.
SERVAL builds a high signal-to-noise template by co-adding all the existing observations, and performs a least-squares minimisation of each observed spectrum against the template, yielding a measure of the Doppler shift and its uncertainty.Each order is fitted separately.The final measurement is a weighted mean of the individual RVs of all the orders.SERVAL automatically masks telluric features deeper than 5% and sky emission lines using a predefined line list.We measured an RV RMS of 6.1 m•s −1 and 8.09 m•s −1 for HARPS and ESPRESSO data, respectively.We measure RV internal uncertainties of 1.2 m•s −1 and 0.66 m•s −1 , respectively.S-BART, similar to SERVAL, builds a high signal-to-noise template by co-adding all the existing observations.Then, unlike usual template-matching algorithm, S-BART uses a single RV shift to describe simultaneously the RV differences between all orders of a given spectrum and the template.The algorithm estimates the posterior distribution for the RV shifts after marginalising with respect to a linear model for the continuum levels of the spectra and template, using a Laplace approximation.From the Laplace approximation to the posterior distribution, the mean is used as the estimated RV and the standard deviation as the estimated RV uncertainty for each epoch.Using S-BART, we measured an RV RMS of 6.4 m•s −1 and 5.65 m•s −1 for HARPS and ESPRESSO data, respectively.We measure RV internal uncertainties of 1.1 m•s −1 and 0.58 m•s −1 , respectively.
Figure A.1 shows the distribution of RV measurements obtained with the three algorithms, as well as the distribution of uncertainties measured, in logarithmic scale.We use the distributions shown in A.1 to reject some poor quality measurements.First we reject the large outlier in the SERVAL data and then the few spectra where the RV uncertainty deviate from the distribution of RV uncertainties in the S-BART measurements, which shows the most compact distribution of uncertainties.The vertical dotted lines in each panel shows the threshold applied.We reject 2 HARPS spectra and 5 ESPRESSO spectra.After the data rejection we measure an RV RMS for the HARPS data of 6.3 m•s −1 , 6.1 m•s −1 and 6.2 m•s −1 , for the CCF, SERVAL and S-BART measurements, respectively.We measure an RV RMS for the ESPRESSO data of 5.30 m•s −1 , 6.56 m•s −1 and 5.55 m•s −1 , for the CCF, SERVAL and S-BART measurements, respectively.The velocities obtained with all three algorithms show similar distributions, with the SERVAL RVs having a few big outlier in the ESPRESSO RVs, which explains the significantly larger RMS.Similarly to the case of Faria et al. (2022), the S-BART measurements show the smaller uncertainties for both instru- ments, while having a very similar distribution of RVs.Based on the improved formal precision of the measurements, we opt for the S-BART RVs as our main RV measurements, although we acknowledge that any of the three methods would provide similar results.

Appendix B: Effect of tellurics in the RV data
We investigated the effect of tellurics in the radial velocity measurements by comparing the measurements we obtained with measurements corrected following Allart et al. (2022) and with uncorrected measurements.S-BART provides RV measurements corrected for the influence of telluric features.To achieve it, S-BART creates a transmittance profile using the night with the highest amount of water content in the atmosphere.From that night, it takes the water content, airmass, temperature to generate a transmittance profile with TELFIT (Gullikson et al. 2014).It then flags the regions of the transmittance profile where the flux drops under 1% of the continuum to create a telluric mask.It then shifts the mask over the complete range of barycentric earth radial velocities for the target and rejects all spectral region that fall within those limits.S-BART's masking is designed to be conservative.The regions rejected are sometimes quite large, which can impact the radial velocity precision, but should remove most of the effect of the tellurics.To calculate uncorrected velocities we switched off the telluric correction of S-BART.Allart et al. (2022) provides a procedure to correct, rather than mask, the telluric features in the spectra.It uses a lineby-line radiative transfer codes assuming a single atmospheric layer.Using the sky conditions and the physical properties of the lines it creates a high-resolution telluric spectrum, which then convolves with the instrumental resolution and samples to the instrumental wavelength grid.It finally uses a subset of selected telluric lines to fit the spectrum and subtract the telluric contribution.This allows to avoid masking regions and losing RV content in the spectrum, thus allowing for potentially better RV precision.A drawback of this algorithm is that, in its current implementation, it is only compatible with ESPRESSO and NIRPS data.We cannot use it to correct the HARPS RVs.We recompute the RVs using the corrected spectra and with S-BART's correction switched off.We obtain a slightly larger RMS of the data (6.0 m•s −1 ) with smaller error bars (0.5 m•s −1 ).
Telluric contamination in the RV measurements often appears as a coherent signal with respect to the barycentric earth radial velocity (BERV) of up to a few m/s (Allart et al. 2022).Figure B.1 shows the uncorrected ESPRESSO RVs as a function of the BERV, the corrected RVs using both methods and the difference, which corresponds to the suspected telluric effect in the RVs.While the difference between the uncorrected and corrected velocities is difficult to see, the residuals after subtracting the corrected RVs show a structure that is similar for both corrections.The biggest difference being that the correction of the algorithm by Allart et al. (2022) seems cleaner.Tellurics in our ESPRESSO data induce an RMS of ∼ 1.3 m•s −1 , large enough to take it into account.Both algorithms produce very similar corrections and provide very similar corrected RVs.The telluric RVs have an Spearman's correlation coefficient (Spearman 1904) of 0.84, while the corrected RVs of 0.99. Figure B.2 show the direct comparison between both sets of RVs.Both corrected RVs follow a 1-1 relationship almost perfectly.The telluric RVs are also well correlated, however there is a small shift between the computations of both algorithms.Nevertheless, the shift is at the level of the precision of the data.The effect of telluric contamination typically produces periodicities in the data related in some way to Earth's year.Both telluric correction algorithms tested in the ESPRESSO data provide very similar results.The algorithm by Allart et al. (2022) seems to provide a cleaner correction and allows us to avoid masking regions of the spectrum, producing higher precision RV measurements.However, the differences are small in our dataset.For consistency between the ESPRESSO and HARPS RV computations, we opted to rely on S-BART's default correction with TELFIT.

Fig. 1 .
Fig. 1.TPF files of TOI-238 for sectors 2, 29 and 69 .TOI-238 is marked with a cross-symbol, and the number 1.The rest of the numbered red symbols show the positions of the closest stars in the field.

Fig. 2 .
Fig. 2. Data used in the global analysis.Panels a, b and c show TESS data of sectors 2, 29, and 69, respectively.Panels d and e show the RV data obtained during the 2019 and 2020-2021 campaigns, respectively.Panels f and g show the FWHM time-series, and panels h and i show the bisector span time-series.
Figure D.2 shows the S-index time-series.

Fig. 3 .
Fig. 3. Analysis of the FWHM of the CCF.a, b: FWHM time-series with the best-model fit.The data is split into two panels because of a large period with no observations between the two campaigns.The shaded area shows the variance of the GP model.c: GLS periodogram of the FWHM data.The grey vertical line highlights the most significant period.d: Relationship between the RV and FWHM data.The best fit is shown when the slope is ≥2σ different from zero.e, f: Residuals of the FWHM after subtracting the best model fit.g: GLS periodogram of the residuals.h: Comparison of the RV and gradient of the FWHM model.

FollowingFig. 4 .
Fig. 4. Analysis of the bisector span of the CCF.Same as Figure 3, only with the bisector span of the CCF instead.

Fig. 5 .Fig. 6 .
Fig. 5. TESS data with the best model fit (GP-only).Panels a, b, and c, show the TESS light curve with the best GP model.Panel d shows the GLS periodogram of the TESS data.Panels e, f, and g, show the detrended TESS data.Panel h shows the TLS periodogram of the detrended light curve.
Fig. 7. TESS data with the best model fit (GP+1p).Panels a, b, and c, show the TESS light curve with the best GP+1p model.The blue lines show the transits of TOI 238.01.Panel d shows the GLS periodogram of the TESS data.Panels e, f, and g, show the detrended TESS data along with the best fit model of TOI 238.01.Panel h shows the TLS periodogram of the detrended light curve.Panels i, j, and k show the residuals after the fit of the best GP+1p model.Panel i shows the TLS periodogram of the residuals.

Figures 9
Figures 9 and 10 show the best-fit model, using the parameters obtained from the posterior distribution of the nested sampling.Table E.1 and figure E.1 show the priors and measured values for all parameters.Additionally, we repeated the modelling process in the RV and TESS data independently.We found consistent results for the planetary parameters in our combined analysis, a RV-only analysis and a TESS-only analysis.Figure11shows the comparison between the periods and times of mid transits obtained in these three models.We tested the effect of different GP Kernels (such as S+LEAF's MEP Kernel or different combinations of SHO Kernels).We found no significant differences in the planetary parameters, however the activity model would be sometimes very different.We found the ESP Kernel to be the most reliable at producing a smooth model that did not show obvious signs of overfitting.After fitting a model with a GP and two planetary signals, the periodograms of the residuals of the TESS light curve and the RV don't show any clear peak that could be attributed to a planetary signal.There is, however, a peak at 14 days in the GLS of the activity-induced component of the RV that does not fully fit as one of the harmonics of the rotation period.

Fig. 8 .
Fig. 8. RV time-series with the best model fit (GP+1p).Panels a and b show the RV time-series with the best GP+1p model.Panel c shows the GLS periodogram of the raw RV data.Panels d and e show the RV data detrended from stellar activity (i.e planetary component).Panel f shows the GLS periodogram of the detrended RV time-series.Panels g and h show the RV data detrended from the planetary component (i.e stellar activity).Panel i shows the GLS periodogram of the activity-induced RVs.Panels j and k show the residuals after the fit.Panel i shows the GLS periodogram of the residuals.

Fig. 9 .
Fig. 9. TESS data with the best model fit (GP+2p).Panels a, b, and c, show the TESS light curve with the best GP+1p model.The blue and green lines show the transits at the periods of 1.27 days, and 8.47 days, respectively.Panel d shows the GLS periodogram of the TESS data.Panels e, f, and g, show the detrended TESS data along with the best fit model of both planets.Panel h shows the TLS periodogram of the detrended light curve.Panels i, j, and k show the residuals after the fit of the best GP+1p model.Panel i shows the TLS periodogram of the residuals.

Fig. 11 .
Fig. 10.RV time-series with the best model fit (GP+2p).Panels a and b show the RV time-series with the best GP+1p model.Panel c shows the GLS periodogram of the raw RV data.Panels d and e show the RV data detrended from stellar activity (i.e planetary component).Panel f shows the GLS periodogram of the detrended RV time-series.Panels g and h show the RV data detrended from the planetary component (i.e stellar activity).Panel i shows the GLS periodogram of the activity-induced RVs.Panels j and k show the residuals after the fit.Panel l shows the GLS periodogram of the residuals.

Fig. 12 .
Fig. 12. FIP periodogram of the different 3-planet models of TOI-238.The blue solid line shows the result of the blind search for a potential third planet.The dotted orange line shows the result of the guided search.The red solid circles show the detected periodicities.The dotted green line shows the 1% FIP threshold.

Fig. 13 .Fig. 14 .
Fig. 13.Phase folded light curves.Panels a and b show the phase folded plot of the TESS light curves for planets b and c, respectively, after subtracting the best fit GP model and, in each case, the other planetary signal.Panels c and d show the residuals after the fit.

Fig. 15 .Fig. 16 .
Fig. 15.Mass-radius diagram of transiting exoplanets with well established mass and radius from the NASA Exoplanet archive along with some theoretical composition curves (Zeng et al. 2019).The purple squares highlight the transiting planets orbiting K-dwarfs.The green blobs show the mass and radius posteriors (1-and 2-σ contours) of TOI-238 b and c.The planets of the solar system are included as reference.

Fig. 17 .
Fig. 17.Mass-density diagram of transiting exoplanets with well established mass and radius from the NASA Exoplanet archive along with some theoretical composition curves.All densities are normalized to the Earth-like density track.The purple squares highlight the transiting planets orbiting K-dwarfs.The green blobs show the mass and radius posteriors (1-and 2-σ contours) of TOI-238 b and c.The shaded regions highlight the high-density and low-density groups of exoplanets.The planets of the solar system are included as reference.

Fig. 18 .
Fig. 18.Ternary plots of the interior structure of the planets of TOI-238 .The top-left panel shows the mass-fraction structure of planet b.The bottom-left panels shows its radius fraction configuration.The top-right panel shows the mass-fraction structure of planet c.The bottom-right panels shows its radius fraction configuration.
TOI-238 b is a super-Earth with f star iron (e r -M e r c u r i e s

Fig. 19 .
Fig. 19.Planet density normalized by the expected density of an Earth-like composition versus the estimated iron-to-silicate mass fraction of the protoplanetary disk for the sample studied in Adibekyan et al. (2021) and for TOI-238 b.The Earth and Mercury are indicated with their respective symbols in black.The black line indicates the original correlation for the super-Earths.
Fig. 20.Individual components of the activity model.Panel a shows the scaled variations of the FWHM model and their effect in RV.Panel b show the scaled variations of the gradient of the FWHM and their effect in RV.Panels c and d show the same for the bisector model.

Fig
Fig. A.1.Histogram of the RV and error measurements for all RV extractions tested.Top panels show the distribution of RV measurements obtained for HARPS (left) and ESPRESSO (right).Bottom panels show the distribution of internal errors for HARPS (left) and ESPRESSO (right), in logarithmic scale.The vertical lines show the thresholds used for the outlier rejection procedure.
Figures A.2 and A.3 show the visual comparison of the three RV datasets.Figure A.2 shows the RV measurements of HARPS and ESPRESSO as a function of time.The three algorithms produce virtually identical measurements except for a few points in which the SERVAL measurements deviate from the rest.This becomes more evident in figure A.3.The S-BART and CCF velocities show good correspondence, with some noise.The S-BART and SERVAL velocities show a tighter correspondence, except for a few clear outliers.The velocities obtained with all three algorithms show similar distributions, with the SERVAL RVs having a few big outlier in the ESPRESSO RVs, which explains the significantly larger RMS.Similarly to the case ofFaria et al. (2022), the S-BART measurements show the smaller uncertainties for both instru-

A
Fig. A.2.Comparison of the velocities coming from the three extraction algorithms as a function of time.Top panel shows the HARPS RV measurements from the CCF method (red), SERVAL (blue) and SBART (green).Bottom panel show the same but for the ESPRESSO measurements.
Fig. B.1.ESPRESSO Radial velocity measurements as a function of the Barycentric Earth Radial Velocity of their measurements.Top panel shows the uncorrected RV measurements.Middle panels show the telluric-corrected RV measurements using the S-BART default correction (left) and the method of Allart et al. (2022) (right).Bottom panels show the telluric-induced RV variations estimated by both methods.
Fig. B.2.Comparison of the RV derived from different telluric correction algorithms.Left panels show comparison between the corrected RVs obtained with the telluric correction of S-BART and with the correction of Allart et al. (2022) Right panels show the comparison between the telluric RVs.The grey dotted lines show the 1-1 relationship.
Figure B.3 shows both sets of telluric RVs as a function of time and their GLS periodogram.The periodicities seen in both datasets seem very similar, with the most prominent peaks in the periodogram arising at 140 and 45 days.

Fig
Fig. B.3.Telluric RVs and their periodogram.Left show the telluric RV measurements obtained with the algorithm of Allart et al. (2022) and with S-BART.Right panel show the GLS periodogram of both datasets.The grey vertical lines highlight the most significant periods, of 45 and 140 days.
Table E.1 shows the priors and measured values for all parameters.

Table 2 .
Parameters of the two planets detected orbiting TOI-238 Exoplanet Follow-up Observation Program (ExoFOP; DOI: 10.26134/ExoFOP5) website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.Funding for the TESS mission is provided by NASA's Science Mission Directorate.Part of the LCOGT telescope time was granted by NOIRLab through the Mid-Scale Innovations Program (MSIP).MSIP is funded by NSF.This research has made extensive use of the SIMBAD database, operated at CDS, Strasbourg, France, and NASA's Astrophysics Data System.This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.This work makes use of observations from the LCOGT network.The manuscript was written using Overleaf, Texmaker and VS Code.Main analysis performed in Python3 (Van Rossum & Drake 2009) running on Ubuntu (Sobell 2015) systems and MS.Windows running the Windows subsystem for Linux (WLS).Extensive use of the DACE platform 11 Extensive usage of Numpy (van der Walt et al. 2011).Extensive usage of Scipy (Virtanen et al. 2020).AstroimageJ (Collins et al. 11 https://dace.unige.ch/2017).TAPIR (Jensen 2013).All figures built with Matplotlib (Hunter 2007).
Table A.1 shows the complete set of statistics of the RV data.
Table A.1.RV dispersion and typical error of the raw RVs and the RVs after the sigma-clip.