TOI-1736 and TOI-2141: two systems including sub-Neptunes around solar analogs revealed by TESS and SOPHIE^††thanks: Based on observations collected with the SOPHIE spectrograph on the 1.93 m telescope at the Observatoire de Haute-Provence (CNRS), France.

The TESS mission observed TOI-1736 with a cadence of 2 minutes in Sectors 18, 19, and 25, and with a cadence of 20 seconds in Sectors 52, 58, and 59. TOI-2141 was observed with a cadence of 2 minutes in Sectors 25, 26, and 52. Table 2 shows the log of TESS observations for these two objects. We obtained the TESS data products from the Mikulski Archive for Space Telescopes (MAST)²²2mast.stsci.edu, where we used the Presearch Data Conditioning (PDC) flux time series (Smith et al., 2012; Stumpe et al., 2012, 2014) processed by the TESS Science Processing Operations Center (SPOC) pipeline (Jenkins et al., 2016) versions listed in Table 2.

The SPOC searches light curves for transiting planets with an adaptive, noise-compensating matched filter (Jenkins, 2002; Jenkins et al., 2010, 2020). The transit signatures of each candidate are fitted with an initial limb-darkened transit model (Li et al., 2019) and subjected to a suite of diagnostic tests (Twicken et al., 2018), all available in the TESS SPOC data validation reports (DVR). The DVR for TOI-1736 reports a candidate planet, hereafter TOI-1736 b, with an estimated radius of R${}_{\rm p}=3.0\pm 1.1$ R_♁ and an orbital period of P$=7.07307(1)$ d. The difference image centroid offsets locate the source of the transit signal within $4.4\pm 4.1$ arcsec of the target star. The TESS Science Office (TSO) issued an alert for TOI 1736.01 on 27 February 2020 based on the DVR associated with the combined light curve for Sectors 18-19 (Guerrero et al., 2021).

The DVR for TOI-2141 reports two candidates, but only one passed validation tests. The unvalidated events have a periodicity of 30.6 d and are likely due to contamination by scattered light. We also note that the odd and even depth test supports a single strong feature folded on top of a much weaker and much less convincing transit-like feature for these events. The validated candidate planet, hereafter TOI-2141 b, has an estimated radius of R${}_{\rm p}=3.2\pm 0.3$ R_♁ and an orbital period of P$=18.26159(7)$ d. The centroid offsets localize the source of the transit signal within $1.7\pm 6.3$ arcsec. The transit signature was first identified in the SPOC search of the Sector 26 light curve, and an alert for TOI 2141.01 was issued by TSO on 7 August 2020. Therefore, both TOI-1736 and TOI-2141 systems host at least one close-in transiting planet candidate with a size consistent with a sub-Neptune. Figures 1 and 2 depict the TESS photometry data for both respective objects.

The TESS pixel scale is $\sim 21\arcsec$ pixel^-1, and photometric apertures typically extend out to roughly 1 arcminute, which generally results in multiple stars blending in the TESS aperture. To attempt to determine the true source of the TESS detection, we conducted ground-based photometric follow-up observations of the field around TOI-1736 as part of the TESS Follow-up Observing Program³³3https://tess.mit.edu/followup Sub Group 1 (TFOP; Collins, 2019).

We observed a full predicted transit window of TOI-1736 b in Pan-STARRS $z$-short band using the Las Cumbres Observatory Global Telescope (LCOGT; Brown et al., 2013) 1.0 m network node at McDonald Observatory on UTC 2020 August 28. If the event detected in the TESS data is indeed on-target, the shallow SPOC reported depth of 380 ppm would not generally be detectable in ground-based observations. Instead, we slightly saturated TOI-1736 to enable the extraction of light curves of nearby fainter stars to attempt to rule out or identify nearby eclipsing binaries (NEBs) as potential sources of the TESS detection. The 1 m telescopes are equipped with $4096\times 4096$ SINISTRO cameras having an image scale of $0\aas@@fstack{\prime\prime}389$ per pixel, resulting in a $26\arcmin\times 26\arcmin$ field of view. The images were calibrated by the standard LCOGT BANZAI pipeline (McCully et al., 2018), and photometric data were extracted using AstroImageJ (Collins et al., 2017).

