© ESO 2022

1 Introduction

The number of discovered exoplanets has grown tremendously over the last two decades thanks to forefront space- and ground-based instruments and the combination of successful techniques such as the transit photometry and the spectroscopic radial velocity (RV) methods. Transiting planetary systems with bright parent stars offer a unique opportunity to determine the bulk density of the planets and to explore the atmospheres of these planets (Madhusudhan 2019). Spectrographs delivering high-precision RV measurements combined with the superb quality of photometric light curves (LCs) obtained from space-based missions are now approaching the golden objective of finding and characterizing planets the size of the Earth (e.g., Luque et al. 2019; Lillo-Box et al. 2020a).

In recent years, the search for new exoplanets has been focused on low-mass stars because they are the most numerousstars in the Galaxy and allow lower-mass planets to be identified for a given level of observational precision (Dressing & Charbonneau 2013, 2015; Reylé et al. 2021). At present, about 800 M-type stars are known to harbor planetary systems, accounting for more than 1000 exoplanets. Considering planetary systems with transiting planets that have a robust mass and radius determination, this number is reduced to about 40–50 exoplanets orbiting M-type dwarfs (Martínez-Rodríguez et al. 2019; Trifonov et al. 2021), 26 of which are in nine different multi-planetary systems. These planets typically have masses in the Earth to super-Earth mass regime in rather circular orbits around their low-mass parent stars with short orbital periods. Very interestingly, a few of these systems (~1%) also harbor additional, more massive planets in the Neptune- to Jupiter-mass regime that pose a challenge to the core accretion theory of planet formation (Laughlin et al. 2004; Howard et al. 2010; Ida et al. 2013; Morales et al. 2019).

Here, we present the confirmation of two transiting planets orbiting the early-M dwarf TOI-1238, whose transit signals were first detected by NASA’s Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015). Their masses and radii are derived through a joint modeling of the TESS photometry and follow-up RV observations with the Calar Alto high-Resolution search for M dwarfs with Exoearths with Near-infrared and optical Échelle Spectrographs (CARMENES; Quirrenbach et al. 2016, 2018; Reiners et al. 2018a), together with extensive ground-based photometric follow-up observations with telescopes at the Joan Oró and Sierra Nevada observatories. TOI-1238 thus becomes the tenth known transiting multi-planetary system around an M dwarf. The signal of a third companion, with no transit detected in the TESS LC, is evident from the CARMENES RV time series. This candidate is located >30 times farther away from its parent star than the outer of the two transiting planets. This work is part of the TESS follow-up program within the CARMENES guaranteed time observations (GTO) survey, which aims to validate the planetary nature of transit events detected around M dwarfs (Luque et al. 2019; Kemmer et al. 2020; Nowak et al. 2020; Dreizler et al. 2020; Bluhm et al. 2020, 2021; Soto et al. 2021).

In Sect. 2 we present all space- and ground-based observations. The properties of the host star, TOI-1238, are introduced in Sect. 3. With the rotation period of the parent star being constrained by our data, we proceeded to analyze the photometric and spectroscopic observations with the aim of characterizing the transiting planetary candidates and searching for new planet candidates. The properties of the newly confirmed and discovered planets orbiting TOI-1238 are given in Sect. 4. A brief discussion on the implications of these findings is presented in Sect. 5, and the conclusions appear in Sect. 6.

2 Observations

2.1 TESS photometry

TESS is an all-sky photometric transit survey with the main objective of finding planets smaller than Neptune orbiting bright stars that are amenable to follow-up observations leading to the determination of planetary masses and atmospheric compositions. In its primary mission, TESS has conducted high-precision photometry of more than 200 000 stars over two years of observations until 4 July 2020. All targets were made available to the community as target pixel files (TPFs) and calibrated LCs. TESS LC files include the timestamps, simple aperture photometry (SAP) fluxes, and pre-search data conditioned simple aperture photometry (PDCSAP) fluxes (Smith et al. 2012; Stumpe et al. 2012, 2014). The SAP flux is the flux after summing 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. 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 second extension header data unit (HDU) consists of an image that contains the registered pixels for the target star and marks the pixels used for the optimal aperture photometry.The TESS detector bandpass spans from 600 to 1000 nm and is centered on the traditional Cousins I band (786.5 nm). This wide red optical bandpass was chosen to reduce photon-counting noise and to increase the sensitivity to small planets transiting cool, red stars.

TOI-1238 (TIC 153951307) was observed by TESS at the 2 min short-cadence integrations in sectors 14, 15, 21, and 22 during the TESS primary mission (see Table 1). All sectors were processed by the Science Processing Operations Center (SPOC) pipeline (Jenkins et al. 2016) and searched for transiting planet signatures with an adaptive, wavelet-based transit detection algorithm (Jenkins 2002; ?). In a first analysis of sectors 14 and 15, conducted separately by the SPOC pipeline, a planet candidate was identified at an orbital period of 3.294 d, and TOI-1238 was announced as a TESS object of interest (TOI) via the dedicated MIT TESS data alerts public website¹ on 17 October 2019. In a subsequent joint analysis of sectors 14 and 15, a second transiting planet candidate at a shorter separation from its parent star with orbital period of 0.764 d was announced on 15 November 2019. Both planet candidates were fitted with limb-darkening transit models (Li et al. 2019) by SPOC and successfully passed all the diagnostic tests (Twicken et al. 2018)performed on them: comparison of even and odd transits to screen against eclipsing binaries and ghost diagnostic tests to helprule out scattered light or background eclipsing binaries. From the difference image centroiding test, the transit source is coincident with the core of the stellar point spread function, so it is clear that the transit events are associated with the target and not with nearby bright stars. The SPOC analysis located the source of the transit signatures to within 1.5 ± 2.3 arcsec, and 4.0 ± 4.0 arcsec, respectively.No further transiting planet signatures were detected in the residual LC by the SPOC. The candidates have estimated radii of 2.29 ± 0.29 R_⊕ and 1.3 ± 3.5 R_⊕ for the outer and inner objects, respectively. The later TESS observations in sectors 21 and 22 did confirm the presence of shallow transit-like signatures at these two different periodicities.

All LCs and TPF files were downloaded from the Mikulski Archive for Space Telescopes, which is a NASA funded project. Our first step was to verify that the SAP and PDCSAP fluxes automatically computed by the pipeline are useful for scientific studies by confirming that no additional bright source is contaminating the aperture photometry.Figure 1 displays the TPFs of TOI-1238 for all four sectors using the publicly available tpfplotter² code (Aller et al. 2020). This code overplots all sources from the Gaia Data Release 2 (DR2) catalog (Gaia Collaboration 2018) with magnitude contrast up to Δm = 6 mag on top of the TESS TPFs. There are no additional Gaia sources within the photometric aperture around TOI-1238 automatically selected by the pipeline. Therefore, we considered the extracted TESS LC to be free of contamination from nearby stars. The SAP and PDCSAP LCs of TOI-1238 are illustrated in the top and middle panels of Fig. 2.

2.2 ASAS-SN photometry

TOI-1238 was photometrically observed by the All-Sky Automated Survey for Supernovae (ASAS-SN) project (Shappee et al. 2014; Kochanek et al. 2017). ASAS-SN currently consists of 24 telescopes distributed around the globe. These telescopes are used to survey the entire visible sky every night down to 18 mag using robotic telescopes with a diameter of 14 cm; the large number of telescopes minimizes the impact of poor weather. The ASAS-SN data of TOI-1238 were taken with two different passbands: V (cameras “bc” and “bd”) and g′ (cameras “bt” and “bs”). A total of 1065 and 1018 measurements in the V and g′ filters were obtained between 6 May 2015 and 28 November 2018, and between 29 November 2018 and 5 September 2020, respectively. TOI-1238 was typically observed on two to four occasions per filter and per day. The ASAS-SN aperture photometry is calibrated using the AAVSO Photometric All-Sky Survey (APASS) catalog (Henden et al. 2015). The ASAS-SN detectors have a pixel scale of 8′′ projected onto the sky; the images of the stars have a typical full-width-half-maximum of about 15′′, so there could be some stellar blending particularly for crowded fields. We downloaded the ASAS-SN V - and g′ -band photometry of TOI-1238. The different cameras were treated separately in order to minimize any possible systematics. The original photometry contained several outliers and we applied a 2.5σ-clipping algorithm to clean the various ASAS-SN LCs. The root mean square (rms) of the data is ~0.4 mag. Both V - and g′ -band photometrictime-series are shown in the four top panels of Fig. 3.

Fig. 1 TPFs of TOI-1238 (cross symbol) for TESS sectors 14, 15, 21, and 22 (left to right). The electron counts are color-coded. The TESS optimal photometric aperture per sector used to obtain the SAP fluxes is marked with red squares. The Gaia DR2 objects with G-band magnitudes down to 6 mag fainter than TOI-1238 are labeled with numbers (TOI-1238 corresponds to number 1), and their scaled brightness, based on Gaia magnitudes, is shown by red circles of different sizes (see figure inset). The pixel scale is 21 arcsec pixel−1.

Fig. 2 TESS LCs. Top and middle panels: SAP and PDCSAP LCs of TOI-1238 as obtained by TESS in sectors 14, 15, 21, and 22 and processed by SPOC (left to right). Bottom panel: PDCSAP LCs with the detrending GP models overplotted (black line). The GP models were used to remove any residual instrumental effect (after masking the position of all transits). As illustrated, only the data of sector 22 required this additional correction to produce a flat LC.

2.3 TJO photometry

We observed TOI-1238 with the 0.8 m Joan Oró Telescope (TJO; Colomé et al. 2010) at the Montsec Astronomical Observatory in Lleida, Spain, from 19 June 2020 through 19 February 2021. A total of 507 images (5–10 images per night) were obtained using the Johnson R filter, an exposure time of 60 s, and theLAIA³ imager, a 4k × 4k charge-coupled device (CCD) with a field of view of 30′ and a plate scale of 0.4′′ pixel⁻¹. Raw frameswere corrected for dark current and bias and were flat-fielded using the ICAT⁴ pipeline (Colome & Ribas 2006) of the TJO. The aperture photometry of TOI-1238 was extracted with the AstroImageJ software (Collins et al. 2017) by using an optimal aperture size that minimized the rms of the resulting relative fluxes. To derive the differential photometry of TOI-1238, we selected the 30 brightest comparison stars in the field that did not show any variability. Then, we employed our own pipelines to remove outliers and measurements affected by poor observing conditions or low signal-to-noise ratio. The rms of the TJO differential photometry after the removal of outliers is 8 mmag. The TJO R-band LC is presented in the fifth panel of Fig. 3.

Fig. 3 Photometric time series of TOI-1238 using ASAS-SN (aperture photometry), TJO (differential photometry), and SNO (differential photometry) data.

2.4 SNO photometry

