© The Authors 2024

1 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 as to 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 (RV) surveys, the number of confirmed detections has surpassed 5500¹ (Christiansen 2022). Despite this 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 core-powered 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 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 min for selected targets and 200 s for full-frame images. TESS image data are reduced and analysed by the Science Processing Operations Center (SPOC) at NASA Ames Research Center. The SPOC pipeline performs photometry on the images to produce light curves (LCs), 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, 2020 and 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.

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.

2 Observations and data

2.1 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 2 yr 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 LC 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. 2012, 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 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 days², and TOI-238 was announced as a TESS object of interest (TOI) via the dedicated MIT TESS data alerts public website³. 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. 2010, 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 Telescopes⁴, which is a NASA funded project.

2.2 LCOGT

The TESS pixel scale is ~21″ pixel⁻¹ 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)⁵.

We observed a full transit window of TOI-238.01 continuously for 128 min 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 a 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 delta-magnitude differences between TESS-band and Pan-STARRS z-short 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 min 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⁶.

2.3 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, 2021; González Hernández et al. 2018). ESPRESSO has a resolving power of approximately R ~ 140000 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⁻¹. 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⁻¹ 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⁷.

ESPRESSO 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 website⁸.

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 min 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⁹.

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). This intervention introduced an RV offset, leading us to consider separate ESPRESSO18 and ESPRESSO19 datasets. We obtained four ESPRESSO spectra before the fibre link update of 2019 and 73 spectra after. Taking into account the expected complexity of the model (including planetary and stellar activity signals), we decided to exclude these four points as they would increase the complexity without adding meaningful new information. Previously, a potential second offset was identified due to a change in calibration lamps in late December 2020. Data reduced with DRS version 2.2.8 showed a ~2 m s⁻¹ jump in the time series of standard stars between 19 December and 26 December 2020 due to an imperfect characterisation of the new lamps. This led Faria et al. (2022) and Suárez Mascareño et al. (2023) to consider an offset before and after the lamp change. DRS version 3.0.0 removed this offset by updating its characterisation of the lamps and corrected the chromatic drift of the Fabry–Perot etalon identified by Schmidt et al. (2022). We used data reduced with DRS 3.0.0 and decided not to split the time series before and after the lamp change. Therefore, we consider a single dataset of 73 ESPRESSO spectra.

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.

2.4 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. Armstrong). 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.

2.5 Radial velocities