To account for possible contamination from the wings of neighboring star point spread functions (PSFs), we searched for NEBs in all known Gaia DR3 and TICv8 nearby stars out to $2\aas@@fstack{\prime}5$ from TOI-1736 that are possibly bright enough in the TESS band to produce the TESS detection (assuming a 100% eclipse and 100% contamination of the TESS aperture). To attempt to account for possible delta-magnitude differences between the TESS band and the follow-up filter band, we checked stars that are an extra 0.5 magnitudes fainter in TESS-band than needed. We find that the RMS of each of the light curves of the 24 stars matching our criteria is more than a factor of 5 smaller than the expected NEB depth in the respective star. We then visually inspected each neighboring star’s light curve to ensure no obvious eclipse-like signal. All of our follow-up light curves and supporting results are available on the EXOFOP-TESS website⁴⁴4https://exofop.ipac.caltech.edu/tess/target.php?id=408618999. Through our process of elimination, we find that the TESS signal must be occurring in TOI-1736 relative to known Gaia DR3 and TICv8 stars.

High-angular-resolution observations can probe close companions within $\sim 1.2$ arcsec that can create a false positive transit signal if that companion is an eclipsing binary, and which dilute the transit signal and thus yield underestimated planet radii (Ciardi et al., 2015). TOI-1736 was observed on October 19, 2021, by the ‘Alopeke dual-channel speckle imaging instrument on Gemini-N (PI: Howell) with a pixel scale of 0.01 arcsec/pixel and a full width at half maximum (FWHM) resolution of 0.02 arcsec. ‘Alopeke provided simultaneous speckle imaging at 562 and 832 nm. The data were processed with the speckle pipeline (Howell et al., 2011), which yielded the 5-sigma sensitivity curves shown in Figure 3. These observations provide a contrast at an angular separation of 0.5 arcsec of 4.2 mag at 562 nm and 6.4 mag at 832 nm. Figure 3 also shows the reconstructed image at 832 nm and its Fourier transform, which shows fringes in the power spectrum, evidence of a stellar companion to TOI-1736. The best-fit binary model for the power spectrum gives a companion at an angular separation of $0.093\pm 0.002$ arcsec and position angle of $26.5\pm 1.0$ deg with a magnitude difference from the primary of $\Delta{\rm mag}=2.43\pm 0.15$ at 832 nm.

TOI-2141 was observed on April 25, 2021, with the 4.1-m SOAR telescope in speckle imaging. These observations provide a contrast of 5 mag at an angular separation of 1.0 arcsec in the I-band (Figure 4). These observations did not reveal any evidence of close stellar companions to TOI-2141.

For spectroscopic reconnaissance, we obtained four spectra of TOI-1736 and two spectra of TOI-2141 using the 1.5 m Tillinghast Reflector Echelle Spectrograph (TRES; Fűrész et al., 2008) located at the Fred Lawrence Whipple Observatory (FLWO) in Arizona. The spectra for TOI-1736 and TOI-2141 were obtained between March 3 - October 4, 2020, and between August 17 - October 20, 2020. TRES is a fiber-fed echelle spectrograph with a wavelength range of 390-910 nm and a resolving power of R=44,000. Spectra were extracted and reduced as described in Buchhave et al. (2010).

We performed an analysis of the broadband spectral energy distribution (SED) of both stars including the Gaia Early Data Release 3 parallaxes (with no systematic offset applied; see, e.g., Stassun & Torres, 2021) in order to determine an empirical measurement of the stellar radius, following the procedures described in Stassun & Torres (2016), Stassun et al. (2017), and Stassun et al. (2018).

We performed a fit using NextGen stellar atmosphere models (Hauschildt et al., 1999), with the effective temperature ($T_{\rm eff}$), metallicity ([Fe/H]), and extinction ($A_{V}$) as free parameters. The surface gravity ($\log g$) has little influence on the broadband SED.