We also monitored TOI-1238 using the 90-cm T90 telescope at the Observatorio de Sierra Nevada (SNO), Spain. The T90 telescope is equipped with a VersArray 2k × 2k CCD camera, which delivers images with a field of view of 13.2′ ×13.2′ (Amado et al. 2021). The observations were collected with both Johnson V and R filters at a total of 31 observing epochs over 88 days in the period January–April 2021. Per observing epoch, we typically obtained 20 individual measurements per filter, each with an integration time of 20 s. We obtained aperture photometry from the unbinned frames, which were bias subtracted and properly flat-fielded beforehand. We explored different aperture sizes and the circular aperture with radius of 12 pixels produced the highest-quality photometric LC for both filters. We identified six stars in the field that were non-variable in the time interval of the SNO observations and could act as reference stars to produce the differential photometry of TOI-1238. In Fig. A.1, the target and the reference stars are marked. These reference stars are 0.2–2.5 mag fainter than TOI-1238. From the final SNO photometric time series, we removed low-quality data obtained under poor observing conditions or at very high air masses. The final SNO LCs have an rms of 4 mmag for the two filters. The V - and R-band SNO photometric time series are shown in the bottom panels of Fig. 3.

2.5 Gemini speckle imaging

TOI-1238 was observed on 17 February 2020 using the ‘Alopeke speckle instrument on the Gemini North 8.1 m telescope⁵. ‘Alopeke provides simultaneous speckle imaging in two bands (562 and 832 nm) with output data products including a reconstructed image with robust contrast limits on companion detections (e.g., Howell et al. 2016). Five sets of 1000 × 0.06 s exposures were collected and subjected to Fourier analysis in our standard reduction pipeline (see Howell et al. 2011). Figure 4 shows our final 5σ contrast curves and the 832 nm reconstructed speckle image.

2.6 CARMENES spectroscopy

After the announcement of TOI-1238 as a TESS object of interest, we obtained a total of 55 high-resolution spectra using the fiber-fed-échelle spectrograph CARMENES from 9 May 2020 through 5 April 2021. The first observing season (9 May 2020 through 1 August 2020) has the largest amount of CARMENES data (48 spectra) and was mainly used to characterize the transiting planets. On three different nights, we acquired two CARMENES spectra per night to address the aliasing issues described in Sects. 3.2.3 and 4.4. The second CARMENES season (23 January 2021 through 5 April 2021) has seven spectra and was intended to monitor the RV long-term trend of TOI-1238. The TJO and SNO photometry and the CARMENES data were thus contemporaneous. CARMENES is installed at the 3.5 m telescope at the Calar Alto Observatoryin Almería (Spain). It was specifically designed to deliver high-resolution spectra at optical (resolving power R = 94 600) and near-infrared (R = 80 500) wavelengths covering the wavelength range from 520 to 1710 nm without any significant gap. CARMENES has two different channels, one for the optical (the VIS channel) and one for the near-infrared (the NIR channel) with the break at 960 nm (Quirrenbach et al. 2016). All data were acquired with integration times of 1800 s (which is the maximum exposure employed for precise RV measurements) and followed the data flow of the CARMENES GTO program. Only the spectrum taken on BJD = 2 459 278.65 was acquired with a lower exposure time, 1200 s. In this work, we used only the CARMENES VIS data because the NIR channel delivers relative RVs with less precision than the VIS channel (Reiners et al. 2018b; Bauer et al. 2020), and the expected RV amplitudes of the transiting planets (≤ 5 m s⁻¹) are notably lower than the median RV precision obtained in the NIR (mean error bar of 24 m s⁻¹). Nevertheless, the CARMENES NIR RVs of TOI-1238 are analyzed briefly in Sect. 4.2.2 in relation to an outer companion candidate with large RV amplitude. From now on, when we mention CARMENES data (without the “VIS”) we refer to the VIS channel, unless specified otherwise.

CARMENES raw data are automatically reduced with the caracal pipeline (Caballero et al. 2016b). Relative RVs are extracted separately for the VIS and NIR channels using the serval software (Zechmeister et al. 2018), where each échelle order is fitted separately. The final RV per epoch is computed as the weighted RV mean over all échelle orders. Our first step was to remove RV outliers from the CARMENES time series by applying a 2.5σ-clipping algorithm and to remove data points with very large error bars (more than three times the mean error bar size, 3.93 m s⁻¹). The final CARMENES data set has 53 RV measurements out of a total of 55 RVs, and it is displayed in the top panel of Fig. 5. The 53 individual CARMENES relative RVs used in our analysis and their associated uncertainties are listed in Table B.2. For completeness, we also give in Table B.3 the deviating velocity measurements not used here. The CARMENES RVs show a marked downward trend, which will be discussed in the next sections.

The CARMENES serval pipeline also provides measurements for a number of spectral features that are considered indicators of stellar activity: pseudo-equivalent widths of He I D₃, He I λ10 833 Å, Hα, Pa_β, Na I D doublet, and the Ca II infrared triplet (IRT), all of which may have a chromospheric component in active M dwarfs. The pipeline also provides measurements of the strength of two photospheric molecular bands: TiO and VO (see also Zechmeister et al. 2018; Schöfer et al. 2019). At high spectral resolution, the profile of the stellar lines may change due to photospheric and chromospheric activity, which has an impact on accurate RV measurements. It is important to disentangle the effects of stellar activity from the Keplerian signals. We used the following CARMENES spectral indices: Hα, Ca II IRT λλλ8498, 8542, 8662 Å, the differential line width (dLW), and the chromatic index (CRX), all of which were defined by Zechmeister et al. (2018). The CRX determines the RV– logλ correlation, and it is used as an indicator of the presence of stellar active regions. The measured CARMENES activity indicators of TOI-1238 are given in Table B.2 and are shown in Fig. 5, where data points deviating significantly (more than 2.5 σ) from the sequence or data with very large error bars have been removed from the time series. We used these σ-clipped activity indices to analyze the properties of the parent star. Typically, one to three data points were rejected per CARMENES activity indicator. These rejected data correspond to the short-integration-time spectrum and the spectra acquired at high air masses. For completeness, they are provided in Table B.3.

Fig. 4 Relative 5σ detection limit magnitude as a function of the angular separation for Gemini (‘Alopeke) at 562 nm (blue line) and 832 nm (red line). The inset figure shows the reconstructed ‘Alopeke high-contrast, high-resolution speckle image at 832 nm.

Fig. 5 CARMENES VIS RV time series (blue dots) and various spectroscopic activity indicators (green dots) of TOI-1238.

3 TOI-1238

3.1 Stellar properties

TOI-1238 (2MASS J13253177+6850106) is not an especially well-studied star, and does not appear very often in the literature except in large catalogs. In Table 2, we provide all optical, near- and mid-infrared photometry extracted from the Gaia, 2MASS⁶, and AllWISE⁷ catalogs (Gaia Collaboration 2016; Skrutskie et al. 2006; Wright et al. 2010), together with the Gaia Early Data Release 3 (EDR3) trigonometric parallax, proper motion, and equatorial coordinates (Gaia Collaboration 2021). TOI-1238 is located at a distance of 70.642 ± 0.061 pc from the Sun and has optical and infrared colors typical of K7–M0 type stars with no obvious evidence of extinction at short wavelengths. The location of TOI-1238 in the Gaia color-magnitude diagram is shown in Fig. 6; it is obvious that this star is not overluminous (i.e., it is not young) and has absolute magnitudes compatible with main sequence K7–M0 star. All stellar sequences shown in the figure were built using Gaia photometry and parallaxes (see Luhman 2018 for the young isochrones and Cifuentes et al. 2020 for the main sequence).

We first estimated the physical parameters of TOI-1238 empirically by fitting the observational spectral energy distribution constructed from cataloged broadband photometric magnitudes, following the prescription of the VOSA⁸ tool (Bayo et al. 2008). The available photometry includes the near-ultraviolet (NUV) magnitude from GALEX⁹ (Bianchi et al. 2017), the Gaia BP and RP magnitudes, the gri magnitudes from the Sloan Digital Sky Survey (Blanton et al. 2017), the BV magnitudes of the APASS catalog (Henden et al. 2015), the JHK_s magnitudes from 2MASS, and the W1, W2, and W3 magnitudes from AllWISE (Wright et al. 2010; Mainzer et al. 2011). We employed the BT-Settl models (Allard et al. 2012), which are computed for a wide range of effective temperatures (T_eff), surface gravities and atmospheric metallicities with steps of 100 K, 0.05 dex and 0.05 dex, respectively. The photometric spectral energy distribution of TOI-1238 covers the wavelength interval 0.23–22.09 μm. We assumed no interstellar extinction due to the proximity of the star and the presence of no obvious reddening in the stellar colors at short wavelengths. From this analysis, we obtained T_eff = 3900 ± 100 K, log g = 5.5 ± 0.5 dex, and [Fe/H] = 0.0–0.5 dex. The high surface gravity is consistent with a main sequence star. The observed spectral energy distribution of TOI-1238 and its best fit are shown in Fig. 7. The GALEX NUV emission of the star at 230 nm is higher than the expected photospheric emission, which is a signpost of the presence of a stellar chromosphere/corona (Cifuentes et al. 2020). From our fit we derived a stellar luminosity of log L/L_⊙ = − 1.0824 ± 0.0024 dex (Table 2).

We also calculated the stellar parameters of TOI-1238 using the CARMENES VIS and NIR spectra and the grid of PHOENIX-ACES synthetic models, following the method of Passegger et al. (2019). The resulting values (listed in Table 2) are T_eff = 4089 ± 54 K, log g = 4.63 ± 0.06 dex, and [Fe/H] = +0.31 ± 0.19 dex; at the 2σ level, they are consistent with those derived from the photometric spectral energy distribution and indicate that the star is likely more metal-rich than the Sun. We used these values (and an age >0.8 Gyr; see below) to determine the radius and mass of TOI-1238 at R = 0.58 ± 0.02 R_⊙ and M = 0.59 ± 0.02 M_⊙, following the mass–radius relation of Schweitzer et al. (2019), which was based on eclipsing binary stars. From the spectral fitting, we also derived an upper limit on the projected rotational velocity of the star at v sin i ≤ 2 km s⁻¹.

The Galactic space velocities UV W of TOI-1238 were derived using the Gaia coordinates, proper motion, and RV with the formulation developed by Johnson & Soderblom (1987). The UV W components in the directions of the Galactic center, Galactic rotation, and north Galactic pole, respectively, are given in Table 2. We note that the right-handed system is used and that we did not subtract the solar motion from our calculations. The uncertainties associated with each space velocity component were obtained from the observational quantities and their error bars. TOI-1238 has kinematics typical of “the field” (it does not appear to belong to any known young stellar moving group) indicating a likely age of >0.8 Gyr. This result was also confirmed by using the BANYAN Σ code (Gagné et al. 2018).