Using the ESPRESSO and HARPS spectra, we extracted RVs using the S–BART algorithm (Silva et al. 2022). 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⁻¹ and 5.55 m s⁻¹ for HARPS and ESPRESSO data, respectively. We measure median RV internal uncertainties of 1.1 m s⁻¹ and 0.58 m s⁻¹, 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 SERVAL (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.

2.6 CCF derived products

The presence of active regions on the stellar disk 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 υs in 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⁻¹, with a median uncertainty of 1.9 m s⁻¹. The dispersion of the bisector span time series is 6.2 m s⁻¹ with a median uncertainty of 1.9 m s⁻¹. 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.

2.7 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.

2.8 Spectral activity indicators

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

2.8.1 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: $$S=\alpha \frac{{\tilde{N}}_{\mathrm{H}}+{\tilde{N}}_{\mathrm{K}}}{{\tilde{N}}_{\mathrm{R}}+{\tilde{N}}_{\mathrm{V}}}+\beta ,$$(1)

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. Figure D.2 shows the S -index time series.

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₁₀ (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 uncertainties being a result of the large uncertainties in the photometric magnitudes shown in Table 1.

2.8.2 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.

2.8.3 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 ID 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.

3 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 $\mathcal{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 $\mathcal{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 $\left(f_{\mathrm{i}\mathrm{f}\mathrm{o}\mathrm{n}}^{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}}\right)$ 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. (2012, Adibekyan 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 $\mathcal{g}$ = 4.52 ± 0.20, [Fe/H] = −0.06 ± 0.04, Age = 4.0 ± 3.8 Gyr, M_* = 0.800 ± 0.020 M_⊙ and R_* = 0.741 ± 0.017 R_⊙. Using the SED, combined with the parallax from Gaia DR3 (Gaia Collaboration 2021), we computed a luminosity L_*/L_⊙ = 0.3297 ± 0.0058. Then, through the Stefan-Boltzmann law, we estimate a R_* = 0.747 ± 0.034 R_⊙ and M_* = 0.785 ± 0.048 M_⊙. Both sets of parameters are consistent with the stellar parameters presented in Table 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⁻¹, V = 2.83 ± 0.10 km s⁻¹, W = 19.50 ± 0.26 km s⁻¹. 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.

4 Analysis

4.1 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 Sec. 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 and Roberts 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)¹⁰, 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: $$k(\tau )=A^2\cdot \exp [-\frac{{\tau }^2}{2L^2}-\frac{{\sin }^2(\pi \tau /P_{\mathrm{r}\mathrm{o}\mathrm{t}})}{2{\omega }^2}]+({\sigma }^2(t)+{\sigma }_j^2)\cdot {\delta }_{\tau },$$(2)

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, 2006) code Dynesty (Speagle 2020; Koposov et al. 2023). All models had seven free parameters (N_free). We went for a number of live points of 10 × N_free. 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⁻¹, with the data having a median uncertainty of 2.1 m s⁻¹. For this time series, we obtain a period of ${24.3}_{-1.1}^{+1.4}$ days, with a damping timescale of ${22.7}_{-4.7}^{+5.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⁻¹, 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⁻¹, and a median uncertainty of 1.9 m s⁻¹. We obtain a period of ${25.28}_{-0.65}^{+0.82}$ days, with a damping timescale of ${52}_{-18}^{+24}$ 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 yr. 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⁻¹, 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.

Fig. 3 Analysis of the FWHM of the CCF. Panels 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. Panel c: GLS periodogram of the FWHM data. The grey vertical line highlights the most significant period. Panel d: relationship between the RV and FWHM data. The best fit is shown when the slope is ≥2σ different from zero. Panels e, f: residuals of the FWHM after subtracting the best model fit. Panel g: GLS periodogram of the residuals. Panel h: Comparison of the RV and gradient of the FWHM model.

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

4.2 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 Eq. (3). $$\begin{array}{l}\mathrm{\Delta }T{S}_1=A_1\cdot G(\tau )+B_1\cdot G^{\mathrm{\prime }}(\tau ),\\ \mathrm{\Delta }T{S}_2=A_2\cdot G(\tau )+B_2\cdot G^{\mathrm{\prime }}(\tau ),\\ \dots \end{array}$$(3)

with Δ TS_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_free.

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 min 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-min bins takes five times longer than in 10-min bins, while it does not offer any advantage given the expected characteristics of the transit.

Following Foreman-Mackey et al. (2017), the k_SHO Kernel is defined as $$k_i(\tau )=C_i^2{\mathrm{e}}^{-\tau /L}\left\{\begin{array}{l}\cosh (\eta 2\pi \tau /P_i)+\frac{P_i}{2\pi \eta \mathrm{L}}\sinh (\eta 2\pi \tau /P_i);\mathrm{i}\mathrm{f}{\mathrm{P}}_{\mathrm{i}}>2\pi \mathrm{L}\\ 2(1+\frac{2\pi \tau }{P_i});\mathrm{i}\mathrm{f}{\mathrm{P}}_{\mathrm{i}}=2\pi \mathrm{L}\\ \mathrm{c}\mathrm{o}\mathrm{s}(\eta 2\pi \tau /P_i)+\frac{P_i}{2\pi \eta \mathrm{L}}\sin (\eta 2\pi \tau /P_i);\mathrm{i}\mathrm{f}{\mathrm{P}}_{\mathrm{i}}<2\pi \mathrm{L}\end{array}\right\},$$(4)

where η = (1 – (2L/P_i)⁻²)^1/2 controls the damping of the oscillator.

This Kernel has a power spectrum density: $$S(\omega )=\sqrt{\frac{2}{\pi }}\frac{Si{\omega }_i^4}{({\omega }^2-{\omega }_i^2{)}^2+{\omega }_i^2{\omega }^2/Q^2},$$(5)

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: $$P_1=P_{\mathrm{r}\mathrm{o}\mathrm{t}},S_1=\frac{C_1}{2\cdot L}{\left(\frac{P_1}{\pi }\right)}^2,Q_1=\frac{\pi L}{P_1}.$$(6)

Our global activity model, including four different time series and three GP kernels, is then defined as follows: $$\begin{array}{rl}\mathrm{\Delta }\mathrm{F}\mathrm{l}\mathrm{u}\mathrm{x}=& G_{\mathrm{T}\mathrm{E}\mathrm{S}\mathrm{S}}(A_{\mathrm{T}\mathrm{E}\mathrm{S}\mathrm{S}},P_{\mathrm{T}\mathrm{E}\mathrm{S}\mathrm{S}},L_{\mathrm{T}\mathrm{E}\mathrm{S}\mathrm{S}},\tau ),\\ \mathrm{\Delta }FWHM=& A_{FWHM}\cdot G_{FWHM}(1,P_{FWHM},L_{FWHM},{\omega }_{FWHM},\tau ),\\ \mathrm{\Delta }\mathrm{B}\mathrm{I}\mathrm{S}=& A_{\mathrm{B}\mathrm{I}\mathrm{S}}\cdot G_{\mathrm{B}\mathrm{I}\mathrm{S}}(1,P_{\mathrm{B}\mathrm{I}\mathrm{S}},L_{\mathrm{B}\mathrm{I}\mathrm{S}},{\omega }_{\mathrm{B}\mathrm{I}\mathrm{S}},\tau ),\\ \mathrm{\Delta }\mathrm{R}\mathrm{V}=& A{11}_{\mathrm{R}\mathrm{V}}\cdot G_{}(\tau )+A{12}_{\mathrm{R}\mathrm{V}}\cdot G_{FWHM}^{\mathrm{\prime }}(\tau )\\ & +A{21}_{\mathrm{R}\mathrm{V}}\cdot G_{\mathrm{B}\mathrm{I}\mathrm{S}}(t)+A{22}_{\mathrm{R}\mathrm{V}}\cdot G_{\mathrm{B}\mathrm{I}\mathrm{S}}^{\mathrm{\prime }}(\tau ),\end{array}$$(7)

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 Fig. 2a).

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 ($\mathcal{L}$$\mathcal{U}$) prior with a range of [0, 20]. P _TESS is the period of the variability of the TESS data and uses a normal ($\mathcal{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 $\mathcal{N}$ prior 4.0 ± 0.3 (54⁺²⁰_-14 days). The zero-point of the TESS data uses a $\mathcal{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 ($\mathcal{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 $\mathcal{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 $\mathcal{N}$ priors 4.0 ± 0.3 (54⁺²⁰_-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 $\mathcal{L}$$\mathcal{U}$ 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 $\mathcal{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 $\mathcal{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 $\mathcal{N}$ priors centred at 0, with sigmas equal to the RMS of their respective dataset. The white noise component of all datasets uses $\mathcal{L}$$\mathcal{N}$ 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. Table E.1 shows the priors and measured values for all parameters.

In the case of the TESS data, we do not gain a lot of information about the stellar activity of the star. However, the GP-model 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}_{-5.8}^{+6.7}$ days, while the BIS+RV shows a periodicity of 26.53 ± 0.46 days, with a timescale of evolution of ${47.8}_{-9.6}^{+10.1}$ days. The omega parameters are ${0.76}_{-0.25}^{+0.36}$, and 0.74 (3σ), respectively. The variability shows a significant RV amplitude related to the FWHM variations (A11 _RV = ${4.7}_{-1.5}^{+2.1}$ m s⁻¹), a nonsignificant RV amplitude correlated with the gradient of the FWHM variations (A12 _RV = ${-1.6}_{-5.0}^{+4.6}$ m s⁻¹), a significant RV amplitude anti-correlated with the BIS variations (A21 _RV = ${-5.1}_{-1.4}^{+1.2}$ m s⁻¹) and an almost-significant RV amplitude correlated with the gradient of the BIS variations (A22 _RV = ${8.6}_{-4.0}^{+3.9}$ m s⁻¹). 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⁻¹ from the original 5.86 m s⁻¹ (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⁻¹ and 4.2 m s⁻¹, 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 Fig. 6f) 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 Figs. 6c,f). 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.

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. 6 RV time series with the best model fit (GP only). Panels a and b show the RV time series with the best GP model. Panel c shows the GLS periodogram of the raw RV data. Panels d and e show the detrended RV data. Panel f shows the GLS periodogram of the detrended RV time series.

4.3 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: $$T{0}_b=x[-1{]}_{\mathrm{T}\mathrm{E}\mathrm{S}\mathrm{S}}+P_b\cdot (P{h}_b-1),$$(8)

where x[-1]_TESS is the last TESS observation, P_b is the orbital period of the planet, and Ph_b is the phase, using a $\mathcal{U}$ prior [0,1]. We sampled the planet radius directly, with $\mathcal{L}$$\mathcal{U}$ prior [-5,5] in R_⊕ (corresponding to ~ 0 to 148 R_⊕). The $\mathcal{L}$$\mathcal{U}$ 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 $\mathcal{U}$ [0,1] prior. We parameterised the eccentricity as e $(\sqrt{e}\text{\hspace{0.17em}cos}(\omega ){)}^2$ + $(\sqrt{e}\text{\hspace{0.17em}sin}(\omega ){)}^2$ and ω = arctan 2$(\sqrt{e}\text{\hspace{0.17em}sin}(\omega )$, $\sqrt{e}\text{\hspace{0.17em}cos}(\omega ){)}^{}$ We then sample $\sqrt{e}\text{\hspace{0.17em}cos}(\omega )$ and $\sqrt{e}\text{\hspace{0.17em}sin}(\omega )$ 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 $\mathcal{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: $$y(t)=K(\mathrm{c}\mathrm{o}\mathrm{s}(\eta +\omega )+e\mathrm{c}\mathrm{o}\mathrm{s}(\omega ))$$(9)

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₀, the argument of periastron ω, and the eccentricity of the orbit e.

The P_orb, mid-point transit time T₀ 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 $\mathcal{L}$$\mathcal{U}$ 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 $\mathcal{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⁻¹. We measure an impact parameter <0.61 (σ). We measure a Δ lnZ of +495.7 with respect to the activity-only 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⁻¹, from the original 5.86 m s⁻¹ (74% reduction), with the ESPRESSO data showing an RMS of the residuals of 0.93 m s⁻¹ and the HARPS data showing 2.11 m s⁻¹.

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.

4.4 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 $\mathcal{N}$ prior of 8.5 ± 0.5 days. The rest of the components of the planet model are the same described in Sect. 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⁻¹. 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⁻¹, from the original 5.86 m s⁻¹ (72% reduction), with the ESPRESSO data showing an RMS of the residuals of 1.50 m s⁻¹ and the HARPS data showing 1.75 m s⁻¹. 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.

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 Fig. 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. 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 do not 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. 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.

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 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.

4.5 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⁻¹, with the ESPRESSO data showing an RMS of 1.50 m s⁻¹ and the HARPS data 1.75 m s⁻¹. These values are significantly larger than the internal errors of the RV of 0.58 m s⁻¹ and 1.10 m s⁻¹ 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 Fig. 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: $$y(t)=-K\sin (2\pi \cdot (t-T_0)/P)$$(10)

which ensures that T₀ 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 $\mathcal{L}$$\mathcal{U}$ between 10 and 100 days, and a guided fit in which we define the prior for the orbital period as $\mathcal{N}$ around 14 ± 2 days. The T₀ and mass are defined as described in Sect. 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). Although 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}_{-3.0}^{+2.2}$ 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 assessed 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.

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. 11 Time of mid-transit against the orbital period for different models. The blue ellipse shows the RV-only solution, the orange ellipse shows the TESS-only solution and the red ellipse shows the combined solution. The orange ellipse in the left panel is completely covered by the red ellipse.

Fig. 12 FIP periodogram of the different three-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.

5 Discussion

5.1 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 2 460 207.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⁻¹). 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 2 460204.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⁻¹). 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, for example 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⁻³, and 1.83 ± 0.66 g cm⁻³, 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, as there was a single identified planet candidate, but the discovery 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 TOI-238 c. Both planets fit comfortably within the usual radius vs. irradiation regions (e.g. Toledo-Padrón et al. 2020).

Figure 16 shows the histogram of the distribution of radii of well characterised 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, normalised 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 Fig. 15, normalised 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 centre, 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. 2024). 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 core-radius 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-mantle-water 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_{\mathrm{i}\mathrm{r}\mathrm{o}\mathrm{n}}^{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}}$) 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_{\mathrm{i}\mathrm{r}\mathrm{o}\mathrm{n}}^{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}}$ (Santos et al. 2015, 2017) to be 32.7 ± 1.8%, which is similar to thevalue determined from the composition of the Sun ($f_{\mathrm{i}\mathrm{r}\mathrm{o}\mathrm{n}}^{\mathrm{s}\mathrm{u}\mathrm{n}}$ = 33.2 ± 1.7%). The scaled density of TOI-238 b, obtained from its radius and mass and normalised 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 H2 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 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.

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. 14 Phase-folded RV time series. Panels a and b show the phase-folded plot of the RV time series 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 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. 16 Radii and density histograms of transiting exoplanets with well-established mass and radius from the NASA Exoplanet archive. Panel a shows the distribution of radii of all planets, with those orbiting K-dwarfs highlighted in purple. Panel b shows the distribution of densities normalised to the Earth-like density track. The teal vertical lines highlight the position of TOI-238 b and c.

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 normalised 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 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.

5.2 Stellar activity of TOI-238

TOI-238 is a K2-dwarf for which we measured a log10 (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 GP-terms: 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}_{-5.8}^{+6.9}$ days, while for the BIS+RV we measure L = ${48.9}_{-8.9}^{+10.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.20}^{+0.29}$, while for the BIS+RV we measure a ω of ${0.12}_{-0.10}^{+0.21}$ 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 irregularly-shaped ~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 ± 0.49 days. This value is fully consistent with the measured log₁₀ (R′_HK) of – 4.74 ± 0.27.

We measure a stellar activity standard deviation (stdev) of ${22.7}_{-4.5}^{+6.9}$ m s⁻¹ in FWHM and a stdev of ${5.5}_{-1.2}^{+1.3}$ m s⁻¹ in the BIS time series. In the RV data, activity-induced RV has stdev of 4.5 m s⁻¹. The value is consistent with what would be expected for a star with a log₁₀ (R′_HK) of -4.74 (Suárez Mascareño et al. 2017), however, the large uncertainty in the log₁₀ (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⁻¹, followed by the variations correlated with the gradient of the BIS, with 2.9 m s⁻¹. The variations linked to the BIS are anti-correlated, which is expected in the case of spotdominated variations. The variations correlated with the FWHM account for 2.1 m s⁻¹. 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 spot-dominated 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₁₀(P_cyc/P_rot) = 3.22 + 0.89 · log₁₀(1/P_rot). We obtain a prediction for the stellar cycle period of 2340 ± 410 days (6.41 ± 1.12 yr). 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.

Fig. 19 Planet density normalised by the expected density of an Earthlike 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.

6 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. TOI-238 b is a super-Earth with a radius of ${1.402}_{-0.086}^{+0.084}$ R_⊕ and a mass of ${3.40}_{-0.45}^{+0.46}$ M_⊕. 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 activity-induced RV signal to have a standard deviation of 4.5 m s⁻¹, 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 yr.

Acknowledgements

We thank Prof. Didier Queloz for providing us the information on the creation of the original HARPS masks. A.S.M thanks Dr. Emily L. Hunt for the advice on using VS Code as Latex editor, and the suggestions on its configuration. A.S.M., V.M.P., J.I.G.H., and R.R. acknowledge financial support from the Spanish Ministry of Science and Innovation (MICINN) project PID2020-117493GB-I00 and from the Government of the Canary Islands project ProID2020010129. The project that gave rise to these results received the support of a fellowship from the “la Caixa” Foundation (ID 100010434). The fellowship code is LCF/BQ/DI23/11990071. A.C.-G. is funded by the Spanish Ministry of Science through MCIN/AEI/10.13039/501100011033 grant PID2019-107061GB-C61. F.P.E. and C.L.O. would like to acknowledge the Swiss National Science Foundation (SNSF) for supporting research with ESPRESSO through the SNSF grants nos. 140649, 152721, 166227, 184618 and 215190. The ESPRESSO Instrument Project was partially funded through SNSF’s FLARE Programme for large infrastructures. Co-funded by the European Union (ERC, FIERCE, 101052347). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. This work was supported by FCT – Fundação para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020 – Programa Operacional Competitividade e InternacionalizaÇão by these grants: UIDB/04434/2020; UIDP/04434/2020; 2022.06962.PTDC. This work was financed by Portuguese funds through FCT – Fundação para a Ciência e a Tecnologia in the framework of the project 2022.04048.PTDC (Phi in the Sky, DOI 10.54499/2022.04048.PTDC). C.J.M. also acknowledges FCT and POCH/FSE (EC) support through Investigador FCT Contract 2021.01214.CEECIND/CP1658/CT0001. A.M.S. acknowledges support from the Fundação para a Ciência e a Tecnologia (FCT) through the Fellowship 2020.05387.BD. R.A. is a Trottier Postdoctoral Fellow and acknowledges support from the Trottier Family Foundation. This work was supported in part through a grant from the Fonds de Recherche du Québec – Nature et Technologies (FRQNT). This work was funded by the Institut Trottier de Recherche sur les Exoplanètes (iREx). This work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation. The authors acknowledge the financial support of the SNSF. D.J.A. acknowledges funding from the UKRI, (Grants ST/X001121/1, EP/X027562/1). K.A.C. and C.N.W. acknowledge support from the TESS mission via subaward s3449 from MIT. This work is based on data obtained via the HARPS public database at the European Southern Observatory (ESO). We are grateful to all the observers of the following ESO projects, whose data we are using: 072.C-0488, 085.C-0019, 183.C-0972, and 191.C-087. We are grateful to the crews at the ESO observatories of Paranal and La Silla. We acknowledge the use of public TESS data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program 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 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 (https://dace.unige.ch/) Extensive usage of Numpy (van der Walt et al. 2011). Extensive usage of Scipy (Virtanen et al. 2020). AstroimageJ (Collins et al. 2017). TAPIR (Jensen 2013). All figures built with Matplotlib (Hunter 2007). The bulk of the analysis was performed on desktop PC with an AMD Ryzen™ 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.9 GHz). 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.

Appendix A Comparison of RV extractions

We performed the RV extraction with three different algorithms. We used the CCF (Fellgett 1953; Baranne et al. 1996) built-in into the ESPRESSO and HARPS DRS (Mayor et al. 2003; Pepe et al. 2021) and two template-matching algorithms (Bouchy et al. 2001; Anglada-Escudé & Butler 2012), SERVAL (Zechmeister et al. 2018) and S–BART (Silva et al. 2022).

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⁻¹ and 5.22 m s⁻¹ for HARPS and ESPRESSO data, respectively. We measure RV internal uncertainties of 1.6 m s⁻¹ and 0.73 m s⁻¹, 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⁻¹ and 8.09 m s⁻¹ for HARPS and ESPRESSO data, respectively. We measure RV internal uncertainties of 1.2 m s⁻¹ and 0.66 m s⁻¹, 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⁻¹ and 5.65 m s⁻¹ for HARPS and ESPRESSO data, respectively. We measure RV internal uncertainties of 1.1 m s⁻¹ and 0.58 m s⁻¹, 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⁻¹, 6.1 m s⁻¹ and 6.2 m s⁻¹, for the CCF, SERVAL and S–BART measurements, respectively. We measure an RV RMS for the ESPRESSO data of 5.30 m s⁻¹, 6.56 m s⁻¹ and 5.55 m s⁻¹, for the CCF, SERVAL and S–BART measurements, respectively. Table A.1 shows the complete set of statistics of the RV data.

Fig. A.1 Histogram of the RV and error measurements for all RV extractions tested. The top panels show the distribution of RV measurements obtained for HARPS (left) and ESPRESSO (right). The 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 of Faria et al. (2022), the S–BART measurements show the smaller uncertainties for both instruments, 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.

Fig. A.2 Comparison of the velocities coming from the three extraction algorithms as a function of time. The top panel shows the HARPS RV measurements from the CCF method (red), SERVAL (blue) and SBART (green). The bottom panel show the same but for the ESPRESSO measurements.

Fig. A.3 Comparison of the velocities coming from SBART with the CCF and SERVAL. The left panel shows the HARPS (orange) and ESPRESSO (teal) RVs computed with SBART compared to the CCF RVs. The right panel shows the HARPS (orange) and ESPRESSO (teal) RVs computed with SBART compared to the SERVAL RVs.

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 line-by-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⁻¹) with smaller error bars (0.5 m s⁻¹).

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⁻¹, large enough to take it into account.

Fig. B.1 ESPRESSO radial velocity measurements as a function of the barycentric Earth radial velocity of their measurements. The top panel shows the uncorrected RV measurements. The middle panels show the telluric-corrected RV measurements using the S–BART default correction (left) and the method of Allart et al. (2022) (right). The bottom panels show the telluric-induced RV variations estimated by both methods.

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.

Fig. B.2 Comparison of the RV derived from different telluric correction algorithms. The left panels show the comparison between the corrected RVs obtained with the telluric correction of S–BART and with the correction of Allart et al. (2022) The right panels show the comparison between the telluric RVs. The dotted grey lines show the one-to-one relationship.

The effect of telluric contamination typically produces periodicities in the data related in some way to Earth’s year. 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. B.3 Telluric RVs and their periodogram. The left panel show the telluric RV measurements obtained with the algorithm of Allart et al. (2022) and with S–BART. The right panel shows the GLS periodogram of both datasets. The vertical grey lines highlight the most significant periods, of 45 and 140 days.

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.

Appendix C Telemetry data of ESPRESSO and HARPS

To prevent the possibility of false-positives or biases due to instrumental effects, we study the telemetry of HARPS and ESPRESSO to search for variations that could be correlated with scientific measurements. Suárez Mascareño et al. (2023) identified the temperature of the optical elements as a source of variations that could be mistaken for an astrophysical effect. The temperature of the Échelle grating is a typical measurement that can be found in most Échelle spectrographs, making it a good proxy to follow the thermal behaviour of the instruments. Another possible source of systematic effects is related to telluric contamination, imperfections on the correction of the Earth’s motion or imperfections of the construction of the detectors (Dumusque et al. 2015). These effects manifest as yearly signals and are usually correlated with the Earth’s barycentric radial velocity (BERV).

To determine whether there is an effect in the science measurements, we compute Spearman’s correlation coefficient (Spearman 1904) between these measurements and the temperature of the Échelle grating and the BERV, independently for each instrument. As science measurements we use the RV measurements, FWHM of the CCF, bisector span of the CCF, Ca II H&K S-index, Hα index and NaI index. Same as in Suárez Mascareño et al. (2023), we find no correlation between the temperature changes and the RV measurements, neither for HARPS nor for ESPRESSO. We do find weak but significant correlations with the FWHM of the CCF, bisector span and S-index, both for HARPS and ESPRESSO, and with Hα in the case of HARPS. We find no correlation between the BERV and the RV or FWHM measurements. We find a correlation with the bisector span of ESPRESSO, the contrast of HARPS, the S-index of ESPRESSO and the NaI index of HARPS. Table C.1 shows the measured correlation coefficients and their p-values for all measurements. We used these information in appendix D, by including a polynomial model linking the science time series and the telemetry time series in the activity model.

Appendix D Extended analysis of stellar activity

Here we include the analysis of stellar activity indicators not explicitly described in section 4.1 to search for periodicities related to stellar activity which could create false-positive detections in RV or bias the amplitude measurements in RV.

Appendix D.1 Contrast of the CCF

Figure D.1 shows the data, model, periodogram and correlation plots of the variations contrast of the CCF. The contrast is a dimensionless measure of the difference between the continuum and the minimum in the normalised CCF. To avoid numerical issues and ease visibility, we chose to work with Δ contrast × 100. We measure an RMS of the data of 24.9, and a median uncertainty of 1.1. We get a period of ${25.1}_{-2.8}^{+2.0}$ days, with an evolutionary timescale of ${29}_{-10}^{+12}$ days. We obtain a fit that is more stochastic than that for the other parameters of the CCF. After subtracting the model, we measure an RMS of the residuals of 11.3 and find no significant signals in their periodogram. We do not find a significant relationship between the contrast, or its gradient, and the RV data.

Fig. D.1 Analysis of the contrast of the CCF. Same as Figure 3, only with the contrast of the CCF instead.

Appendix D.2 S_MW index (Ca II H&K)

Figure D.2 shows the data, model, periodogram and correlation plots of the variations of the Mount Wilson S-index. To avoid numerical issues and aid the reader in visualisation, we opted to work with Δ S_MW × 100. We measure an RMS of the data of 0.064. The GLS periodogram of the data shows significant structure at periods longer than the rotation period. We attempted to include different strategies to account for this, but none of them gave good results. Among the tested models are low-order polynomials, sinusoidal models with several different priors for the periods and low-order polynomials against the measurements of the telemetry of the instruments. With the GP model, we obtain a period of 29¹⁵_-10 days, with an evolutionary timescale of 28⁺²²_-10 days. We obtain a fit that is mostly stochastic. After subtracting the model, we find no significant signals in their periodogram. We do not find a significant relationship between the S_MW, or its gradient, and the RV data.

Appendix D.3 Hα index

Figure D.3 shows the data, model, periodogram and correlation plots of the variations of the Hα index. To avoid numerical issues and ease visibility, we opted to work with Δ Hα × 100. The GLS periodogram of the data shows significant peak at the same period detected in FWHM and bisector span. With the GP model, we obtain a period of ${24.8}_{-1.4}^{+1.6}$ days, with an evolutionary timescale of ${27.6}_{-6.0}^{+7.7}$ days. We obtain a fairly stable fit. After subtracting the model. We find significant positive relationships between the variations in Hα and their gradient, and the RV data.

Appendix D.4 Na I index

Figure D.4 shows the data, model, periodogram and correlation plots of the variations of the Na I index. To avoid numerical issues and ease visibility, we opted to work with Δ Na I × 1000. The GLS periodogram of the data shows significant structure at periods longer than the rotation period. We attempted to include different strategies to account for this, but none of them gave good results. Among the tested models are low-order polynomials, sinusoidal models with several different priors for the periods and low-order polynomials against the measurements of the telemetry of the instruments. With the GP model, we obtain a period of ${34.0}_{3.8}^{+4.2}$ days, with an evolutionary timescale of ${27.3}_{-7.4}^{+10.9}$ days. After subtracting the model, we find no significant signals in the periodogram of the residuals. We find a significant relationship between the NaI data and the RVs, but not with its gradient.

Fig. D.2 Analysis of the SMW Ca II H&K index. Same as Figure 3, only with the SMW Ca II H&K index instead.

Fig. D.3 Analysis of the Hα index. Same as Figure 3, only with the SMW Hα index instead.

Appendix D.5 GP timescales

The previous models provided a variety of measurements of the timescales of stellar activity of the system. Figure D.5 shows the posterior distributions for the GP period and damping timescale. All activity indicators roughly agree with the periodicity, with the FWHM, bisector span and Hα index providing the tightest determinations of the period. In those three cases, the period determination is 1σ consistent. In these three indicators we also measured damping timescales significantly different from zero. The specific timescale varies, but the posteriors are very broad. For the case of the contrast, S-index and Na I index, we obtain a timescale that sits at the edge of the prior (2 days), which indicates a mostly stochastic behaviour of the models.

Fig. D.4 Analysis of the Na I index. Same as Figure 3, only with the Na I index instead.

Fig. D.5 GP timescales measured in the activity indicators. Top panel: Posterior distribution of the periodicity measured for each activity indicator. Bottom panel: Posterior distribution of the damping timescale measured for each activity indicator. The y-axis is cut at 0.1 to ease visibility of the distributions with a peak different from 2 (edge of the prior).

Appendix E Parameters of the models

Fig. E.1 Posterior distribution (blue distributions) of every parameter sampled in the definitive model. The red distributions show the prior used in the Nested Sampling. The yellow dashed and dotted lines show the median value and 1σ-range, respectively.

Appendix F Additional figures

Fig. F.1 TSM and ESM of TOI-238 compared with the rest of transiting exoplanets with precise measurements of mass and radius. Panels a, b and c show the TSM of all planets against their mass, density and equilibrium temperature. Panels d, e and f show the equivalent plots for the ESM.

References

All Tables

All Figures

	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.
In the text

	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.
In the text

Fig. 3 Analysis of the FWHM of the CCF. Panels 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. Panel c: GLS periodogram of the FWHM data. The grey vertical line highlights the most significant period. Panel d: relationship between the RV and FWHM data. The best fit is shown when the slope is ≥2σ different from zero. Panels e, f: residuals of the FWHM after subtracting the best model fit. Panel g: GLS periodogram of the residuals. Panel h: Comparison of the RV and gradient of the FWHM model.

In the text

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

	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.
In the text

	Fig. 6 RV time series with the best model fit (GP only). Panels a and b show the RV time series with the best GP model. Panel c shows the GLS periodogram of the raw RV data. Panels d and e show the detrended RV data. Panel f shows the GLS periodogram of the detrended RV time series.
In the text

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.

In the text

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.

In the text

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.

In the text

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.

In the text

	Fig. 11 Time of mid-transit against the orbital period for different models. The blue ellipse shows the RV-only solution, the orange ellipse shows the TESS-only solution and the red ellipse shows the combined solution. The orange ellipse in the left panel is completely covered by the red ellipse.
In the text

	Fig. 12 FIP periodogram of the different three-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.
In the text

	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.
In the text

	Fig. 14 Phase-folded RV time series. Panels a and b show the phase-folded plot of the RV time series 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.
In the text

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.

In the text

	Fig. 16 Radii and density histograms of transiting exoplanets with well-established mass and radius from the NASA Exoplanet archive. Panel a shows the distribution of radii of all planets, with those orbiting K-dwarfs highlighted in purple. Panel b shows the distribution of densities normalised to the Earth-like density track. The teal vertical lines highlight the position of TOI-238 b and c.
In the text

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 normalised 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.

In the text

	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.
In the text

	Fig. 19 Planet density normalised by the expected density of an Earthlike 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.
In the text

	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.
In the text

	Fig. A.1 Histogram of the RV and error measurements for all RV extractions tested. The top panels show the distribution of RV measurements obtained for HARPS (left) and ESPRESSO (right). The 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.
In the text

	Fig. A.2 Comparison of the velocities coming from the three extraction algorithms as a function of time. The top panel shows the HARPS RV measurements from the CCF method (red), SERVAL (blue) and SBART (green). The bottom panel show the same but for the ESPRESSO measurements.
In the text

	Fig. A.3 Comparison of the velocities coming from SBART with the CCF and SERVAL. The left panel shows the HARPS (orange) and ESPRESSO (teal) RVs computed with SBART compared to the CCF RVs. The right panel shows the HARPS (orange) and ESPRESSO (teal) RVs computed with SBART compared to the SERVAL RVs.
In the text

Fig. B.1 ESPRESSO radial velocity measurements as a function of the barycentric Earth radial velocity of their measurements. The top panel shows the uncorrected RV measurements. The middle panels show the telluric-corrected RV measurements using the S–BART default correction (left) and the method of Allart et al. (2022) (right). The bottom panels show the telluric-induced RV variations estimated by both methods.

In the text

	Fig. B.2 Comparison of the RV derived from different telluric correction algorithms. The left panels show the comparison between the corrected RVs obtained with the telluric correction of S–BART and with the correction of Allart et al. (2022) The right panels show the comparison between the telluric RVs. The dotted grey lines show the one-to-one relationship.
In the text

	Fig. B.3 Telluric RVs and their periodogram. The left panel show the telluric RV measurements obtained with the algorithm of Allart et al. (2022) and with S–BART. The right panel shows the GLS periodogram of both datasets. The vertical grey lines highlight the most significant periods, of 45 and 140 days.
In the text

	Fig. D.1 Analysis of the contrast of the CCF. Same as Figure 3, only with the contrast of the CCF instead.
In the text

	Fig. D.2 Analysis of the SMW Ca II H&K index. Same as Figure 3, only with the SMW Ca II H&K index instead.
In the text

	Fig. D.3 Analysis of the Hα index. Same as Figure 3, only with the SMW Hα index instead.
In the text

	Fig. D.4 Analysis of the Na I index. Same as Figure 3, only with the Na I index instead.
In the text

	Fig. D.5 GP timescales measured in the activity indicators. Top panel: Posterior distribution of the periodicity measured for each activity indicator. Bottom panel: Posterior distribution of the damping timescale measured for each activity indicator. The y-axis is cut at 0.1 to ease visibility of the distributions with a peak different from 2 (edge of the prior).
In the text

	Fig. E.1 Posterior distribution (blue distributions) of every parameter sampled in the definitive model. The red distributions show the prior used in the Nested Sampling. The yellow dashed and dotted lines show the median value and 1σ-range, respectively.
In the text

	Fig. F.1 TSM and ESM of TOI-238 compared with the rest of transiting exoplanets with precise measurements of mass and radius. Panels a, b and c show the TSM of all planets against their mass, density and equilibrium temperature. Panels d, e and f show the equivalent plots for the ESM.
In the text