As a test to investigate whether the companion to TOI-1736 detected by the high contrast imaging results presented in Section 2.3 is real, we first fit a single-component model to the SED, as shown in the top panel of Figure 5, where we obtained a reduced $\chi^{2}$ of 1.1, with the best-fit parameters $T_{\rm eff}=5800\pm 75$ K, [Fe/H]=$0.2\pm 0.1$, and $A_{V}=0.15\pm 0.03$. According to the tool Stilism⁶⁶6https://stilism.obspm.fr/ (Lallement et al., 2014), the expected extinction for the distance of $\sim 90$ pc, galactic latitude and longitude from TOI-1736 is $A_{V}=0.003$ mag, therefore an abnormal reddening (or an infrared excess flux) was detected in our model. Then we consider a two-component stellar model to fit the SED of TOI-1736, assuming that the secondary has the same extinction ($A_{V}=0.003$) and metallicity as the primary. The bottom panel of Figure 5 shows these results, which give a $\chi^{2}$ of 1.0 with the same fit parameters for the primary and $T_{\rm eff}=4800\pm 200$ K for the secondary. Hence, the SED of TOI-1736 provides additional support to the binary hypothesis. For TOI-2141, we adopt a single-component model as illustrated in Figure 6, which gives a reduced $\chi^{2}$ of 1.3 and the best fit parameters $T_{\rm eff}=5750\pm 75$ K, [Fe/H]=$-0.1\pm 0.1$, and $A_{V}=0.06\pm 0.04$.

Integrating the SED model gives the bolometric flux on Earth ($F_{\rm bol}$), which can be combined with the $T_{\rm eff}$ and the Gaia parallax to derive the stellar radius. For TOI-1736, the single component model gives a radius of $R_{\star}=1.390\pm 0.037$ R_⊙ and the two-component model gives a $3\sigma$ smaller radius of $R_{\star}=1.241\pm 0.032$ R_⊙ for the primary and $R_{\star}=0.748\pm 0.073$ R_⊙ for the secondary. For TOI-2141, we obtained a radius of $R_{\star}=0.938\pm 0.027$ R_⊙. Using the mass-radius relations of Torres et al. (2010), we calculate the stellar masses of $1.15\pm 0.07$ M_⊙ for TOI-1736 and $0.98\pm 0.06$ M_⊙ for TOI-2141.

To obtain the stellar parameters from the SOPHIE data, we constructed an average spectrum for both stars, using all SOPHIE spectra that were unaffected by the Moon’s pollution, and performed a standard spectral analysis on them. The average spectra of TOI-1736 and TOI-2141 have an S/N at 649 nm of 1100 and 850, respectively. We use the methods presented in Santos et al. (2004) and Sousa et al. (2008) to derive the effective temperatures $T_{\rm eff}$, the surface gravities $\log g$, and the metallicities ${\rm[Fe/H]}$. Using these spectroscopic parameters as input, the stellar masses were derived from the Torres et al. (2010) calibration with a correction following Santos et al. (2013). The errors were calculated from 10,000 random realizations of the stellar parameters within their error bars and assuming Gaussian distributions. The parameters obtained from this analysis are presented in Table 4 labeled as “standard.”

We use the reduced TRES spectra to obtain the stellar parameters using the Stellar Parameter Classification tool (SPC; Buchhave et al., 2012). In short, SPC correlates each observed spectrum against a grid of synthetic spectra based on Kurucz (Kurucz, 1992) atmospheric models and derives the effective temperature, surface gravity, metallicity, and rotational velocity of the star. The parameters obtained from this analysis are presented in Table 4.

We take advantage of the spectroscopic similarity between our targets and the Sun to apply a strictly differential analysis to determine atmospheric parameters and elemental abundances (e.g., Ramírez et al., 2014; Bedell et al., 2014). The differential analysis in such cases can potentially mitigate the Fe I/II excitation-ionization balance internal errors caused by possible differences in the micro-physics prescriptions, in the local determination of the continuum, or in the equivalent width (EW) measurements. Our reference solar spectrum was obtained from observations of the Moon with SOPHIE on 2023-04-01 with an S/N of 245 at 649 nm, maintaining the same instrumental configuration as the one used to observe the scientific targets.

We perform a line-by-line analysis with respect to the solar spectrum, which is based on the following steps. Given a narrow spectral range around each selected line, we normalize the stellar and solar spectrum using the same local continuum regions and then fit a Gaussian function to the solar line profiles, computing all resulting parameters (e.g., line depth, centroid, FWHM, and continuum offset). Then, we repeat the line fitting procedure for the locally normalized stellar spectrum, now fixing the continuum offset and line region information obtained from the respective solar spectrum. We employ the ionization-excitation equilibrium method where we adopted a line list from Meléndez et al. (2014), which contains excitation potentials, oscillator strengths, laboratory $\log{gf}$ values, and hyperfine structure corrections⁷⁷7albeit it has a small impact on stellar abundances of solar analogs, in the Y II lines we corrected for this effect.. We included iron-peak, s- and alpha-capture elements, such as Fe, Mg, Al, Si, Ti, Y, and Ca to perform this spectroscopic analysis.