Fig. 6 Location of TOI-1238 (red circle) in the Gaia color-magnitude diagram. The young stellar sequences of the β Pic (blue) andTucana-Horologium (green) moving groups and the Pleiades star cluster (purple) are taken from Luhman (2018) and the main sequence of field stars (black) from Cifuentes et al. (2020). The gray area represents the dispersion observed among stars on the main sequence.

Fig. 7 Photometric spectral energy distribution of TOI-1238 (red dots) and photometry of the best BT-Settl model (Teff = 3900 ± 100 K, log g = 5.5 ± 0.5 dex, and [Fe/H] = 0.0–0.5 dex, Allard et al. 2012). The reddest data point (W4) is an upper limit.

3.2 Stellar variability and rotation period

M dwarfs can exhibit a large range of stellar activity levels from inactivity (Jeffers et al. 2020) to active states with intensities orders of magnitude greater than what is commonly observed in the Sun. The most active stars show inhomogeneities on their surface, such as dark starspots corotating with the star (Leto et al. 1997; Osten et al. 2005). It is well known that stellar activity can cause line asymmetries that hinder the very precise measurement of the line center, and consequently induce an apparent RV shift, which may mimic a Keplerian signal or hide the presence of a real planet orbiting the star (Barnes et al. 2015). Therefore, before exploring the CARMENES RVs in search for or confirming planets, we analyzed all photometric and spectroscopic index time series of established activity indices available to us in order to identify the characteristic frequencies of the stellar variations and, if possible, the rotation period of TOI-1238. The spectrophotometric spectral type of TOI-1238 is K7–M0. Following previous works (e.g., McQuillan et al. 2014; Suárez Mascareño et al. 2016, 2018; González-Álvarez et al. 2019), early-M dwarfs typically exhibit measurable rotation periods ranging from 20 d through ~100 d, and magnetic activity cycles of several hundred to thousand days.

3.2.1 TESS light curves

We analyzed the SAP fluxes (uncorrected for instrumental features; Twicken et al. 2010; Morris et al. 2020) of all four TESS sectors. We assumed that each sector has a different flux offset in order to bring all sectors to a common reference. In addition, we also assumed intra-sector flux offsets to account for possible drifts in the photometry of the intervals before and after the TESS data downlink time. Therefore, in practice we studied the TESS LCs divided into eight chunks. The top panel of Fig. 2 shows that the SAP fluxes do not suffer from strong instrumental effects and some low-amplitude variability becomes evident from the data. We modeled the TESS SAP fluxes using two sinusoidal functions with four free parameters: two different amplitudes, one period at the rotation period of the star (P_rot) and other period at approximately half of the stellar rotation period (P_{1∕2 rot}). We consider that this approach can take into account the stellar variability produced by spots located at different latitudes of the stellar surface and also accounts for any possible differential rotation from cycle to cycle. As for the main parameter, we allowed P_rot to vary between 20 and 60 d with a uniform distribution; all other parameters (LC amplitudes and offsets) were explored from initial uniform distributions with a wide range of possible values. The fit was performed using the juliet python package (Espinoza et al. 2019), which uses a nested sampling algorithm in the framework of Bayesian analysis. The MultiNest library (Feroz et al. 2009) was employed via the PyMultiNest wrapper (Buchner et al. 2014) to explore the parameter space and to efficiently compute the Bayesian model log-evidence, $\ln \mathcal{Z}$.

The simulations converge on rotation periods of P_rot = 40 ± 8 d and P_{1∕2 rot} = 18 ± 2 d. The model with the rotation period and its half value has a higher statistical significance than the one-sinusoid function model (Δ log-evidence ≈800). The flux amplitude of the variability associated with P_{1∕2 rot} is about the half of that of the main periodicity. Figure 8 illustrates the TESS SAP fluxes together with the best model, and Fig. 9 shows the resulting corner plot of the fit (the flux offsets were excluded for the clarity of the figure). We remark that all priors were uniform distributions but the distributions of the posteriors are quite Gaussian-like. From Fig. 8, the best model nicely reproduces the flux variability from sector to sector, with a few exceptions: the beginning of the second parts of sectors 21 and 22 (BJD ~2 458 888 and ~2 458 916). This is likely due to instrumental effects that remain uncorrected in the SAP fluxes. We can make an association ascribing the 40 d periodicity to the rotation period of TOI-1238, although as discussed in the next sections, this value is likely affected by differential rotation.

Fig. 8 TESS photometry. Top panels: TOI-1238 TESS SAP photometric fluxes for the four different sectors (each half sector is shown with a different color) normalized to a common reference by fitting sinusoidal functions. The best model fit is plotted asthe solid black line, and the binned photometric data are indicated with gray dots. Bottom panels: residuals shown as a function of time.

Fig. 9 Posterior distributions of the sinusoidal fit to the TESS SAP fluxes. Vertical dashed lines correspond to the 16 and 84% quantiles (1σ uncertainty). The red line stands for the mean values of the distributions.

Fig. 10 GLS periodograms of the ASAS-SN V and g, TJO R, and SNO V and R LCs shown in the time space between 0.5 and 1000 d. In all panels, the horizontal dashed lines indicate FAP levels of 10 (blue), 1 (orange), and 0.1% (green). The orbital periods of the two transiting planets are marked with vertical blue lines. The highest peak of the CARMENES RV GLS periodogram is shown with a violet line. The greenish band indicates the region where most of the spectroscopic activity indicators have their highest GLS peaks. In the fourth panel, the vertical dashed line indicates the periodicity of the signal removed from the data (see text).

3.2.2 ASAS-SN, TJO, and SNO light curves

We computed the generalized Lomb-Scargle (GLS) periodograms (Zechmeister & Kürster 2009) for all ASAS-SN, TJO, and SNO LCs. The periodograms as a function of time and frequency are shown in Figs. 10 and A.3, respectively. We included three different false alarm probability (FAP) levels and the location of the two transiting planets at 0.764 and 3.294 d as a reference. We also indicate the location of the highest peak of the CARMENES RV GLS periodogram for completeness (see Sect. 4.2). In the time line, the ASAS-SN V -band photometry is the oldest data set, then ASAS-SN g-band photometry was taken with about 1year overlap with the V -band photometry. The TJO R-band photometrictime series followed with a small overlap of 1 month with the ASAS-SN g-band data. Finally, the SNO photometry consists of the most recent data and they overlap for a period of ~1 month with the TJO observations. For all LCs we evaluated the GLS periodograms from the Nyquist frequency down to 0.001 d⁻¹. For plottingpurposes, we show the GLS periodograms in the time interval between 0.5 d (2 d⁻¹) and 1000 d (0.001 d⁻¹). The peak at1 d, although visible in all periodograms of Fig. 10, is ascribed to the effects of the observing window (typical of ground-based observations) and will not be considered further.

The ASAS-SN data of the same filter (V and g bands) but taken with different cameras were combined by applying a relative offset so that any possible systematic effect is removed. The highest peak of the ASAS-SN V -band GLS periodogram (top panel of Fig. 10) occurs at ~70 d with a second significant peak at ~40 d. The ASAS-SN g-band data (second panel of Fig. 10) has the strongest peak at around 90 d and a second strong peak at around 45 d. Given the ratio 2:1 between the two peaks for each filter, it is likely that one of the peaks is the first harmonic of the other. The GLS periodogram of the TJO R-band photometry is rather puzzling (third panel of Fig. 10): the highest, significant peak is located at ~133 d, which we ascribed to a sampling effect. In order to see whether this strong peak is veiling other peaks at shorter periodicities (thus in better agreement with the TESS and ASAS-SN data), we subtracted the 133 d signal from the original data and performed a new GLS analysis on the residuals. The resulting GLS periodogram is shown in the fourth panel of Fig. 10: there is one significant peak at around 30 d and another peak at around 19 d, which are near the periods found in our previous analysis of the TESS LC (the TJO data were acquired immediately after TESS sector 22). The GLS periodograms of the SNO V - and R-band LCs are very much alike (two bottom panels of Fig. 10), indicating that a period of ~43 d and a forest of smaller peaks between 10 and 30 d are present in both data sets. The rms values of all ground-based photometric time series are summarized in Table 3.

For completeness, we combined the TJO and SNO R-band LCs and searched for the characteristic period of TOI-1238 by fitting a sinusoidal function, following the same procedure as for the TESS data. The priors on the stellar rotation period are uniform in the interval 30–150 d. The best model has a period of 40.3 ± 0.4 d, which agrees with all values previously derived from the various photometric time series (TESS, ASAS-SN, and SNO V band). The amplitude of the R-band stellar variability is 3.0 mmag, which is small and suggests that TOI-1238 is a quiet star. Figure A.2 displays the TJO and SNO R-band photometry folded in phase with the periodicity of 40.3 d.

Our conclusion is that from all ASAS-SN, TJO, and SNO photometric data sets covering a total of 6.8 yr of regular monitoring, it is not possible to extract one single, characteristic rotation period for TOI-1238. Instead, we constrained the rotation period of this star to be in the interval between 30 and 45 d. This is likely explained by the presence of differential rotation. In the following subsection, we demonstrate that the spectroscopic activity indicators also point to a rotation period within this range.

Fig. 11 Window function of the CARMENES data (top panel) and the GLS periodograms of the CARMENES RVs and various spectroscopic indices (panels 2–10) shown in the time space between 0.5 and 1000 d. The horizontal and vertical lines and bands have the same meaning as in Fig. 10.

3.2.3 CARMENES activity indicators and radial velocities

We computed the GLS periodograms of some stellar activity indicators included in the CARMENES serval pipeline. All “spectroscopic” GLS periodograms are shown in panels 5–10 of Figs. 11 and A.4. The top panel of Fig. 11 displays the window function of the CARMENES observations, where the most significant peaks are located at 1 d and ~1 yr.

There are three Ca II IRT indices, one per atomic line; their corresponding GLS periodograms are depicted together in the fifth panel of Fig. 11. The dLW index periodogram is shown in the ninth panel of the same figure. The periodograms of the Ca II IRT and dLW all show peaks at ~40–45 d reaching FAP levels above the 0.1% (Ca II IRT) and the 10% (dLW) significance levels. In both cases, the peak is quite broad with a pedestal extending from ~30 through ~60 d for the dLW index (marked with a greenish band in all panels of Figs. 10 and 11). The 40 d feature of the Ca II IRT indices is asymmetric. In order to explore the presence of a signal at half of the stellar rotation period, we removed the 40 d signal (by using a sinusoid function) from the original time series (shown in the sixth panel of Fig. 11). Albeit with less significance, these residuals have a peak at ~19 d (about half of the main peak), a pattern that was already seen in the photometric LCs and which hints at the presence of active regions on different locations of TOI-1238’s surface.

The GLS periodogram of the Hα index (bottom panel of Fig. 11) shows no significant peak, which is consistent with TOI-1238 being an Hα inactive star. This is further corroborated by the fact that the Hα line is observed in absorption in TOI-1238 (Jeffers et al. 2018). The CRX times series shows a linear trend in the first 2.5 months of observations (Fig. 5), which is reflected as a long-term periodicity in the GLS periodogram of the seventh panel of Fig. 11. After removing this trend from the data, the resulting GLS periodogram of the CRX residuals (eighth panel of Fig. 11) has no significant peaks above any of the three defined FAP levels. We also investigated correlations between the measured CARMENES RVs and all the individual activity indices provided by the serval pipeline (including He I D₃, He I λ10 833 Å, Pa_β, the Na I D doublet, and the TiO and VO bands) using Pearson’s r and p coefficients to detect correlations and to access the significance of the correlation. We found that there are no strong, or even moderate correlations between the RVs and the activity indices. The lack of correlation or anticorrelation is also confirmed when considering the data of the first season only.

In conclusion, the CARMENES dLW and Ca II IRT indices suggest that the characteristic frequency for the stellar variations is around ~40–45 d, which agrees with the results obtained from the analysis of the coeval TJO and SNO LCs (Sect. 3.2.2). Interestingly, there is no significant peak at this periodicity in the GLS periodogram of the RVs illustrated in Fig. 5, but the RVs show a strong decreasing trend. We proceeded to remove this trend by fitting a sinusoidal function with a long period of ~800 d (see Sect. 4.2.3). The GLS periodogram of the RV residuals is illustrated in the third panel of Fig. 11, where the orbital periods of the two transiting planets become the most prominent peaks. Clearly, the Keplerian signals are present in our CARMENES RV data, thus allowing us to characterize the TOI-1238 planetary system (see next sections). The stellar activity signal expected at ~30–45 d, although present, it is not significant. Only when the planetary imprint is removed from the RV data, the signal at ~45 d gains power (fourth panel of Fig. 11) and has a broad structure with a width similar to that of the other spectroscopic indices. Because thisperiodicity coincides with that observed from the photometric data, we ascribe it to the rotation period of TOI-1238. To summarize, taking into account all photometric and spectroscopic indices, we conclude that the rotation period of TOI-1238 is P_rot = 40 ± 5 d.

4 TOI-1238 planetary system

4.1 TESS transit analysis

We applied the box least squares (BLS) periodogram (Kovács et al. 2002; Hartman & Bakos 2016) to the TESS time series data to verify that the announced planets are there and to search for additional transits that may have been missed by the TESS pipeline. Weemployed the BLS algorithm programmed in the astropy.timeseries python package and made the process iterative: once a transiting planet candidate is identified, it is masked out and the algorithm is run again to search for additional transiting planets, from the strongest to the weakest signals. The BLS periodograms of the PDCSAP fluxes of the four sectors computed for periodicities between 0.5 and 100 d are depicted in Fig. 12, although in the diagram we only show up to 20 d for clarity. The first transiting planet candidate to be detected in the TESS data has a BLS peak at 3.29 d, and the second candidate appears at 0.76 d. No other significant BLS peaks were identified after these two candidates were properly masked out from the LCs. The python package also delivered first estimates of the periastron passage (t₀) and the transit depth (t_depth) for the two candidates, all of which agree with the TESS alerts on TOI-1238. From now on, we name the 0.76 d and the 3.29 d transiting planets TOI-1238 b and TOI-1238 c, respectively.

An inspection of the bottom panels of Fig. 2 shows some structure in the TESS LCs, particularly in Sector 22. This may have an impact on the determination of the planetary parameters, particularly for shallow transits. We therefore flattened the TESS LCs using Gaussian process (GP) regressors with the celerite Mátern kernel (Foreman-Mackey et al. 2013)of the form $$\underset{ϵ\to 0}{\lim }{k}_{i,j}(\tau )={\sigma }^2\left(1+\frac{\sqrt{3}\tau }{\rho }\right)\exp \left(\frac{-\sqrt{3}\tau }{\rho }\right),$$(1)

where τ = |t_i − t_j| is the time-lag, σ is the amplitude of the GP, and ρ is the time/length-scale of the GP. All planetary transits were masked out. The fit was done using the juliet package, which also required thefollowing parameters: a dilution factor to account for possible contaminating sources in the photometric aperture, and a jitter to be added to the nominal flux error bars. We adopted a dilution factor of 1.0, that is, no contamination in the TESS photometry that may mimic a possible planetary transit (see Sect. 2.1 and Fig. 1). The priors on the GP σ and ρ parameters were based on log-uniform distributions with wide open intervals. We modeled the LCs of all four sectors finding that the models were essentially flat for sectors 14, 15, and 21 (Fig. 2). We used the GP model to “detrend” the data of sector 22. This flattened sector was added to the other three sectors for the characterization of the transiting planets.

In a next step, we modeled the planetary transits from the flattened photometric data to measure the orbital periods (P) and relative planet-to-star sizes (R_p ∕R_*), the time-of-transit center of the planet (t₀), the inclination of the planetary orbital planes (i), and the star–planet distances (a∕R_*). The julietpackage uses the batman package (Kreidberg 2015) to this end. The stellar limb-darkening coefficients (quadratic law) were parameterized following Kipping (2013). We also used the r₁ and r₂ parameterization (Espinoza et al. 2019) instead of determining the relative radii and the impact parameters (b) of the planets directly: r₁ and r₂ can vary between 0 and 1 and are defined to explore all the physically meaningful ranges for R_p ∕R_* and b. We defined a prior on the stellar density, ρ_⋆ (in kg m⁻³) instead of the scaled semimajor axis of the planets, a. In this way, a single value of ρ is defined for the system. We assumed that the two transiting planets candidates are in circular orbits. The median of the posterior distributions for each fitted parameter is given in the upper part of Table 4, while the bottom part presents the derived transit parameters. Quoted error bars correspond to 1σ uncertainties. The outer transiting planet candidate is twice the size of the inner planet candidate and is located only 2.6 times farther away from their parent star. Within the error bars, the two objects are in coplanar orbits. The TESS LCs folded in phase with orbital periods are illustrated in Fig. 13.

Fig. 12 BLS periodograms of the four TESS sectors. Top panel: detection of the 3.29 d transiting planet (strongest peak marked with an orange vertical line). Various harmonics are also identified with vertical thin red lines. Bottom panel: detection of the 0.76 d transiting planet (strongest peak) after masking out the previous signal. The only harmonic visible in the plot is marked with a thin red line.

4.2 CARMENES radial velocity analysis

4.2.1 Pre-whitening

The CARMENES RVs show a long-term linear trend whose origin is discussed in Sect. 4.2.2. After removing this trend, the GLS periodogram of Fig. 11 shows strong peaks at the expected planetary orbital periods. This confirms that the TESS transit-like features are planetary in nature.

To identify the presence of possible aliasing phenomena in the CARMENES RV data caused by the gaps in the time sampling of the observations (e.g., Dawson & Fabrycky 2010; Stock et al. 2020), we took into account the spectral window displayed in the top panel of Fig. 11. The strong peaks of the window function may introduce alias peaks in the RV periodogram at frequencies according to the expression f_alias = f_true ± mf_window, where m is an integer, f_true is the frequency identified in the RV periodogram and f_window the frequency from the window function (Deeming 1975). Typical aliases affecting ground-based observations are those associated with the year, synodic month, sidereal day, and solar day. In our spectroscopic window function, the highest peaks occur at ~1 d and at periodicities beyond 300 d or 1 yr. We checked that the latter does not introduce any peak in the RV periodogram that can be misinterpreted. However, the orbital periods of TOI-1238 b and c are, coincidentally, the ~1 d alias of eachother. This implies that, without the assistance of any other data free of this aliasing effect, using the RV time series alone will not deliver reliable masses for any of the two transiting planets. As an example, to produce the GLS periodogram of the fourth panel of Fig. 11, where the planetary signal is removed from the RV data using a simple sinusoidal function, we only subtracted the peak at 3.29 d and the 0.76 d was automatically subtracted, too. Fortunately, the combined TESS–CARMENES data analysis of Sect. 4.3 will overcome this issue.

We verified that the RV signals at 0.76 d and 3.29 d are stable and coherent over the entire observational time baseline by producing the stacked Bayesian generalized Lomb-Scargle (s-BGLS) periodogram (Mortier et al. 2015) shown in Fig. 14. For this, we employed the CARMENES RVs free of the long-term trend. The significance or probability of both planetary signals increases with time until a stability is reached at a certain number of observations. Then, the signals become narrower, which is a behavior expected for signals with a Keplerian origin. This provides further support for the exoplanetary nature of TOI-1238 b and c.

Fig. 13 TESS LCs folded in phase with the orbital periods of TOI-1238 b (top panel) and TOI-1238 c (bottom panel). Binned data are plotted as white circles. The best transit model for each planet is shown with a black line (associated uncertainty is illustrated with a whitish color).

Fig. 14 Evolution of the s-BGLS periodogram of the CARMENES RV data around 0.764 d (left) and 3.294 d (right). We used the RVdata after subtracting the long-term signal (see text).

4.2.2 Nature of the long-term RV signal

TOI-1238 is a K7–M0 star with a projected rotational velocity of ≤ 2 km s⁻¹, which is consistent with our inferred rotation period of 40 ± 5 d typical for field stars of this type. The decreasing, long-term trend of the CARMENES RVs displayed in the top panel of Fig. 5 may have two possible origins: a companion moving around the central star at a long orbital period or a magnetic activity cycle intrinsic to the star. Our spectroscopic data yield a lower limit on the period of such variation, P ≥ 600 d (twice the CARMENES time coverage), and a minimum RV amplitude of 70 m s⁻¹ (half the observed RV peak-to-peak variation).

The GLS periodogram of the combined ASAS-SN V and g LCs (covering more than 5 yr of observations) is not conclusive (the high dispersion of the photometry is an issue). We also combined all TJO and SNO photometry by applying suitable offsets to the individual LCs and thus covering about 300 d of continuous observations, but this time baseline is insufficient to sample one full cycle length of the long-term RV signal. The CARMENES CRX and Hα indices and all Ca II IRT indices after subtraction of the signal at the stellar rotation period appear to show a shallow slope with time, but it is not possible to relate these variations to the long-term trend seen in the RV time series. Also, the GALEX strong NUV emission suggests significant chromospheric activity, thus leaving open the possibility that the P ≥ 600 d RV feature could be due to a strong magnetic activity cycle of TOI-1238.

However, the large RV amplitude of the long-term variation (≈140 m s⁻¹ from peak to peak) suggests that the most likely scenario is that of the companion. With the exception of young stars (less than several hundred Myr), which show RV variations of tens to hundreds of m s⁻¹ due to rather strong and rapid stellar activity (Lillo-Box et al. 2020b; Suárez Mascareño et al. 2021), no main sequence star with slow rotation period appears to have such high RV variability caused by atmospheric changes. TOI-1238 is likely an old star based on its rotation period and Galactic kinematics. The age of TOI-1238 inferred from gyrochronology relations is likely >1 Gyr (see Fig. 13 of Barnes 2007). Furthermore, Hα is not seen in emission, and the photometric amplitude of the stellar variations in the TESS, TJO, and SNO data is on the order of a few to several mmag.