In Figure 7 we show the relationship between the solar Fe I and Fe II EWs and those of TOI-1736 and TOI-2141. The low scatter ($\sigma{\rm EW}\sim 1-2\%$) for both stars indicates that the differential spectroscopic equilibrium analysis with respect to the Sun is sufficiently precise for the determination of the stellar atmospheric parameters. While not an exact match to our Sun, exercising caution and maintaining a conservative approach in considering uncertainties is warranted.

We use the code q2 (Ramírez et al., 2014) to model both the stellar and solar EWs of Fe I and Fe II from our line-by-line analysis. The input-relevant line data (e.g., oscillator strength, excitation potential, etc.) are also provided. The q2 was configured to employ the Kurucz ODFNEW model atmospheres (Castelli & Kurucz, 2003) and the 2019 version of the local thermodynamic equilibrium (LTE) code MOOG (Sneden, 1973; Sneden et al., 2012). The stellar parameters are obtained through an iterative convergence process, where q2 iterates until it finds: 1) no systematic dependence of the derived Fe I and Fe II abundances as a function of excitation potential and reduced EW, constraining the $T_{\rm eff}$ and micro-turbulence; 2) Fe I and Fe II yielding the same abundance value, constraining the $\log{g}$. Once convergence is achieved, the final atmospheric parameters are estimated together with their uncertainties, which include systematic errors, as detailed in Appendix B. For TOI-1736 we obtained $T_{\rm eff}=5807\pm 13$ K, $\log{g}=4.35\pm 0.04$ dex, $[{\rm Fe}/{\rm H}]=+0.138\pm 0.014$ dex, and $v_{\rm t}=1.26\pm 0.02$ km s^-1. For TOI-2141 we obtained $T_{\rm eff}=5659\pm 12$ K, $\log{g}=4.42\pm 0.04$ dex, $[{\rm Fe}/{\rm H}]=-0.120\pm 0.013$ dex, and $v_{\rm t}=0.94\pm 0.02$ km s^-1.

Stellar elemental abundances are obtained by q2 using local thermodynamic equilibrium (LTE) model interpolations based on the previously derived atmospheric parameters. The adopted errors for each element account for the errors from atmospheric parameter measurements and the scatter from the line-to-line abundance estimates. Table 3 summarizes the abundances obtained in our analysis.

We performed an evolutionary analysis using the derived spectroscopic parameters (T_eff, $\log{g}$ and abundances), where we adopted the Yonsei-Yale evolutionary tracks (Yi et al., 2001; Demarque et al., 2004) that account for the contribution of the alpha enhancement elements to compute the ages, masses, and radii of TOI-1736 and TOI-2141. The derived atmospheric parameters are compared with those predicted by stellar models using the Bayesian framework described in Grieves et al. (2018) and Yana Galarza et al. (2021). Table 4 shows the derived values of stellar mass and radius and the derived ages are shown in Table 5.

As we have employed different methods and used data from different instruments to obtain the stellar parameters, in this section we make a comparison between these results. At the time of writing, MacDougall et al. (2023) has also published independent results for the parameters of TOI-1736 ($T_{\rm eff}$ and $[{\rm Fe}/{\rm H}]$), which we also include here for comparison.

For TOI-1736, the mean values of $T_{\rm eff}$ range from 5636 to 5807 K, where the SED, standard, and solar diff. agree within $1\,\sigma$, whereas TRES and MacDougall et al. (2023) agree within $3\,\sigma$ with respect to the other values. All surface gravity values also agree within $3\,\sigma$, although the solar diff. gives a higher value than other methods. The GAIA eDR3 trigonometric value is slightly more precise than the solar diff, but it can be contaminated by the close companion in the same way as the SED value. All metallicity values agree within $1\,\sigma$. For TOI-2141, all but the SED value of $T_{\rm eff}$ agree within $1\,\sigma$, and all agree within $2\,\sigma$. All surface gravity and metallicity values also agree within $1\,\sigma$. The parameters of this star are more similiar to our Sun, being poorer in metals and more active than normal for its age.