We also inspected the CARMENES NIR data. These RVs, despite their large error bars and dispersion, are useful for studyingthe large-amplitude, long-term trend. The downward slope of the RV time series is present in the NIR as well. A joint RV analysis of the CARMENES VIS and NIR data to fit this feature only (following the prescription described in Sect. 4.2.3) revealed that both RV curves share the same amplitude and period. The constant amplitude with wavelength is expected for a Keplerian signal, thus supporting the companion scenario.

Based on the spectroscopic constraints on the orbital period and RV amplitude of the putative outer companion of TOI-1238 and using Kepler’s laws, we imposed the following minimum values on the companion’s mass and orbital semimajor axis: $M\ge 2\sqrt{1-e^2}{M}_{\text{Jup}}$ (e is orbital eccentricity) and a ≥ 1.1 au. In what follows, we used all available photometric and astrometric observations to further constrain the companion’s properties, particularly its maximum mass. Our results are graphically summarized in Fig. 15, where the constraints derived from the various observables are illustrated in the form of exclusion areas (colored regions) in the mass versus semimajor axis diagram of the figure. The axis of mass goes from the minimum possible mass of the companion in a circular orbit, 2 M_Jup, up to the stellar mass of TOI-1238 (0.59 M_⊙). The axis of orbital separation starts at 1.1 au and extends to 1700 au, which is the value of the cross between the restriction imposed by Kepler’s laws to satisfy K ≥ 70 m s⁻¹ (plotted as orange triangle in Fig. 15) and another restriction described next.

From the Gemini high spatial resolution image (Fig. 4), TOI-1238 has no stellar companion with Δm ≈ 6 mag (contrast) at projected angular separations in the interval 0′′.2–1′′.2 (i.e., 14–85 au at the distance of the system). To explore closer separations, we used the precise Gaia photometry shown in Fig. 6, where it becomes evident that TOI-1238 is not overluminous with respect to the main sequence of early-M type stars: we thus set Δm = 0.05 mag as the maximum possible flux/luminosity contribution of the putative companion to the observed photometry of TOI-1238. Following the mass–magnitude relations of Pecaut & Mamajek (2013) and Cifuentes et al. (2020) valid for main sequence stars with solar metallicity, both the Gemini flux contrast and the Gaia maximum flux contribution translate into masses of ~0.16 M_⊙ and 0.18 M_⊙ (or spectral types in the interval M4–M5). That is, companions more massive than ~0.18 M_⊙ at a < 14 au and companions more massive than ~0.16 M_⊙ at a = 14–85 au can be discarded. These exclusion regions are marked with blue and green colors in Fig. 15.

Gaia imaging and astrometry can also be used to exclude resolved companions to TOI-1238. Gaia EDR3 data have an effective spatial resolution of 0′′.4 and a limiting G-band magnitude of 21 (or ≈28 au and M(G) = 16.7 at the distance of the system) in the sky area of TOI-1238 (Gaia Collaboration 2021). We explored all Gaia detections in a radius of 2′ around our star and none has the same proper motion and trigonometric parallax as TOI-1238, thus indicating that there are no proper motion companions identified down to the Gaia detectability limit, corresponding to a mass of ≈0.077 M_⊙ (close to thedefinition of the substellar mass limit for solar metallicity). This constraint is plotted as a pink exclusion region in Fig. 15.

Gaia EDR3 contains data spanning 34 months (1.4 yr) of observations out of a nominal five-year duration mission. TOI-1238 has a Gaia astrometric excess noise of 98 μas, which does not stand out for being very high compared to that of other Gaia objects of similar brightness in the same region of the sky, but it is not small either. Actually, this astrometric excess noise is detected with a σ =13, which is significant in terms of Gaia measurements according to Lindegren et al. (2021). Because the time coverage of Gaia EDR3 and our minimum limit on the orbital period of the putative external companion are alike, a measurable astrometric excess noise may hint at a relatively short orbital period on the order of a few to several years. We will need to wait for the completion of the Gaia mission before we can conclude on the nature of the companion to TOI-1238 using precise astrometric measurements. Also, additional RV measurements will contribute to find a final solution for the mass and other orbital parameters of the putative external companion to TOI-1238. This mass is expected to be contained within the white area of Fig. 15: the outer companion of TOI-1238 could be a very low-mass star with M ≤ 0.18M_⊙, a brown dwarf or a massive planet.

4.2.3 Radial velocity models

Despite the aliasing issue described in Sect. 4.2.1 that prevented us from determining the masses of the transiting planets from the RV data only, we modeled the CARMENES RV time series in order to constrain the most critical ingredients for a subsequent, more sophisticated, joint photometric and spectroscopic analysis of the TOI-1238 planetary system. This step was necessary tosave computing time. The models were produced with the juliet code that calls the radvel (Fulton et al. 2018) package to model Keplerian RV signals. The stellar activity signals were modeled by means of a GP with an exp-sin-squared kernel provided by george (Ambikasaran et al. 2015) python library, which is suited for learning periodic functions: $$k(\tau )={\sigma }_{\text{GP}}^2\exp \left(-{\alpha }_{\text{GP}}{\tau }^2-\Gamma {\sin }^2\left(\frac{\pi \tau }{P_{\text{rot}}}\right)\right),$$(2)

where σ_GP is the amplitude of the GP component given in m s⁻¹, Γ is the amplitude of the GP sine-squared component and is dimensionless, α is the inverse length-scale of the GP exponential component given in d⁻², P_rot the period of the GP quasi-periodic component given in days, and τ is the time lag.

We based the selection of the best model on the rules defined by Trotta (2008) based on the Bayesian model log-evidence, $\ln \mathcal{Z}$: if $\Delta \ln \mathcal{Z}\le 3$ the two models are indistinguishable and neither is preferred, while if $\Delta \ln \mathcal{Z}>3$ the model with the larger Bayesian log-evidence is favored. We performed four different approaches that are summarized in Table 5. All include one base model (BM) consisting of a RV offset and a RV jitter. All planetary orbits are assumed to be circular (with one exception indicated below). The other ingredients are as follows. The first is two Keplerian signals at 0.764 d (TOI-1238 b) and 3.294 d (TOI-1238 c) plus a linear trend to account for the RV long-term trend (BM+2pl+LT model). Thesecond is two Keplerian signals at 0.764 d (TOI-1238 b) and 3.294 d (TOI-1238 c) plus the periodic GP with an open P_rot prior value to account for the RV long-term trend (BM+2pl+GP model).

The third is three Keplerian signals (TOI-1238 b, TOI-1238c, and the external companion), that is, we assumed that the RV long-term trend is due to an external companion (BM+3pl model). This model is justified by the fact that the stellar rotation signal in the RV data appears to be lower in amplitude than the contribution of the two transiting planets (see Sect. 3.2.3). This model was explored with both a null eccentricity and a free eccentricity for planets TOI-1238 b and c. The computations indicate that an eccentricity near zero is preferred.

The final one is the same as above, with the addition of a GP model with P_rot in the interval 30–60 d to simulate the stellar activity due to stellar rotation (BM+3pl+GP model).

The resulting log-evidence for each model is provided in Table 5. Both the BM+3pl and BM+3pl+GP models are equally valid. We interpret this as follows: the former model is successful because the impact and amplitude of the stellar variability at timescales of 30–45 d in the RV time series are smaller than those of the planetary components. However, the BM+3pl+GP is physically more plausible because it includes all components observed in the data (two confirmed planets, a long-term signal, and stellar rotation). The GP model does not improve the fit in RVs significantly, and the increase in the number of free parameters is mathematically penalized in the log-evidence parameter. In any case, after the comparison of the models with the two and three Keplerians, the log-evidence criteria suggest that the RV long-term feature is better reproduced by a Keplerian function than by a periodic GP function or a single linear trend, which favors the presence of a third companion in the TOI-1238 system with relatively short orbital periods on the order of years.

Fig. 15 Constraints on the mass and semimajor axis of the external companion to TOI-1238 based on photometric, astrometric,and spectroscopic observables. Colored areas represent regions excluded by the various observations (see legend and text). Possible masses and orbital separations of this companion are contained within the white area. The limits of theordinate axis are given by the possible minimum mass of the companion (in a circular orbit) and the stellar mass, while the limits of the abscissa axis are determined by the minimum and maximum orbital separations.

4.3 Masses of the transiting planets

To determine the true masses of TOI-1238 b and c, we performed a combined photometric and spectroscopic analysis using the TESS and CARMENES VIS data. Our final model consists of two transiting planets at 0.76 and 3.29 d and a long-term ≥ 400 d signal, all three of which are modeled by non-eccentric Keplerian orbits, plus the stellar rotation component at 20–80 d simulated by the exp-sine-squared GP kernel of Eq. (2). On the one hand, we used the results obtained from the transit-only analysis (Sect. 4.1) to define normal priors on the orbital period and time of periastron passage of the transiting planets with Gaussian distributions centered at the median values of the posteriors. This is fully justified because these parameters are mainly constrained by the LCs with the RV data adding little information (Kemmer et al. 2020). On the other hand, the RV amplitudes, K, are fit by adopting a prior with a uniform distribution for each Keplerian signal. We also fit a jitter for the TESS and CARMENES data. All priors are summarized in Table B.1.

As illustrated by the median posterior parameters of our combined fit presented in Table 6, the derived GP period is ${47.8}_{-8.7}^{+11.5}$ d, which agrees within the quoted uncertainties with the stellar rotation period determined from independent LC data and from the CARMENES spectroscopic indices. The combined fit also confirmed that TOI-1238 displays a weak activity at optical wavelengths and timescales typical of the stellar rotation period. Figure 16 shows the four TESS sectors and the location of all transits of TOI-1238 b and c captured by the observations. On three different occasions (BJD = 2 458 713.9 – sector 15, BJD = 2 458 895.1 – sector 21, and BJD = 2 458 905.0 – sector 22), the transits of TOI-1238 b and c overlap. The LCs folded in phase with the orbital periods of the transiting planets per sector are shown in Fig. A.5. For completeness, the corner plot depicting all the posterior distributions of the planetary parameters as obtained from the joint fit is shown in Fig. A.6.

The resulting RV model is depicted in the top panels of Fig. 17, whereas the RV curves folded in phase for the two transiting planets are shown in the bottom panels of the figure. The rms of the RV residuals (that is, observed RVs minus the best fit) is 4.66 m s⁻¹, which is slightly higher than the mean value of the CARMENES VIS RV errors (see Sect. 2.6). This suggests that no other component is detectable at the quality of our data. With RV amplitudes of 3.74${}_{-0.99}^{+1.03}$ m s⁻¹ (TOI-1238 b) and 5.10${}_{-1.06}^{+1.02}$ m s⁻¹ (TOI-1238 c), these transiting planets have true masses of 3.76${}_{-1.07}^{+1.15}$ M_⊕ and 8.32${}_{-1.88}^{+1.90}$ M_⊕, respectively, with a significance of 3.7–5 σ. They are moving around their parent star in orbits that are 28 and 11 times smaller than the Sun–Mercury orbital size. TOI-1238 b and TOI-1238 c are a small super-Earth and a mini-Neptune, each one located on either side of the radius gap for planetary systems around M-type stars (Cloutier & Menou 2020). As for the possible third companion, we can only impose a minimum mass of $M\ge 2\sqrt{1-e^2}{M}_{\text{Jup}}$ as explained above.

We also explored whether the two inner planets have eccentric orbits by leaving this parameter free in a new combined photometricand spectroscopic analysis while maintaining a zero eccentricity for the outer companion. The results revealed that TOI-1238 b and c are very likely in circular orbits and we derived an upper limit on the orbital eccentricities of e = 0.25 and 0.15 at the 1σ level for planets b and c, respectively.

In summary, TOI 1238 is the host of a super-Earth planet (TOI-1238 b) with a very short orbital period (0.76 d), a mini-Neptune planet (TOI-1238 c) also at a close-in orbit (3.29 d), and a likely more massive and distant companion with an orbital period ≥600 d. We also derived the bulk densities of the transiting planets ($\rho ={11.7}_{-3.4}^{+4.2}$ and 4.9${}_{-1.8}^{+2.5}$ g cm⁻³ for planets b and c, respectively), which are very distinct despite their proximity to the parent star, thus suggesting dissimilar atmosphericconditions.

Fig. 16 Location of the transits of TOI-1238 b (weakest, black line) and TOI-1238 c (strongest, orange line) across all four TESS sectors.

Fig. 17 TOI-1238 CARMENES RVs (blue dots) and the best model (black line) from the combined photometric and spectroscopic fit. The top panels shows the entire RV time series, where the short-period planets and the long-term feature are clearly visible. The inset displays an enlargement of the first part of the data for a zoomed-in view of the two transiting planets’ contribution to the modeled RV curve. The middle panel illustrates the RV residuals of the fit. The bottom panels present the CARMENES RVs folded in phase with the orbital period of the inner Keplerian components (from left toright: TOI-1238 b and TOI-1238 c). Binned data are plotted as open black symbols. All error bars include the quoted CARMENES uncertainties and the RV jitter as obtained from the model added in quadrature.

4.4 Aliasing

The 1 d aliasing issue between the 0.764 d and 3.294 d planets in the RV data may cause spectral leakage making the RV amplitude of one signal to be affected by the other. The work by Dawson & Fabrycky (2010) established that the aliases of an RV signal also includes phase information that depends on the true underlying signal. This information can be used to disentangle the alias signal of the short orbital period planet from the true signal of the second planet and vice versa if the phases are not equal by coincidence. The success in disengaging the two alias signals also depends on the time of the day at which the measurements were performed, that is, the time sampling. In our case we obtained a few RV measurements on different occasions during the same nights, which helped break this degeneracy. The additional phase information obtained when performing the combined photometric and spectroscopic analysis delineated in Sect. 4.3 contributed to distinguish the two planetary signals. Another argument indicating that the phases of TOI-1238 b and c differ enough and that they are well constrained to distinguish the two signals comes from the relation between both RV amplitudes (K). If no strong correlation between K_b and K_c is found, the two signals can be safely untangled.

To confirm that the degeneracy is broken in this way, we performed the following test: we removed the planetary signals from the RV curve using the results of Sect. 4.3 and injected two fake signals with the same planetary orbital periods and phases but different RV amplitudes. We used the same observing CARMENES epochs to produce the fake RV time series. This provides a useful test for checking potential correlations between K_b and K_c. The exact same method (including the priors) described in Sect. 4.3 was then applied to the combined TESS photometry and fake RVs. We managed to recover the injected amplitudes within the error bars in all tests, thus supporting the combined analysis for deriving reliable planetary masses of the TOI-1238 system. Furthermore, the box displaying the RV amplitude of TOI-1238 b versus the RV amplitude of TOI-1238 c in the corner plot of Fig. A.6 shows a nicely uncorrelated relation resulting from the photometric and spectroscopic joint analysis of the real data. The Pearson’s r and Spearman’s ρ coefficients obtained for the planet amplitude relation have a value of − 0.4, and the p coefficient is 0.0 for both correlations. Thanks to the precise phase and period information from the joint analysis described in Sect. 4.3 and an adequate time sampling of the RV observations, we managed to break down the alias issue for planets TOI-1238 b and c. This resulted in a robust and precise solution for the RV amplitudes of the two transiting planets, and therefore, the masses of TOI-1238 b and c are uniquely measured.

5 Discussion

5.1 Planetary system stability

We investigated whether the planetary system is dynamically stable using the angular momentum deficit (AMD) stability criterion (Laskar & Petit 2017). The AMD can be interpreted as a measure of the excitation of the orbits that limits the close encounters among the planets and ensures long-term stability. The main ingredients are the semimajor axes, the masses and the orbital eccentricities. The result of the analysis rejected high eccentric planetary orbits yielding stable AMD solutions only for those posterior distributions with low eccentricity values: e ≤ 0.2 for the TOI-1238 b and c and e ≤ 0.1 for the outer companion, in agreement with our findings.

5.2 Outermost companion to TOI-1238

The putative, external, massive companion of TOI-1238 has an estimated radius between ~1.2 and 1.7 R_Jup for masses in the interval 2 M_Jup and 0.18 M_⊙ (Baraffe et al. 2003; Chen & Kipping 2017; Cifuentes et al. 2020). Despite its large size, the probability for this possible companion to cross the disk of its parent star is very low, on the order of 0.2 ± 0.2%, and it is thus compatible with it being a non-transiting companion. However, in the case where all companions have coplanar orbits, the transit probability of the external companion significantly increases up to ~20%. Under the low probability hypothesis of a supposed transit with an inclination of 90 deg, the depth of the transit would be 43456 ± 1242 ppm (adopting 2 R_Jup for the external companion). This is an extremely intense phenomenon that would be visible even in poor-quality LCs. Although statistically not very abundant, planetary systems with unresolved, cold Jupiter-mass companions orbiting G, K, and M dwarfs are known to exist (e.g., Rivera et al. 2010; Zhu & Wu 2018; Morales et al. 2019; Bakos et al. 2020; Trifonov et al. 2020; Lindor et al. 2021; Kim et al. 2021). If the external companion of TOI-1238 is a substellar object (planet or brown dwarf), TOI-1238 would add to the increasing number of these systems. What is more peculiar about the TOI-1238 system is, however, is that it also contains warm super-Earth and mini-Neptune planets in close-in orbits because massive companions are expected to influence their environment due to their large masses. Planetary systems with the architecture similar to that of TOI-1238 must be described by the theories of planet formation and dynamical evolution. Using the Generation III Bern Model of planet formation and evolution (Emsenhuber et al. 2021), where most planets are formed via the core accretion mechanism and moved to their final location in protoplanetary disks, Schlecker et al. (2021) investigated the mutual interaction between cold giants and warm super-Earths orbiting solar-type stars. Their findings indicate that there is a positive correlation between the occurrences of inner super-Earths and cold Jupiters, which may indicate that these systems are born from intermediate-mass disks and that the efficiency of the planetary inward migration of the giants is somehow overestimated in planetary formation and evolution models. The presence of this massive outer companion is also compatible with formation by disk gravitational instability leading to the collapse of dense regions in the protoplanetary disk (Boss 1997; Kratter & Lodato 2016, and references therein).

5.3 TOI-1238 b and c

Kopparapu et al. (2014) calculated a conservative estimation of the inner habitable zones (HZs) around stars with stellar effective temperatures in the range 2600–7200 K for planetary masses between 0.1 and 5 M_⊕. According to these authors, the inner edge of the HZ around TOI-1238 is located at 0.28 au with an effective stellar flux incidence of S_eff = 1.02S_⊙. TOI-1238 b and c lie between the star and the inner boundary of the HZ. The theoretical equilibrium temperature (T_eq) of these two planets was derived by using the Stefan–Boltzmann equation, the stellar parameters given in Table 2, and two very different values of the planetary albedo (A = 0.65 and 0.0 for the high- and non-reflectance cases). The resulting T_eq ranges are 965–1300 K and 590–800 K for TOI-1238 b and c, respectively (Table 6). These temperature intervals include the uncertainties in the stellar T_eff and semimajor axis of the planets.

Planets with orbital periods of less than one day are known as ultra-short-period planets (USPs; Sanchis-Ojeda et al. 2014; Winn et al. 2018). The number of known USPs is low, and about 30% of them have been found around M-type stars. To place the TOI-1238 planetary system in context with other known systems, Fig. 18 shows the planetary mass versus the orbital period diagram for all planets discovered by any observing method (from the NASA Exoplanet Archive¹⁰). All of these systems orbit around F-, G-, K-, and M-type stars. TOI-1238 b is an USP low-mass planet in a hierarchical planetary system where the inner-most planet is the least massive and the outermost likely companion is the most massive. As shown in the left panel of Fig. 19, TOI-1238 b is one of the few super-Earths orbiting M dwarfs with the highest stellar irradiation index and highest T_eq (in the interval 1000–1258 K; see Table 6). Therefore, the day side of TOI-1238 b is presumably very hot.

If TOI-1238 b had an atmosphere, it would be subject to intense erosion (Uzsoy et al. 2020; Millholland & Spalding 2020), making it a potentially interesting candidate for planetary atmospheric studies. The stellar extended UV (XUV) irradiation ionizes H atoms, giving rise to atmospheric inflation in the planet in the first place, and to mass loss afterward (Sanz-Forcada et al. 2011, and references therein). In order to test this possibility, we would need to know the XUV irradiation flux from the star. This flux can be calculated using coronal models, but no X-ray data of TOI-1238 are available in the literature. Instead, we can estimate the current X-ray emission of the star by using the relation between X-rays and the stellar rotational period, as displayed in Fig. 2 of Wright et al. (2011). According to this relation, using P_rot = 40  ± 5 d, we would expect an X-ray luminosity of L_X = 3.3 × 10⁺²⁸ erg s⁻¹ in the 5–100 Å spectral range. The stellar extreme ultraviolet (EUV; 100–920 Å) luminosity can be calculated using the scaling law calculated by Sanz-Forcada et al. (2011), resulting in L_EUV = 2.0 × 10⁺²⁹ erg s⁻¹. Using these quantities, and the planetary system parameters from Table 6, we can calculate the mass loss rate in the energy limited regime, as in Eq. (7) of Sanz-Forcada et al. (2011), resulting in 3.2 M_⊕ Gyr⁻¹ for TOI-1238 b and 0.9 M_⊕ Gyr⁻¹ for TOI-1238 c. This formula assumes a planet atmosphere dominated by H, which is probably not the case presently for TOI-1238 b, given its high density. We can instead estimate what was the expected mass loss rate at some point in the past, assuming that TOI-1238 planets weremore massive with an H-dominated atmosphere. Typical densities of planets with heavy XUV irradiation can be as low as 0.3 g cm⁻³. Even assuming that the star had the same XUV emission in the past, we would have a mass loss rate of 125 M_⊕ Gyr⁻¹ for TOI-1238 b, and 15 M_⊕ Gyr⁻¹ for TOI-1238 c, assuming the density of the planets reported in Table 6. The stellar XUV emission decreases as their rotation rate slows down. Therefore, we adopted a conservative assumption, yet we inferred planet photo-evaporation in the past was very strong, likely leading to the loss of a significant fraction of the H atmosphere. This might explain the observed planet density, inconsistent with a gas-rich atmosphere. This scenario is also consistent with an early migration of the planets through the protoplanetary disk.