Different instruments/methods may have different offsets for derived parameters. For this reason, we performed an absolute calibration of our differential analysis of the SOPHIE data using other stars and the Sun as a reference (see Appendix B). Since the determination is more precise and since the other methods generally agree with it, we adopt the parameters of the solar differential analysis as final values. Star masses and radii are derived from spectroscopic parameters, so they carry more or less the same level of discordance as those discussed above. For consistency, we adopted the final mass and radius values as those obtained from the solar differential analysis. As pointed out by (Tayar et al., 2022), we recognize that even the differential analysis has systematic errors in the models of stellar evolution that were not taken into account in our analysis. Therefore, we use Tayar et al. (2022)’s methods ⁹⁹9https://github.com/zclaytor/kiauhoku/blob/v1.4.0/notebooks/model_offsets.ipynb to calculate the additional errors in stellar parameters due to systematic differences in the models, as shown in Table 4. We add these errors in quadrature to the uncertainties obtained by our differential analysis, giving a more realistic uncertainty in the final values.

Turbulence velocity and rotational velocity are two quantities that depend on the adopted model and can absorb some instrumental broadening. Thus, to avoid any bias toward a specific instrument/model, we adopted their mean values.

The CCF analysis performed on the SOPHIE spectra provides three quantities that are sensitive to stellar activity: the FWHM, bisector span (BIS), and RVs (e.g., Queloz et al., 2001; Boisse et al., 2009), as shown in Appendix A. In addition, we measured two well-known spectral indices that are proxies for chromospheric and coronal activity, the S-index, which relies on the emissions in the cores of the H and K Ca II lines, and H$\alpha$ (e.g., Cincunegui et al., 2007). In Appendix C we present the details of our measurements of the S-index and H$\alpha$ from the SOPHIE spectra. In particular, we present a calibration of the S-index to the Mount-Wilson system. The temporal variability of these activity indicators can be used to identify systematic errors in the RVs caused by spurious signals of stellar activity, thus allowing to improve the determination of planetary orbits, as will be explored in Section 4.

We convert the Mount-Wilson S-index to the chromospheric activity index $\log{\rm R^{\prime}}_{\rm HK}$ using the Czesla et al. (2019) recipe that includes the photospheric correction (e.g., Mittag et al., 2013), where we use the $B-V$ values from the magnitudes listed in Table 1 and T_eff from our analysis. The result was $\log{\rm R^{\prime}}_{\rm HK}(B-V)=-4.94^{+0.15}_{-0.25}$ for TOI-1736 and $\log{\rm R^{\prime}}_{\rm HK}(B-V)=-4.78^{+0.13}_{-0.17}$ for TOI-2141. Using the activity-rotation empirical calibration of Mamajek & Hillenbrand (2008) we obtained $P_{\rm rot}=28\pm 5$ d for TOI-1736 and $P_{\rm rot}=21\pm 5$ d for TOI-2141.

Also using the calibration of Mamajek & Hillenbrand (2008) we obtained gyrochronological ages of $4.6^{+1.7}_{-1.3}$ Gyr for TOI-1736 and $2.7^{+1.6}_{-1.0}$ Gyr for TOI-2141. Applying the age-chromospheric activity relation for solar analogs of Lorenzo-Oliveira et al. (2016) we estimated the ages in $5.4^{+1.4}_{-1.1}$ Gyr for TOI-1736 and $4.2^{+1.6}_{-1.1}$ Gyr for TOI-2141, which are consistent with those from gyrochronology. We adopted the activity ages as the latter, since Lorenzo-Oliveira et al. (2018) take mass and metallicity into account in their empirical relationship, proving to be more accurate for the range of parameters of our targets.

Table 5 shows the derived age values of the evolutionary analysis of stellar parameters. As an additional check, we calculate the stellar ages from chemical clocks (e.g., da Silva et al., 2012; Tucci Maia et al., 2015; Delgado Mena et al., 2019). We calculated the yttrium abundance ratios ([Y/X], where X=Si I, Ti I, Ti II, or Al I) from the values reported in Table 3, where we estimated the stellar ages using the Delgado Mena et al. (2019) relations that account for T_eff and [Fe/H] effects on chemical ages. The average calibration errors are about 1.4 Gyr and the impact of T_eff-[Fe/H]-[YII/X] errors are respectively 0.5 and 0.3 Gyr for TOI-1736 and TOI-2141. The average chemical ages are shown in Table 5.