Following the metrics proposed by Kempton et al. (2018) to identify the TESS transiting planets most amenable for atmosphericcharacterization via transmission and/or emission spectroscopy with the James Webb Space Telescope (JWST), we derived TSM = 2.98 and 33.80, and ESM = 4.63 and 5.12 for TOI-1238 b and c, respectively, where TSM and ESM stand for the metrics for transmission and emission spectroscopy. These factors depend on the stellar brightness, the host star’s temperature, and the planetary properties, including the planet equilibrium and day-side temperatures. As a rule of thumb, for a given host star, the greater radius and temperature of the planet, the larger the TSM and ESM factors are, and the higher the probability of detecting the planetary atmosphere is. Both TOI-1238 b and c have TSM and ESM factors below the threshold metric values defined for their respective planetary categories (Earth-size and small sub-Neptune) by Kempton et al. (2018), thus indicating that the study of their atmospheres, if they had any, is not an easy task even for the JWST.

Various studies have revealed that the radius distribution of planets slightly larger than the Earth is bimodal (Fulton et al. 2017; Van Eylen et al. 2018; Petigura 2020; Cloutier et al. 2020). Super-Earth planets of up to ~1.5 R_⊕ are relatively common, as are mini-Neptunes in the range of ~2–4 R_⊕. But there appears to be a deficit of planets in between these sizes. The location of the radius gap depends on the orbital period and/or the planet insolation and/or the spectral type of the parent stars (Van Eylen et al. 2018, 2021; Cloutier et al. 2020). The right panel of Fig. 19 illustrates the planet radius versus orbital period for all known planets from the NASA Archive with parameters measured to better than 30%. We marked the location of the planetary radius valley for solar-type stars and M-type stars following Van Eylen et al. (2021). TOI-1238 b and c have very different bulk densities and are located at the two sides of the planetary radius gap of close-in planets, which is believed to account for the distinction between small rocky planets and bigger planets with volatile rich envelopes (Dressing & Charbonneau 2015; Fulton et al. 2017).

Using the same subsample of well characterized planetary systems as in Fig. 19, we produced the mass-radius diagram shown in the left panel of Fig. 20, where we overplotted various theoretical models planets with multilayer interior composition (Zeng et al. 2016). The comparison of the loci of TOI-1238 b and c with these models yields that the inner, hotter planet is mainly made of iron (thus, very likely solid) whereas the composition of planet c is compatible with 50% Earth-like rocky core and 50% H₂O-gaseous layer. TOI-1238 c has a planetary density that resembles those of Venus and the Earth. Interestingly, TOI-1238 b appears to be denser than the Earth making it one of the densest planets known to date, particularly among the subsample of planets orbiting M dwarfs. This property is more apparent in the right panel of Fig. 20, where planetary radius is shown as a function of density. It is worth noting that the number of very dense close-in planets (ρ > 10 g cm⁻³) is quite small as compared to the number of planets with Earth-like densities. In their theoretical study, Schlecker et al. (2021) predicted that systems with close-in, high-density super-Earths are more likely to host an outer giant planet. TOI-1238 matches this theoretical prediction, although we caution that these models are based on several poorly constrained assumptions such as efficient planetary core formation via planetesimal accretion and a simplified disk chemistry.

Fig. 18 Mass versus orbital period diagram for all discovered planets around F-, G-, and K-type (gray dots) and M-type (blue dots) stars from the NASA Exoplanet Archive. The red star symbols represent TOI-1238 b and c. The shaded area on the left corresponds to the location of USPs. The red shaded area on the right represents the possible locations of the outercompanion to TOI-1238 (circular orbits). Venus, the Earth, Mars, Jupiter, and Uranus are also shown for comparison purposes.

Fig. 19 Radius of all discovered planets around F-, G-, and K-type stars (gray dots) and M-type stars (blue dots) with planetary mass and radius error bars of 30% or smaller, plotted as a function of planetary insolation (left panel) and orbital period (right panel). These data were extracted from the NASA Exoplanetary Archive. The two transiting planets of the TOI-1238 system are marked with a red star symbol. In the right panel, the dashed gray line indicates the radius valley for planets around solar-type stars, whereas the dashed blue line shows the location of the radius valley of planetary systems with M-type parent stars.

Fig. 20 Planetary radius shown as a function of mass (left panel) and density (right panel) for all planets with measured error bars smaller than 30%. Gray dots stand for F-, G-, and K-type primary stars, and blue dots represent planets orbiting M-type stars. The two transiting planets of TOI-1238, with radius and mass uncertainties of ~5% and ~25%, respectively,are shown with red stars. The various planetary composition models of Zeng et al. (2016) are shown with solid colored lines.

6 Summary

We have presented the discovery and validation of two transiting planets, TOI-1238 b and TOI-1238 c, orbiting a nearby K7–M0 star using photometric TESS and spectroscopic CARMENES data. The RVs also revealed the presence of an outer, more massive companion with a long-period orbit. The intrinsic stellar variability was analyzed using the spectroscopic activity indicators provided by CARMENES and the optical photometric monitoring of the ASAS-SN, TJO, SNO, and TESS data. Our findings suggest that TOI-1238 is a rather quiet star with a likely rotation period of 40±5 d.

The joint analysis of the TESS and CARMENES data yields orbital periods of ${0.764597}_{-0.000011}^{+0.000013}$ d and ${3.294736}_{-0.000036}^{+0.000034}$ d, masses of 3.76${}_{-1.07}^{+1.15}$ M_⊕ and 8.32${}_{-1.88}^{+1.90}$ M_⊕, and radii of ${1.21}_{-0.10}^{+0.11}$ R_⊕ and ${2.11}_{-0.14}^{+0.14}$ R_⊕ for the transiting planets TOI-1238 b and TOI-1238 c, respectively, both of which likely have nearly circular orbits with eccentricity below 0.2. From our data, we imposed the following constraints for the outer, massive companion: an orbital period longer than 600 d and an RV amplitude above 70 m s⁻¹, which translates into a minimum mass of $M\ge 2\sqrt{1-e^2}$ M_Jup and orbital separations ≥1.1 au from its parent star.

The TOI-1238 system contains all the fascinating ingredients of low-mass planetary systems, having two close-in low-mass planets at both sides of the planetary radius valley and a likely massive companion on an outer orbit. One of the inner planets is a very dense super-Earth that has very likely suffered strong irradiation and erosion. The characteristics of the TOI-1238 planetary system offer a compelling laboratory to test innovative models for explaining theories of planet formation and evolution.

Acknowledgements

CARMENES is an instrument for the Centro Astronómico Hispano-Alemán de Calar Alto (CAHA, Almería, Spain). CARMENES is funded by the German Max-Planck- Gesellschaft (MPG), the Spanish Consejo Superior de Investigaciones Científicas (CSIC), the European Union through FEDER/ERF funds, and the members of the CARMENES Consortium (Max-Planck-Institut für Astronomie, Instituto de Astrofísica de Andalucía, Landessternwarte Koönigstuhl, Institut de Ciències de l’Espai, Insitut für Astrophysik Göttingen, Universidad, Complutense de Madrid, Thüringer Landessternwarte Tautenburg, Instituto de Astrofísica de Canarias, Hamburger Sternwarte, Centro de Astrobiología and Centro Astronómico Hispano-Alemán), with additional contributions by the Spanish Ministry of Economy, the state of Baden-Wüttemberg, the German Science Foundation (DFG), the Klaus Tschira Foundation (KTS), and by the Junta de Andalucía. This work was based on data from the CARMENES data archive at CAB (CSIC-INTA). 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. Funding for the TESS mission is provided by NASA’s Science Mission Directorate. This paper includes data collected by the TESS mission that are publicly available from the Mikulski Archive for Space Telescopes. 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. Some of the observations in the paper made use of the High-Resolution Imaging instrument ‘Alopeke obtained under Gemini LLP Proposal Number: GN/S-2021A-LP-105. ‘Alopeke was funded by the NASA Exoplanet Exploration Program and built at the NASA Ames Research Center by Steve B. Howell, Nic Scott, Elliott P. Horch, and Emmett Quigley. ‘Alopeke was mounted on the Gemini North (and/or South) telescope of the international Gemini Observatory, a program of NSF’s OIR Lab, which is managed bythe Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. on behalf of the Gemini partnership: the National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). Data were partly collected with the 90 cm telescope at the Sierra Nevada Observatory (SNO) operated by the Instituto de Astrofí fica de Andalucí a (IAA-CSIC). We acknowledge the telescope operators from Observatori Astronómic del Montsec, Sierra Nevada Observatory, and Centro Astronómico Hispano-Alemán de Calar Alto (CAHA). E.G.A., M.R.Z.O., J.A.C., J.S.F, and D.M. acknowledge financial support from the Spanish Ministry of Science and Innovation through project PID2019-109522GBC5[1:4]. E.G.A also acknowledges support from AEI Project No. MDM-2017-0737 Unidad de Excelencia “María de Maeztu” – Centro de Astrobiología (CSIC-INTA). V.M.P. acknowledges financial support from NASA through grant NNX17AG24G. S.V.J. acknowledges the support of the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets (JE 701/5-1)”. M.J.L.-G., E.R., C.R.-L., and P.J.A. acknowledge financial support from the Agencia Estatal de Investigación of the Ministerio de Ciencia e Innovación through projects PID2019-109522GB-C52, PID2019-107061GB-C64, PID2019-110689RB-100 and the Centre of Excellence Severo Ochoa Instituto de Astrofí sica de Andalucí a (SEV-2017-0709). G.M. received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 895525. S.S. and S.R. acknowledgesupport by the DFG Research Unit FOR 2544 Blue Planets around Red Stars, project no. RE 2694/4-1.

Appendix A Figures

Fig. A.1 SNO R-band image of TOI-1238 (in the center). The reference stars used to obtain the differential photometry are marked with circles. The field of view is 13.2′×13.2′.