Lithium (Li) is also known to be a sensitive indicator of stellar evolution (e.g., Carlos et al., 2016; Lyubimkov, 2016). To obtain Li abundance, we performed a spectral synthesis analysis using the code iSpec¹⁰¹⁰10https://www.blancocuaresma.com/s/iSpec (Blanco-Cuaresma et al., 2014; Blanco-Cuaresma, 2019) and the radiative transfer code MOOG(Sneden, 1973; Sneden et al., 2012), where we use the lines in the spectral range around Li from Meléndez et al. (2012). Figure 8 illustrates our results, where we show the SOPHIE spectra and the best-fit synthetic models for TOI-1736 and TOI-2141, and also for a solar spectrum obtained from observations of the Moon. We calculated a Li abundance for the solar spectrum of ${\rm A(Li)}=1.1\pm 0.2$, which is consistent with the expected value for the Sun. For TOI-2141, our results are compatible with a null value, giving an upper limit of A(Li)$<0.5$. For TOI-1736, we obtained an abundance of ${\rm A(Li)}=1.56\pm 0.16$. Using the relationship between A(Li) and age from Carlos et al. (2019) we obtained an age of $4.4\pm 1.0$ Gyr for TOI-1736.

All the different age determinations for both systems are presented in Table 5. There is general agreement between the ages of TOI-1736, which are in the range of 4-6 Gyr. For TOI-2141, both the chemical and evolutionary ages agree around 7.5 Gyr, which is consistent with the absence of Li, but the activity-derived age appears younger. For a final determination of the ages of both stars, we adopted the average values, where we added Tayar et al. (2022)’s model uncertainties in quadrature to the measured errors, as in Section 3.5, giving an age of $4.9\pm 1.3$ Gyr for TOI-1736 and $6.4\pm 1.8$ Gyr for TOI-2141.

To further characterize the transits detected by TESS, we analyzed the PDCSAP flux data by applying the methods described in Martioli et al. (2021), where we fit a transit model plus a baseline polynomial in selected windows around each transit in such a way that the window size is twice the duration of the transit and with a minimum of 300 data points. We selected windows for a total of 19 transits of TOI-1736 b and four transits of TOI-2141 b. Our transit model is calculated using the BATMAN toolkit (Kreidberg, 2015) and the posterior distribution of transit parameters is sampled using a Bayesian Monte Carlo Markov Chain (MCMC) framework with the package emcee (Foreman-Mackey et al., 2013). We use uninformative priors for the transit parameters, as presented in Table 12 of Appendix D, with initial values given by those reported in the TESS DVRs and assuming circular orbits.

We divided the TESS light curve by the best-fit transit model and binned the data by the weighted average with bin sizes of 0.1 d to fit a baseline flux using Gaussian processes (GP) regression with a quasi-periodic (QP) kernel as in Martioli et al. (2022). This kernel was chosen to detect a possible variability modulated by the star’s rotation period. The priors of the GP parameters are listed in Table 12 and the posteriors are listed in Table 13. We opted to fix the smoothing factor and decay time values at 0.1 and 10 days, respectively, as altering these values does not impact the other parameters. Unfortunately, the periodicities obtained in our analysis are not well constrained, and their values change significantly depending on the initial guess adopted for our GP model. The low level of intrinsic photometric variability in both targets, coupled with the rotation periods expected from stellar activity (16 – 32 days) being closely aligned with the TESS sector duration, strongly biases the detection of periodicity. This bias is primarily due to systematic variations that occur during TESS observations, particularly near the edges of each sector, where the amplitude of these variations is greater. While there may be some intrinsic stellar variability captured by TESS photometry, the GP model obtained here may account for both systematic and intrinsic stellar phenomena in unknown proportions.