Fig. A.2 TJO and SNO R-band LCs folded in phase with the period Prot = 40.3 d (top panel). Residuals (observations minus model) are shown in the bottom panel. The white dots correspond to the binned data, where the x-error bars illustrate the interval of the bins and the y-error bars take into account the uncertainties of the binned data.

Fig. A.3 GLS periodograms of the ASAS-SN V and g, TJO R, and SNO V and R LCs in the frequency space between 2 and 0.001 d−1. In all panels, the horizontal dashed lines indicate FAP levels of 10% (blue), 1% (orange), and 0.1% (green). The orbital periods of the two transiting planets are marked with vertical blue lines. The greenish band indicates the region where most of the spectroscopic activity indicators have their highest GLS peaks. In the fourth panel, the vertical dashed line indicates the frequency of the signal removed from the data (see text).

Fig. A.4 Window function of the CARMENES data (top panel) and the GLS periodograms of the CARMENES RVs and various spectroscopic indices (panels 2–10) in the frequency space between 2 and 0.001 d−1. The horizontal and vertical lines and bands have the same meanings as in Figs. 10 and 11.

Fig. A.5 TESS LCs (blue points) folded in phase with the orbital periods of the transiting planets per sector. The best joint fit solution is plotted as a black line. The white dots correspond to the binned photometric data. The x axis represents the time computed from the mid-transit times as derived from the best joint fit.

Fig. A.6 Posterior distributions of some of the fitted planetary parameters of the TOI-1238 system as obtained from the combined photometric and spectroscopic fit. The vertical dashed lines indicate the 16, 50, and 84 % quantiles that were used to define the optimal values and their associated 1σ uncertainty. The blue line stands for the median values of each fitted parameter.

Appendix B Tables

References

All Tables

All Figures

Fig. 1 TPFs of TOI-1238 (cross symbol) for TESS sectors 14, 15, 21, and 22 (left to right). The electron counts are color-coded. The TESS optimal photometric aperture per sector used to obtain the SAP fluxes is marked with red squares. The Gaia DR2 objects with G-band magnitudes down to 6 mag fainter than TOI-1238 are labeled with numbers (TOI-1238 corresponds to number 1), and their scaled brightness, based on Gaia magnitudes, is shown by red circles of different sizes (see figure inset). The pixel scale is 21 arcsec pixel−1.

In the text

Fig. 2 TESS LCs. Top and middle panels: SAP and PDCSAP LCs of TOI-1238 as obtained by TESS in sectors 14, 15, 21, and 22 and processed by SPOC (left to right). Bottom panel: PDCSAP LCs with the detrending GP models overplotted (black line). The GP models were used to remove any residual instrumental effect (after masking the position of all transits). As illustrated, only the data of sector 22 required this additional correction to produce a flat LC.

In the text

	Fig. 3 Photometric time series of TOI-1238 using ASAS-SN (aperture photometry), TJO (differential photometry), and SNO (differential photometry) data.
In the text

	Fig. 4 Relative 5σ detection limit magnitude as a function of the angular separation for Gemini (‘Alopeke) at 562 nm (blue line) and 832 nm (red line). The inset figure shows the reconstructed ‘Alopeke high-contrast, high-resolution speckle image at 832 nm.
In the text

	Fig. 5 CARMENES VIS RV time series (blue dots) and various spectroscopic activity indicators (green dots) of TOI-1238.
In the text

	Fig. 6 Location of TOI-1238 (red circle) in the Gaia color-magnitude diagram. The young stellar sequences of the β Pic (blue) andTucana-Horologium (green) moving groups and the Pleiades star cluster (purple) are taken from Luhman (2018) and the main sequence of field stars (black) from Cifuentes et al. (2020). The gray area represents the dispersion observed among stars on the main sequence.
In the text

	Fig. 7 Photometric spectral energy distribution of TOI-1238 (red dots) and photometry of the best BT-Settl model (Teff = 3900 ± 100 K, log g = 5.5 ± 0.5 dex, and [Fe/H] = 0.0–0.5 dex, Allard et al. 2012). The reddest data point (W4) is an upper limit.
In the text

	Fig. 8 TESS photometry. Top panels: TOI-1238 TESS SAP photometric fluxes for the four different sectors (each half sector is shown with a different color) normalized to a common reference by fitting sinusoidal functions. The best model fit is plotted asthe solid black line, and the binned photometric data are indicated with gray dots. Bottom panels: residuals shown as a function of time.
In the text

	Fig. 9 Posterior distributions of the sinusoidal fit to the TESS SAP fluxes. Vertical dashed lines correspond to the 16 and 84% quantiles (1σ uncertainty). The red line stands for the mean values of the distributions.
In the text

Fig. 10 GLS periodograms of the ASAS-SN V and g, TJO R, and SNO V and R LCs shown in the time space between 0.5 and 1000 d. In all panels, the horizontal dashed lines indicate FAP levels of 10 (blue), 1 (orange), and 0.1% (green). The orbital periods of the two transiting planets are marked with vertical blue lines. The highest peak of the CARMENES RV GLS periodogram is shown with a violet line. The greenish band indicates the region where most of the spectroscopic activity indicators have their highest GLS peaks. In the fourth panel, the vertical dashed line indicates the periodicity of the signal removed from the data (see text).

In the text

	Fig. 11 Window function of the CARMENES data (top panel) and the GLS periodograms of the CARMENES RVs and various spectroscopic indices (panels 2–10) shown in the time space between 0.5 and 1000 d. The horizontal and vertical lines and bands have the same meaning as in Fig. 10.
In the text

Fig. 12 BLS periodograms of the four TESS sectors. Top panel: detection of the 3.29 d transiting planet (strongest peak marked with an orange vertical line). Various harmonics are also identified with vertical thin red lines. Bottom panel: detection of the 0.76 d transiting planet (strongest peak) after masking out the previous signal. The only harmonic visible in the plot is marked with a thin red line.

In the text

	Fig. 13 TESS LCs folded in phase with the orbital periods of TOI-1238 b (top panel) and TOI-1238 c (bottom panel). Binned data are plotted as white circles. The best transit model for each planet is shown with a black line (associated uncertainty is illustrated with a whitish color).
In the text

	Fig. 14 Evolution of the s-BGLS periodogram of the CARMENES RV data around 0.764 d (left) and 3.294 d (right). We used the RVdata after subtracting the long-term signal (see text).
In the text

Fig. 15 Constraints on the mass and semimajor axis of the external companion to TOI-1238 based on photometric, astrometric,and spectroscopic observables. Colored areas represent regions excluded by the various observations (see legend and text). Possible masses and orbital separations of this companion are contained within the white area. The limits of theordinate axis are given by the possible minimum mass of the companion (in a circular orbit) and the stellar mass, while the limits of the abscissa axis are determined by the minimum and maximum orbital separations.

In the text

	Fig. 16 Location of the transits of TOI-1238 b (weakest, black line) and TOI-1238 c (strongest, orange line) across all four TESS sectors.
In the text

Fig. 17 TOI-1238 CARMENES RVs (blue dots) and the best model (black line) from the combined photometric and spectroscopic fit. The top panels shows the entire RV time series, where the short-period planets and the long-term feature are clearly visible. The inset displays an enlargement of the first part of the data for a zoomed-in view of the two transiting planets’ contribution to the modeled RV curve. The middle panel illustrates the RV residuals of the fit. The bottom panels present the CARMENES RVs folded in phase with the orbital period of the inner Keplerian components (from left toright: TOI-1238 b and TOI-1238 c). Binned data are plotted as open black symbols. All error bars include the quoted CARMENES uncertainties and the RV jitter as obtained from the model added in quadrature.

In the text

Fig. 18 Mass versus orbital period diagram for all discovered planets around F-, G-, and K-type (gray dots) and M-type (blue dots) stars from the NASA Exoplanet Archive. The red star symbols represent TOI-1238 b and c. The shaded area on the left corresponds to the location of USPs. The red shaded area on the right represents the possible locations of the outercompanion to TOI-1238 (circular orbits). Venus, the Earth, Mars, Jupiter, and Uranus are also shown for comparison purposes.

In the text

Fig. 19 Radius of all discovered planets around F-, G-, and K-type stars (gray dots) and M-type stars (blue dots) with planetary mass and radius error bars of 30% or smaller, plotted as a function of planetary insolation (left panel) and orbital period (right panel). These data were extracted from the NASA Exoplanetary Archive. The two transiting planets of the TOI-1238 system are marked with a red star symbol. In the right panel, the dashed gray line indicates the radius valley for planets around solar-type stars, whereas the dashed blue line shows the location of the radius valley of planetary systems with M-type parent stars.

In the text

Fig. 20 Planetary radius shown as a function of mass (left panel) and density (right panel) for all planets with measured error bars smaller than 30%. Gray dots stand for F-, G-, and K-type primary stars, and blue dots represent planets orbiting M-type stars. The two transiting planets of TOI-1238, with radius and mass uncertainties of ~5% and ~25%, respectively,are shown with red stars. The various planetary composition models of Zeng et al. (2016) are shown with solid colored lines.

In the text

	Fig. A.1 SNO R-band image of TOI-1238 (in the center). The reference stars used to obtain the differential photometry are marked with circles. The field of view is 13.2′×13.2′.
In the text

	Fig. A.2 TJO and SNO R-band LCs folded in phase with the period Prot = 40.3 d (top panel). Residuals (observations minus model) are shown in the bottom panel. The white dots correspond to the binned data, where the x-error bars illustrate the interval of the bins and the y-error bars take into account the uncertainties of the binned data.
In the text

Fig. A.3 GLS periodograms of the ASAS-SN V and g, TJO R, and SNO V and R LCs in the frequency space between 2 and 0.001 d−1. In all panels, the horizontal dashed lines indicate FAP levels of 10% (blue), 1% (orange), and 0.1% (green). The orbital periods of the two transiting planets are marked with vertical blue lines. The greenish band indicates the region where most of the spectroscopic activity indicators have their highest GLS peaks. In the fourth panel, the vertical dashed line indicates the frequency of the signal removed from the data (see text).

In the text

	Fig. A.4 Window function of the CARMENES data (top panel) and the GLS periodograms of the CARMENES RVs and various spectroscopic indices (panels 2–10) in the frequency space between 2 and 0.001 d−1. The horizontal and vertical lines and bands have the same meanings as in Figs. 10 and 11.
In the text

	Fig. A.5 TESS LCs (blue points) folded in phase with the orbital periods of the transiting planets per sector. The best joint fit solution is plotted as a black line. The white dots correspond to the binned photometric data. The x axis represents the time computed from the mid-transit times as derived from the best joint fit.
In the text

	Fig. A.6 Posterior distributions of some of the fitted planetary parameters of the TOI-1238 system as obtained from the combined photometric and spectroscopic fit. The vertical dashed lines indicate the 16, 50, and 84 % quantiles that were used to define the optimal values and their associated 1σ uncertainty. The blue line stands for the median values of each fitted parameter.
In the text