Since the periodic kernel is not essential, the primary purpose of employing the GP in this analysis is for detrending. To ensure that the choice of GP does not impact the planet parameters, we conducted experiments using both a squared exponential (SE) kernel and a model without GP. The posterior distributions of key planet parameters ($T_{c}$, $a/R_{\star}$, $P$, and $R_{p}/R_{\star}$) are consistent within 1-sigma for all three choices of baseline model. The results of our best-fit transit models and GP analysis using a QP kernel are illustrated in Figures 1 and 2.

The TOI-1736 RVs have a median value of -25408.95 m s^-1 with an rms of 158 m s^-1 and a median individual error of 2 m s^-1. These RVs present a clear signal with an amplitude of hundreds of m s^-1, as can be seen in Figure 13. An initial analysis of the generalized Lomb-Scargle (GLS) periodogram (Zechmeister & Kürster, 2009) of this data show that there is a dominating signal with periodicity around 570 d and a long-term linear trend was also detected by visual inspection. To detect the TESS candidate planet, TOI-1736.01, at a period of 7.073 d, we fit the SOPHIE RV data using the online tool Data & Analysis Center for Exoplanets (DACE¹²¹²12https://dace.unige.ch), where we adopted a two-planet Keplerian model plus a linear trend. As illustrated in Figure 9, the GLS of the RVs after subtracting the linear trend and each planet’s orbit shows a strong peak at 7.0731 d and 570.07 d for the inner and outer planet, respectively. The former agrees well with the periodicity of the transits detected by TESS, showing the RV detection of the inner planet candidate, TOI-1736 b, and an outer giant planet, TOI-1736 c. Moreover, this confirms that transits detected by TESS occur in the primary star, as the dominant flux in the SOPHIE spectra comes from the primary and not from the K3 dwarf companion. As a way to test if the SOPHIE spectra are contaminated by the flux from the companion, we performed a CCF analysis using a K0, K5, or M0 masks, which should favor the detection of spectral lines of the cooler star. Those tests resulted in slightly noisier CCF with parameters in agreement as when using the G2 mask, showing that indeed, the flux contribution of the companion to the RV measurements in the SOPHIE spectra is negligible.

The TOI-2141 RVs have a median value of -19860.8 m s^-1 with an rms of 7 m s^-1 and a median individual error of 2.4 m s^-1. The RVs of this object do not appear to have any long-term trends or variability. Figure 10 shows the GLS periodogram for this data, which shows the highest peak at a period of 18.259 d, again in agreement with the periodicity of the transits detected by TESS. This shows the detection of the planet candidate TOI-2141 b in our SOPHIE data. In the following sections, we present an analysis to characterize these two planetary systems.

For both systems, we found that there is no correlation between the RVs and the four stellar activity indicators: FWHM, BIS, H_α, and S-index. Therefore, the signals detected in RV are probably due to Doppler shift and not to scenarios that imply changes in the line profile, for example, blended binaries. We have also inspected possible correlations between the RV residuals (after subtracting the fit models) and activity indicators. Only BIS shows a slight correlation of 0.29 and 0.15 for TOI-1736 and TOI-2141, respectively. We tried to de-correlate the RVs as in Boisse et al. (2009), but this did not have a significant impact on the values of the fit parameters, so we preferred to keep the final analysis without applying this correction. In the Appendix C we show the time series of each activity indicator, an analysis of the GLS periodogram, and its correlations with the RV data.

To obtain the physical parameters of the planets, we perform a joint analysis of the photometry and RV data sets as in Martioli et al. (2022). We adopt the priors listed in Table 12 and the initial values are those obtained by the analysis performed on the SOPHIE RVs and TESS photometry data independently.

We first fit the RV orbits and transit models jointly using the scipy.optimize.leastsq optimization tool. Then, we use the same tool to fit the RV jitter, which is a term quadratically added to the RV errors. Finally, we explore the full range covered by the priors with a Bayesian MCMC that uses 50 random walkers, with 10000 iterations, discarding the first 3000 burn-in samples. The posterior distributions of sampled parameters are illustrated in Appendix D and the resulting best-fit parameters, obtained from the mode of the posterior and 1$\sigma$ uncertainties (34% on each side of the central value) are given in Table 6. Some parameters converged to values close to the edge of their physical domains. For example, it happened for the orbital inclination ($i_{b}\sim 90$ deg) and the linear limb-darkening coefficient ($u_{0}\sim 0$) for both transiting planets, TOI-1736 b and TOI-2141 b. We also adopted circular orbits ($e=0$, $\omega=90$ deg) for both inner planets, as noncircular solutions resulted in higher values of the Bayesian Information Criterion (BIC), suggesting that our data cannot constrain the low values of the eccentricities of these planets.

Figures 12 and 11 illustrate the TESS photometry data for all transits detected by TESS and the best-fit transit models obtained by our joint analysis. Figure 13 and 14 illustrate the TOI-1736 RV data and the best-fit models that include the orbit of the two planets, TOI-1736 b and TOI-1736 c, and the linear trend. Figures 15 and 16 illustrate the RV data and best-fit model for TOI-2141.

We consider the best-fit parameters from our analysis to calculate the habitable zone for these solar analogs, using the equations and data from Kopparapu et al. (2014). We obtained an optimistic lower limit (recent Venus) at 0.87 au for TOI-1736 and 0.71 au for TOI-2141, and an upper limit (early Mars) at 2.05 au for TOI-1736 and 1.68 au for TOI-2141. The runaway greenhouse limits ($M_{\rm p}=1$ M_♁) ranges between 1.10 au and 1.94 au for TOI-1736, and between 0.90 au and 1.59 au for TOI-2141.

To estimate the level of radiation that each planet receives from its host star, we calculate the effective insolation flux $S_{\rm eff}$ relative to the flux incident on Earth $S_{\oplus}$ as $S_{\rm eff}/S_{\oplus}=L_{\star}/a_{p}^{2}$, where $L_{\star}$ is the star luminosity in solar units and $a_{p}$ is the semi-major axis of the planet’s orbit in au. Figure 19 shows the insolation for exoplanets as a function of star effective temperature. TOI-1736 b and TOI-2141 b receive an insolation of $246\pm 26\,S_{\oplus}$ and $50\pm 4\,S_{\oplus}$, respectively, showing that these planets reside in a hot environment. We estimate the equilibrium temperature for the planets as in Heng & Demory (2013), assuming a uniform heat redistribution and an arbitrary geometric albedo of 0.1, which gives $T_{\rm eq}=1076\pm 39$ K and $T_{\rm eq}=722\pm 23$ K for TOI-1736 b and TOI-2141 b, respectively.

On the other hand, TOI-1736 c resides at a larger distance from the star, in an eccentric orbit with $e_{c}=0.362\pm 0.003$, and a semi-minor axis of $1.287\pm 0.016$ au and a semi-major axis of $1.381\pm 0.017$ au. At these distances, TOI-1736 c receives an insolation of $0.81\pm 0.09\,S_{\oplus}$ and $0.71\pm 0.08\,S_{\oplus}$ for which we estimate an equilibrium temperature of $T_{\rm eq}=258\pm 9$ K and $T_{\rm eq}=249\pm 9$ K, respectively. This planet resides within the conservative habitable zone around the star.

From mass and radius measurements of the transiting planets TOI-1736 b and TOI 2141 b, we derive their mean densities of $4.9\pm 1.3$ g cm^-3 and $4.6\pm 1.3$ g cm^-3, respectively, showing that these sub-Neptune sized planets have densities similar to the rocky planets of the Solar System. To explore further the possible chemical composition and structure of these planets we inspect the mass-radius (M-R) diagram illustrated in Figure 20, which shows the exoplanet data from exoplanets.eu, the data from the two hot sub-Neptunes studied here, and a set of models from Zeng et al. (2019), for comparison. The M-R location of TOI-1736 b matches two models: (1) a 50% Earth-like rocky core plus a 50% H₂O layer; (2) an almost pure Earth-like rocky core with a thin 0.1% H₂ envelop. Therefore, this planet probably has a large dense rocky core surrounded by either a thick water-rich atmosphere or a thin H₂ atmosphere. TOI-2141 b is compatible with three scenarios: (1) a 50% Earth-like rocky core with a 50% H₂O layer; (2) a 49.95% Earth-like rocky core + 49.95% H₂O layer and a thin 0.1% H₂ envelop; and (3) a pure H₂O. Thus, TOI-2141 b is likely a water-rich planet. A further investigation of the internal structure of these planets can be done by including simulations constrained by the star metallicity and abundances of volatiles, which is out of the scope of this paper.