Mass determinations of the three mini-Neptunes transiting TOI-125

ABSTRACT

The Transiting Exoplanet Survey Satellite, TESS, is currently carrying out an all-sky search for small planets transiting bright stars. In the first year of the TESS survey, a steady progress was made in achieving the mission’s primary science goal of establishing bulk densities for 50 planets smaller than Neptune. During that year, the TESS’s observations were focused on the southern ecliptic hemisphere, resulting in the discovery of three mini-Neptunes orbiting the star TOI-125, a V = 11.0 K0 dwarf. We present intensive HARPS radial velocity observations, yielding precise mass measurements for TOI-125b, TOI-125c, and TOI-125d. TOI-125b has an orbital period of 4.65 d, a radius of 2.726 ± 0.075 R_E, a mass of 9.50 ± 0.88 M_E, and is near the 2:1 mean motion resonance with TOI-125c at 9.15 d. TOI-125c has a similar radius of 2.759 ± 0.10 R_E and a mass of 6.63 ± 0.99 M_E, being the puffiest of the three planets. TOI-125d has an orbital period of 19.98 d and a radius of 2.93 ± 0.17 R_E and mass 13.6 ± 1.2 M_E. For TOI-125b and d, we find unusual high eccentricities of 0.19 ± 0.04 and |$0.17^{+0.08}_{-0.06}$|⁠, respectively. Our analysis also provides upper mass limits for the two low-SNR planet candidates in the system; for TOI-125.04 (R_P = 1.36 R_E, P = 0.53 d), we find a 2σ upper mass limit of 1.6 M_E, whereas TOI-125.05 (⁠|$R_P=4.2^{+2.4}_{-1.4}$| R_E, P = 13.28 d) is unlikely a viable planet candidate with an upper mass limit of 2.7 M_E. We discuss the internal structure of the three confirmed planets, as well as dynamical stability and system architecture for this intriguing exoplanet system.

Planets and satellites: detection, Planets and satellites: individual: (TOI-125, TIC 52368076)

1 INTRODUCTION

The Transiting Exoplanet Survey Satellite (TESS – Ricker et al. 2015) is more than halfway through a survey of about 85 per cent of the sky. More than 1000 planet candidates have been announced so far. The Level-1 mission goal of TESS is to measure the masses and radii of at least 50 exoplanets smaller than 4 R_E. Among the first planets that meet the Level-1 requirement are HD 15337b and c (TOI-402, Dumusque et al. 2019; Gandolfi et al. 2019), HD 21749b (TOI-186, GJ 143, Dragomir et al. 2019; Trifonov, Rybizki & Kürster 2019), GJ 357 b (TOI-562, Luque et al. 2019), LTT 1445Ab (Winters et al. 2019), HD 23472 b and c (TOI-174, Trifonov et al. 2019), and π Men c (HD 39091, Gandolfi et al. 2018; Huang et al. 2018).

TESS is building on top of a great legacy from Kepler (Borucki et al. 2010), which detected numerous multiplanet systems for which system architecture has been studied in detail, e.g. Lissauer et al. (2011). The identification of the distinct populations of super-Earths and mini-Neptunes separated by a valley caused by stellar irradiation evaporating the planet atmosphere (Fulton et al. 2017; Owen & Wu 2017; Fulton & Petigura 2018) is also owed to Kepler. This process can potentially strip a planet down to its core. Multiplanet systems provide prime target for testing both bulk composition models and atmospheric evaporation, and are thus crucial for advancing exoplanet science.

We present the confirmation and precise mass measurements of three mini-Neptunes orbiting the bright (V = 11.0 mag) K0 dwarf star TOI-125; see Table 1 for a full summary of the stellar properties. This work builds largely on intensive radial velocity (RV) follow-up observations with HARPS (Mayor et al. 2003). The three planets all fall within the TESS Level-1 mission goal, with similar radii but quite different masses. The system was previously validated by Quinn et al. (2019), so the main focus of this paper is the mass characterization presented in Section 3, analysis of the system architecture presented in Section 4, and internal structure presented in Section 5. Finally, we explore future possibilities for atmospheric characterization in Section 6.

Table 1.

Stellar properties for TOI-125.

Property	Value	Source
Other names
2MASS ID	J01342273-6640328	2MASS
Gaia ID	4698692744355471616	Gaia DR2
TIC ID	52368076	TESS
TOI	TOI-125	TESS
Astrometric properties
R.A.	01:34:22.43	TESS
Dec	−66:40:34.8	TESS
μR.A. (mas yr−1)	−119.800 ± 0.066	Gaia DR2
μDec. (mas yr−1)	−122.953 ± 0.080	Gaia DR2
Parallax (mas)	8.9755 ± 0.0356	Gaia DR2
Distance (pc)	111.40 ± 0.44	Gaia DR2
Photometric properties
V (mag)	11.02 ± 0.01	Tycho
B (mag)	11.72 ± 0.12	Tycho
G (mag)	10.718 ± 0.020	Gaia
T (mag)	10.1985 ± 0.006	TESS
J (mag)	9.466 ± 0.021	2MASS
H (mag)	9.112 ± 0.025	2MASS
Ks (mag)	8.995 ± 0.021	2MASS
W1 (mag)	8.945 ± 0.030	WISE
W2 (mag)	9.006 ± 0.030	WISE
W3 (mag)	8.944 ± 0.030	WISE
W4 (mag)	8.613 ± 0.262	WISE
AV	|$0.032 ^{+0.032}_{-0.023}$|	Section 3.3
Bulk properties		This work
Teff (K)	5320 ± 39	Sections 3.1 and 3.3
Spectral type	K0V	Sections 3.1 and 3.3
log g (cm s−2)	4.516 ± 0.024	Section 3.3
ρ (g cm−3)	|$1.99^{+0.13}_{-0.11}$|	Section 3.3
[Fe/H]	−0.02 ± 0.03	Sections 3.1 and 3.3
v sin i (km s−1)	<1.0 ± 0.5	Section 3.1
Age (Gyr)	6.8 ± 4.3	Section 3.3
Radius (R⊙)	0.848 ± 0.011	Section 3.3
Mass (M⊙)	|$0.859^{+0.044}_{-0.038}$|	Section 3.3

Note.Tycho (Høg et al. 2000); 2MASS (Skrutskie et al. 2006); WISE (Wright et al. 2010); Gaia (Gaia Collaboration 2018).

2 OBSERVATIONS

2.1 TESS photometry

TOI-125 (TIC 52368076) was observed by TESS in Sectors 1 and 2 from 2018 July 25 to September 20. It appeared on CCD1 of camera 3 in Sector 1 and CCD2 of camera 3 in Sector 2.

The data are available with 2 min time sampling (cadence) and were processed by the Science Processing Operations Center (SPOC – Jenkins et al. 2016) to produce calibrated pixels, and light curves. Based on the Data Validation report produced by the transit search conducted by the SPOC (Twicken et al. 2018; Li et al. 2019), two TESS objects of interest, TOI-125b and TOI-125c, were announced by the TESS Science Office (TSO) from Sector 1. This was the first multiplanet-candidate system announced by the TSO. With data from Sector 2, a third planet candidate, TOI-125d, was revealed with one transit observed in each sector.

For transit modelling, we used the publicly available Simple Aperture Photometry flux, after the removal of artefacts and common trends with the Pre-search Data Conditioning (PDC-SAP) algorithm (Twicken et al. 2010; Smith et al. 2012; Stumpe et al. 2014) provided by SPOC. The light-curve precision in both the sectors is 125 ppm, averaged over half an hour, consistent with the value predicted by Sullivan et al. (2015) for a star with apparent TESS magnitude of 10.2. Fig. 1 shows the full 2 min cadence TESS light curve, with data points binned to 10 min overplotted, along with the phase-folded light curves for TOI-125b, TOI-125c, and TOI-125d.

Figure 1.

The TESS data for TOI-125 spanning Sectors 1 and 2. Top panel: full light curve with the 2 min cadence data in light grey and the same data binned to 10 min in dark grey. The binned data surrounding the transits are highlighted in red, yellow, and green. The light curve from the two sectors consists of four segments that each correspond to one TESS orbit of 13.7 d. After each orbit, the spacecraft interrupts observations to downlink the data to the Earth, causing gaps in the data coverage. Furthermore, there are features in the light curve from the momentum dumps of the satellite, which take place approximately every 2.5 d. None of the detected transits occurred during momentum dumps. Bottom panel: phase-folded TESS light curves for TOI-125b, TOI-125c, and TOI-125d, again with 2 min cadence data in grey and binned to 10 min in the same colours as the top panel.

The TOI-125 system was vetted by Quinn et al. (2019) using ground-based photometry, high-angular-resolution imaging, and reconnaissance spectroscopy. TOI-125b and TOI-125c were statistically validated as planets, while TOI-125d (then called TOI-125.03) remained a high-SNR planet candidate, based on only two observed transits. Two additional low-SNR candidates were identified: TOI-125.04, with a period of 0.53 d, making it an ultra-short-period (USP) planet candidate, and TOI-125.05 at 13.28 d. Quinn et al. (2019) stressed that these two candidates are marginal detections, and did not attempt to validate them statistically.

2.2 High-resolution spectroscopy with HARPS

TOI-125 was observed intensively with the HARPS spectrograph (Mayor et al. 2003) on the ESO 3.6 m telescope at La Silla Observatory, Chile, from 2018 September 21 to 2019 January 8. In total, 122 spectra were obtained under programmes 1102.C-0249 (PI: Armstrong), 0101.C-0829/1102.C-0923 (PI: Gandolfi), 0102.C-0525 (PI: Díaz), 0102.C-0451 (PI: Espinoza), and 60.A-9700 (technical time). HARPS is a stabilized high-resolution spectrograph with a resolving power of |$R\sim {115\, 000}$|⁠, capable of sub-m s⁻¹ RV precision. We used the instrument in high-accuracy mode with a 1 arcsec science fibre on the star and a second fibre on sky to monitor the sky-background during exposure. We used a nominal exposure time of 1800 s, which on occasion was adjusted within a range of 800–2100 s depending on sky condition and observation schedule.

RVs were determined with the standard (offline) HARPS data reduction pipeline using a K0 binary mask for the cross-correlation (Pepe et al. 2002), and a K3 template for flux correction to match the slope of the spectra across Echelle orders. We performed the data reduction uniformly for all the data from the six programmes under which data had been acquired, to mitigate any possible RV offsets induced by different data reduction parameters and catalogue coordinates in the FITS headers. With a typical signal-to-noise ratio (SNR) of 55, we achieved an RV precision of 1.5 m s⁻¹. The RV data have been made publicly available through The Data & Analysis Center for Exoplanets (DACE¹) hosted at the University of Geneva. For each epoch the bisector span (BIS), contrast, and full width at half-maximum (FWHM) of the CCF were calculated, as well as the chromospheric activity indicators Ca II H&K, H α, and Na.

In our RV analysis, we excluded data taken on the nights starting 2018 November 25, 26, and 27. On these dates, the ThAr lamp used for wavelength calibration of HARPS was deteriorating and subsequently exchanged on 2018 November 28.² The changing flux ratio between thorium and argon emission lines of the dying ThAr lamp induced a |$2\, \mathrm{m\, s^{-1}\, d^{-1}}$| drift in the wavelength solution of HARPS over 5 d. The problematic data were confirmed by comparing unpublished data from the HARPS-N solar telescope (Dumusque et al. 2015; Collier Cameron et al. 2019) and Helios on HARPS, which also observes the Sun daily. The Helios RVs show a clear drift away from the RVs from the HARPS-N solar telescope on the dates of November 25–27 2018, before returning to a nominal level after the change of the ThAr lamp. We still include spectra taken on those dates in our spectral analysis described in the following Section 3.1.

We clearly detect RV signals for TOI-125b, TOI-125c, and TOI-125d. The top panel in Fig. 2 shows a Lomb–Scargle periodogram of the raw RVs in which there are clear signals at 4.65 and 19.98 d for TOI-125b and TOI-125d with a false-alarm probability (FAP) <0.1 per cent. A hint can be seen for TOI-125c at 9.15 d, possibly interfering with a P/2 alias from TOI-125d. The residuals of a two-planet fit are shown in Fig. 2 below the raw RVs, and show a significant peak at the period of TOI-125c. Peaks at 1.27 and 0.82 d in the raw RVs are aliases of TOI-125b. No signal is found for TOI-125.04 (P = 0.53 d) or TOI-125.05 (P = 13.28 d).

Figure 2.

Lomb–Scargle periodograms, from the top: raw RVs, residuals of a two-planet fit (including TOI-125b and d), FWHM, |$\log R^{\prime }_{\mathrm{ HK}}$|⁠, and bisector span. 1 and 0.1 per cent FAPs are indicated as horizontal lines. Orbital periods for TOI-125b, TOI-125c, and TOI-125d are marked as red, yellow, and green dashed lines, respectively. The expected rotational period of the star is highlighted in grey.

3 DATA ANALYSIS AND RESULTS

3.1 Spectral classification and stellar chemical abundances

The 122 1-D HARPS spectra were stacked to produce a high-fidelity spectrum with SNR per resolution element ∼500 at 5500 Å for spectral analysis. Retrieving stellar parameters from the observed spectrum can be done using several different methods. In the case of TOI-125, stellar atmospheric parameters (T_eff, [Fe/H], and log g) and relative abundances of refractory material were derived using two different methodologies: (a) as described in Sousa et al. (2008) and Santos et al. (2013) using equivalent widths (EWs) of chosen lines while assuming ionization and excitation equilibrium, and (b) with the Spectroscopy Made Easy (SME) code (Valenti & Piskunov 1996; Valenti & Fischer 2005; Piskunov & Valenti 2017) as applied to a grid of model atmospheres.

For the first method, T_eff, [Fe/H], and log g were calculated using the EW of 237 Fe i and 33 Fe ii lines. A grid of Kurucz model atmospheres (Kurucz 1993) and the radiative transfer code moog (Sneden 1973) were used to model the stellar atmosphere. For the derivation of abundances of refractory elements, we used the approach from Adibekyan et al. (2015). TOI-125 shows typical abundances for a main-sequence star, comparable to the ensemble of HARPS GTO stars.

As a second approach, we used SME version 5.22 applied to a grid of MARCS model atmospheres. These are 1D-LTE plane-parallel and spherically symmetric model atmospheres applicable to solar-like stars (Gustafsson et al. 2008). Synthetic spectra were then calculated based on the model grid and fitted to the observed spectral features, focusing on those that are especially sensitive to different photospheric parameters, including T_eff, [Fe/H], log g, micro- and macro-turbulence, and rotational velocity (vsin i). Here, one is changing one or more input parameters and then iteratively using a χ² minimization procedure to arrive at the actual stellar parameters. We used the calibration equation of Bruntt et al. (2010) and Doyle et al. (2014) to estimate the micro- and macro-turbulent velocities, based on the derived values on T_eff and log g. We also fitted 45 isolated and unblended metal lines to determine the projected stellar rotation velocity (vsin i), which was found to be 1.0 ± 0.5 km s⁻¹.

The derived parameters and abundances for both methods are presented in Table 2. It should be noted that the uncertainties were derived from internal errors only, and thus do not include uncertainties inherent to the models themselves. While the abundances and surface gravity log g agree as a whole, there is a 2σ discrepancy between the effective temperature, T_eff, obtained with the two methods. The [Fe/H] measurements also differ slightly between the two methods, but are consistent to 1σ. We have investigated the impact of this on the final set of system parameters, and found less than |$5{{\ \rm per\ cent}}$| difference in stellar and planetary masses and radii. For the final modelling of the system, we used the average of T_eff and [Fe/H] as Gaussian priors in the MCMC. The errors were inflated to encompass both values at a 1σ level, in order to reflect the model dependence of the stellar atmospheric parameters.

3.2 Stellar rotation and activity

The average value of the Ca ii H&K chromospheric activity indicator for TOI-125 is |$\log R^{\prime }_\mathrm{HK}\, =\, -5.00\, \pm \, 0.08$|⁠, indicating a low activity level that would introduce an RV signal on the scale of 0.4 m s⁻¹ (Suárez Mascareño et al. 2017). According to Suárez Mascareño et al. (2015), the expected rotation period of an early K-type dwarf with |$\log R^{\prime }_\mathrm{HK}\, =\, -5.00\, \pm \, 0.08$| is |$P_\mathrm{rot}\, =\, 32^{+5}_{-4} ~\mathrm{d}$|⁠. This is in good agreement with the classical empirical relation from Noyes et al. (1984), which gives |$P_\mathrm{rot}\, =\, 31\, \pm \, 6~\mathrm{d}$|⁠. Assuming that the star is seen equator-on, the projected rotational velocity vsin i = 1 ± 0.5 km s⁻¹ and stellar radius imply a rotation period of ≲43 d. This could be indicative of the stellar spin and the planetary orbits being aligned.

We searched the RVs and activity indicators for a signal matching the expected P_rot. Fig. 2 shows Lomb–Scargle periodograms derived for the raw RVs, RV residuals to a two-planet fit, and RV residuals to a three-planet fit including an additional term fitting a possible 35 d period. We also include periodograms of, FWHM, |$\log R^{\prime }_{\mathrm{ HK}}$| and BIS. FAP thresholds have been computed analytically for levels of 1 and 0.1 per cent. For both the RVs and FWHM, there is an FAP > 0.1 per cent signal close to 40 d, highlighted in grey in Fig. 2. This is in reasonable agreement with the expected stellar rotation period based on |$\log R^{\prime }_{\mathrm{ HK}}$|⁠. The periodogram for |$\log R^{\prime }_{\mathrm{ HK}}$| has signal at 25.5 d, which could be a P_rot/2 alias with FAP 1 per cent. BIS shows no significant signals, though the main peak at 47 d somewhat matches the ones found in FWHM and |$\log R^{\prime }_{\mathrm{ HK}}$|⁠. As a test, we fit three planets along with a 35 d modulation mimicking a signal induced by stellar rotation. The periodogram of the residuals is presented in Fig. 2. It is evident that residual signals at longer periods are still present, including a long-term (P > 100 d) signal we later model as quadratic drift in the RVs.

We searched both the SAP and PDC-SAP light curves for photometric modulation from stellar rotation, but found no convincing signal. This is not too surprising as the baseline of the TESS observation is short (2 × 27 d) compared to the expected rotational period (≥30 d).

Based on the signal seen in both FWHM and RV measurements at 40 d, we attempted to model our RVs with a Gaussian process (GP) trained on the FWHM using a quasi-periodic kernel. The GP had problems converging and the planets’ parameters were unchanged from a classic RV fit. Given the low SNR of both the FWHM and |$\log R^{\prime }_{\mathrm{ HK}}$| indicators, combined with the small expected effect of stellar activity on the RVs, we proceeded to model our data without a GP. Since the period of TOI-125d is about half the stellar rotation period, this might affect the mass measurement of that planet, but we expect this to be a small offset. We can, however, not exclude that the RV semi-amplitude of TOI-125d is slightly affected by stellar activity. Our RV data span several stellar rotations, which to some degree helps mitigate this as we average over epochs with different activity levels.

3.3 Joint modelling with EXOFASTv2

The planetary and stellar parameters were modelled self-consistently through a joint fit of the HARPS RVs and TESS photometry with EXOFASTv2 (Eastman, Gaudi & Agol 2013; Eastman et al. 2019). EXOFASTv2 can fit any number of transits and RV sources for a given number of planets while exploring the vast parameter space through a differential evolution Markov Chain coupled with a Metropolis–Hastings Monte Carlo sampler.

The local χ² minimum in parameter space is identified with AMOEBA, which is a non-linear minimizer using a downhill simplex method (Nelder & Mead 1965). The starting point of the MCMC is set to be within 1σ of the best-fitting value. Hereafter, the full parameter space is explored with a Monte Carlo sampler in numerous steps. At each step the stellar properties are modelled, and limb-darkening coefficients for this specific star are calculated by interpolating tables from Claret & Bloemen (2011). The analytic expressions from Mandel & Agol (2002) are used for the transit model. The eccentricity is parametrized as |$e^{\frac{1}{4}} \cos (\omega _*)$| and |$e^{\frac{1}{4}} \sin (\omega _*)$| to impose uniform eccentricity priors and mitigate Lucy–Sweeney bias of final measurement (Lucy & Sweeney 1971). EXOFASTv2 rejects any solutions where the planetary orbits cross.

At each step χ² is evaluated and assumed to be proportional to the likelihood, which is true for fixed uncertainties. The Metropolis–Hastings algorithm is invoked and 20 per cent of all steps with lower likelihood are kept in the chain. The MCMC thus samples the full posterior distribution.

The size and direction of the next step in the MCMC are determined by the differential evolution Markov Chain method (Ter Braak 2006), where several chains (twice the number of fitted parameters) are run in parallel. The step is determined by the difference between two random chains. In EXOFASTv2, a self-adjusting step size scale is implemented to ensure optimal sampling across the orders of magnitude difference in scales of uncertainty. This is crucial to effectively sample all parameters (e.g. from the orbital period that can be determined to 10⁻⁴ d for transiting planets to the RV semi-amplitude, which commonly can have 10 per cent uncertainty).

The first part of the chains with χ² above the median χ² is discarded as the ‘burn-in’ phase, so as not to bias the final posterior distributions towards the starting point. A built-in Gelman–Rubin statistic (Gelman & Rubin 1992; Gelman et al. 2003; Ford 2006) is used to check the convergence of the chains. When modelling RVs and transit photometry simultaneously, each planet has seven free parameters and up to four additional RV terms for the systemic velocity, drift of the system, and jitter. For the transit light curve, two limb-darkening coefficients for the TESS band are fitted, along with the baseline flux and variance of the light curve.

Another four parameters are fitted for the star: T_eff, [Fe/H], log M_*, and R_*. We applied Gaussian priors on T_eff and [Fe/H] from the spectral analysis, presented in Section 3.1. The mean stellar density is determined from the transit light curve. The Gaia DR2 parallax was used, along with SED fitting to constrain the stellar radius further. We include the broad-band photometry presented in Table 1 in our analysis, apart from the very wide Gaia G band. We set an upper limit on the V-band extinction from Schlegel, Finkbeiner & Davis (1998) and Schlafly & Finkbeiner (2011) to account for reddening along the line of sight. Combining spectroscopic T_eff and [Fe/H] with broad-band SED fitting allow us to perform detailed modelling of the star with the MESA Isochrones and Stellar Tracks (MIST; Choi et al. 2016; Dotter 2016), which are evaluated at each step in the MCMC.

We ran EXOFASTv2 with 50 000 steps on the HARPS RVs and TESS photometry with a quadratic drift in the RVs, with and without eccentricities for TOI-125b, TOI-125c, and TOI-125d. TOI-125b and TOI-125d have significant eccentricity. Fig. 3 displays the HARPS RVs with the final model and Fig. 7 shows a sample of the posterior distribution for the eccentricity of TOI-125b. For simplicity, we fit eccentricities for all three planets in the system.

Figure 3.

HARPS RVs for TOI-125 with a three-planet model including eccentric orbits and a quadratic drift. The residuals to the best fit are shown right below the RV time series. The bottom panel shows the data phase folded and binned for each planet.

The final median values of the posterior distributions and their 1σ confidence intervals for the stellar and planetary parameters are listed in Table 3. We find that TOI-125b has an orbital period of 4.65 d, a radius of 2.726 ± 0.075 R_E, and a mass of 9.50 ± 0.88 M_E, yielding a mean density of 2.57 g cm⁻³. It has the highest orbital eccentricity of the three planet in the system, |$e_\mathrm{ b} = 0.194 ^{+0.041}_{-0.036}$|⁠. With an orbital period of 9.15 d, TOI-125c is near the 2:1 mean motion resonance (MMR) with its inner companion. It has a radius of 2.759 ± 0.10 R_E and a mass of 6.63 ± 0.99 M_E, implying a mean density of 1.73 g cm⁻³. TOI-125d is thus the least dense of the three. Its orbital eccentricity is consistent with zero, |$e_\mathrm{ c} = 0.066^{+0.070}_{-0.047}$|⁠. The outer transiting planet, TOI-125d, has an orbital period of 19.98 d and eccentricity |$e_\mathrm{ d} = 0.168^{+0.088}_{-0.062}$|⁠. With a radius of 2.93 ± 0.17 R_E and mass 13.6 ± 1.2 M_E, it is the densest of the three planets, with ρ_P = 2.98 g cm⁻³.

TOI-125b, TOI-125c, and TOI-125d are thus all mini-Neptunes with similar radii, but different masses yielding a high–low–higher density pattern outwards in the system. The planets straddle the gap identified in the mass-period plane by Armstrong et al. (2019). All three planets have the same orbital inclination to within a degree. The high orbital eccentricities detected for TOI-125b and d are unusual for such a compact system of mini-Neptunes (Van Eylen et al. 2019).

TOI-125 is found to be a main-sequence K0-star with a mass |$0.859^{+0.044}_{-0.038}$| M_⊙, radius 0.848 ± 0.011 R_⊙, and T_eff = 5320 ± 39 K. This is in reasonable agreement with the properties reported in the Gaia Data Release 2: R_* = 0.90 ± 0.03 R_⊙ and T_eff = 5150 ± 84 K (Gaia Collaboration 2018). The quadratic drift found in the RVs might indicate the existence of an additional massive companion in the system, at a long period P ≳ 100 d. We obtained a few RV points in 2019 July with low precision to rule out a stellar companion. More high-precision RVs would be needed to determine the nature of this long-term signal.

3.3.1 Marginal planet candidates TOI-125.04 and TOI-125.05

Fig. 4 shows the residuals from the three-planet fit. We see no hint of any signal from TOI-125.04 or TOI-125.05. The strong peak at P = 0.49 d is an alias of the residual signal at 50 d.

Figure 4.

Periodogram of the RV residuals after fitting three planets with eccentric orbits and a quadratic trend. The horizontal black line is the 1 per cent FAP.

We derive upper mass limits for the two-planet candidates by running EXOFASTv2 on the HARPS RVs while only including priors on the orbital period and transit depth from Quinn et al. (2019). We do not include the TESS photometry, to save computational time. Fitting three, four, or five planets has little impact on the final parameters for TOI-125b, TOI-125c, and TOI-125d. For the marginal USP candidate TOI-125.04 (P = 0.53 d, |$R_P = 1.36^{+0.14}_{-0.16}$|R_E), we find an RV semi-amplitude of |$K=0.56^{+0.4}_{-0.3}$| m s⁻¹ corresponding to a 2σ upper mass limit of 1.6 M_E. Our measurement is compatible with no planet and we cannot validate this candidate. The highest bulk density allowed by the data (based on the upper mass limit and 1σ lower radius 1.20 R_E) is ρ_{P, max} = 5.10 g cm⁻³. For highly irradiated super-Earth candidates such as TOI-125.04, we expect highly irradiated rocky cores with high densities. More observations either with an HARPS-like or more precise instrument such as ESPRESSO (Pepe et al. 2010) would be required to confirm the existence and mass of TOI-125.04.

For TOI-125.05 (P = 13.28 d), we find an RV semi-amplitude consistent with zero, |$K=0.2^{+0.4}_{-0.18}$| m s⁻¹ corresponding to a 2σ upper mass limit of 2.7 M_E. The posterior distribution for the planetary radius presented by Quinn et al. (2019) is bi-modal and peaks at 4.2 and 13.5 R_E. The 1σ median for the whole distribution is |$8.8^{+4.7}_{-4.4}$|R_E, which does not reflect the true nature of the posterior. The RV data presented by Quinn et al. (2019) and this study both exclude the upper part of the distribution, meaning that if the planet is real its radius will most likely be similar to that of TOI-125b, TOI-125c, and TOI-125d. We thus only consider the lower part of the radius posterior distribution with 68 per cent confidence intervals |$4.2^{+2.2}_{-1.4}$|R_E. The highest bulk density allowed by the data (based on the upper mass limit and 1σ lower radius 2.8 R_E) is ρ_{P, max} = 0.38 g cm⁻³. This is a very low density close to being unphysical for a mini-Neptune. We thus conclude that TOI-125.05 is unlikely a viable planet candidate.

4 DYNAMICAL STABILITY AND SYSTEM ARCHITECTURE

The period ratios in the TOI-125 system are interesting. If we assume that the low-SNR USP candidate TOI-125.04 is a planet, then the orbital period ratios of adjacent pairs are (beginning from the outside) 2.183, 1.966, and 8.806. The period ratio between planet d and planet c, 2.183, lies at the second most prominent peak in the period ratio distribution (Lissauer et al. 2011; Fabrycky et al. 2014; Steffen & Hwang 2015) of known exoplanets, close to the 2:1 MMR. The origin of this peak is unknown, though it appears both in systems with known intermediate planets (as we see here) and in systems with no observed intermediate planets.

Next, between planets c and b the period ratio is 1.966 – sufficiently interior to the 2:1 MMR to be consistent with the observed gap in planet pairs interior to such resonances (Lissauer et al. 2011; Fabrycky et al. 2014). There are multiple explanations for this gap interior to the first-order MMR, though none have been demonstrated as the primary cause (Delisle et al. 2012; Lithwick & Wu 2012; Rein 2012; Batygin & Morbidelli 2013; Petrovich, Malhotra & Tremaine 2013; Chatterjee & Ford 2015). Further study of systems like TOI-125 may shed additional light on its origin. Finally, the innermost planet candidate has an orbital period that is less than 1 d. With its neighbour this pair has the largest period ratio in the system. This is consistent with the observed trend that when one member of an adjacent pair of planets has an orbital period less than 1 d, the period ratio is unusually large (Steffen & Farr 2013; Sanchis-Ojeda et al. 2014; Steffen & Coughlin 2016). The origin of the USP planets remains unknown (Winn, Sanchis-Ojeda & Rappaport 2018) though a number of hypotheses have been proposed ranging from stripped cores of giant planets (Valsecchi, Rasio & Steffen 2014; Königl, Giacalone & Matsakos 2017) to various dynamical effects coupled with stellar tides (Muñoz, Lai & Liu 2016; Lee & Chiang 2017; Petrovich, Deibert & Wu 2019; Pu & Lai 2019). The nearby presence of additional small planets would seem not to support the stripped-cores possibility, since hot Jupiter planets tend to be alone with few exceptions (Wright et al. 2009; Steffen et al. 2012; Becker et al. 2015). Moreover, (Winn et al. 2017) showed that the metallicity trends of these USP planets do not match those of hot Jupiters – implying that if USP planets are stripped cores, they must be from smaller, sub-Neptune planets.

The masses of the planets are sufficiently large that in situ formation is unlikely (see e.g. Schlichting 2014). Thus, formation at larger distances in a protoplanetary disc and migration inwards are a possibility. Planets in resonance are a clear indication of planet migration. Furthermore, if the planets formed in the same location in a protoplanetary disc, it would be expected that they would have formed out of similar disc material and thus have similar densities. The fact that neighbouring planets have significantly different densities is also indicative that they formed in different locations and migrated inwards, as investigated for the Kepler 36 system (Carter et al. 2012) by Bodenheimer et al. (2018) and Raymond et al. (2018).

4.1 N-body simulations

We attempted to refine the orbital parameters and planet masses for the TOI-125 system by requiring the system parameters to be compatible with dynamical stability. For this purpose, we considered the three-planet model for TOI-125,³ as illustrated in Table 3. We used several thousand draws uniformly selected over the full EXOFASTv2 MCMC posterior as sets of initial conditions.

Each set was integrated over a time span of 5000 yr, corresponding to approximately 91 000 revolutions of the outer planet TOI-125d. The simulations were performed with an adaptive time-stepping using the N-body 15th-order integrator IAS15 (Rein & Spiegel 2015), available from the software package REBOUND⁴ (Rein & Liu 2012). The general relativity correction was included following Anderson et al. (1975), via the python module REBOUNDx. Then, the stability of each system was explored using the NAFF chaos indicator (Laskar 1990, 1993). The latter consists in estimating precisely the average of the mean motion n of each planet over the first half of the simulation, and repeating this procedure over the second half. The bigger the variation in this average, the more chaotic the system is. Most often, this leads to escapes or close encounters between bodies, defining the system as unstable. Finally, we define a new posterior distribution by keeping only the stable systems. Linking the MCMC exploration of the parameter space with fast chaos indicators is particularly efficient (Stalport et al., in preparation).

The coupled photometric and RV observations give constraints on all the orbital parameters except the longitudes of the nodes of the planets Ω. As a result, this parameter is absent from the EXOFASTv2 MCMC posterior. Therefore, we performed a first series of 5000 numerical simulations in which the initial values for the Ω parameters were selected randomly from a uniform distribution between −π and π. The new, dynamically stable posterior distribution strikingly selects only the systems in which the planets have aligned or anti-aligned lines of nodes. This result is illustrated in Fig. 5. It is explained by the fact that, in these configurations, the mutual inclinations between the adjacent planets are minimal.⁵ Let us note that no information is provided regarding the individual value of Ω for each planet. However, the dynamical constraints allow us to state that Ω_k − Ω_j = 0 or π, for j and k denoting the planets.

Figure 5.

Dynamically stable posterior distribution projected on to Ω_c − Ω_b in green, and Ω_d − Ω_c in blue. The peaks at around 0 and ±180 deg strongly favour the aligned or anti-aligned configurations for the lines of nodes of the planets.

Projected on to the other orbital parameters and planetary masses, the dynamically stable posterior distribution does not bring more information. It mimics the original MCMC posterior distribution. This poor refinement can be explained by the aforementioned observation about the lines of nodes. Indeed, many systems turned out to be unstable only because of the unfavourable configurations given by Ω, and the real constraints on the observations were hidden.

To overcome this bias, we launched a second set of 10 000 numerical simulations. This time, the longitudes of the nodes of the planets were selected randomly in windows around the alignment or anti-alignment, as illustrated by the vertical lines on Fig. 5. An interesting result of this process is shown in Fig. 6. The posterior distribution is projected on to the plane of two parameters, the eccentricity and argument of periastron of the outer planet (e_d and ω_d). As seen in the figure, a branch of solutions at ω_d ∼ 60° explores high values of e_d. However, this region is disfavoured, as expressed by the decrease in the median of e_d.

Figure 6.

Sample of the posterior distribution from the EXOFASTv2 MCMC of the three-planet model, projected on the parameters e_d and ω_d. In red, the full sample is projected. The dots are coloured in blue if the corresponding systems are qualified as stable by the NAFF indicator. The black horizontal lines denote the median values of the distributions of e_d. The dashed line is associated with the full sample, while the plain line corresponds to the dynamically stable sample. The same applies for ω_d and the vertical lines.

Another result concerns the relatively high eccentricity of the inner planet, which has a best-fitting value of e_b ∼ 0.194. In Fig. 7, we show the posterior distribution projected on to this parameter in red. The observations are inconsistent with zero eccentricity. A slight displacement towards lower eccentricities is observed in the dynamically stable distribution. Indeed, with the stability constraint, the median of the distribution shifted from med(e_b) ∼ 0.188 (red histogram) to med(e_b) ∼ 0.177 (blue histogram). However, many systems with large eccentricities remain stable. Therefore, such large eccentricities do not seem incompatible with stability.

Figure 7.

Similar plot as Fig. 6. The posterior distribution is now projected on to the single parameter e_b.

4.2 Tidal interactions

The high eccentricity of planet b also raises questions concerning the tidal evolution of the system. To investigate those aspects, we also performed N-body integrations taking into account the tidal forces and torques. To perform those simulations, we used Posidonius⁶ (Blanco-Cuaresma & Bolmont 2017), which allows us to take into account tides, as well as rotational flattening and general relativity using the same prescriptions as in Bolmont et al. (2015).

For tides, Posidonius uses an equilibrium tide model (Mignard 1979; Hut 1981; Eggleton, Kiseleva & Hut 1998), for which the tidal dissipation of the different bodies is quantified by the product k₂Δτ of the constant time lag Δτ and the Love number of degree 2 k₂ (the bigger this quantity, the bigger the dissipation and the faster the evolution). As the underlying assumption of this constant time lag model is that the planet is made of a weakly viscous fluid, it is appropriate for the low-density planets of TOI-125. We use a constant time lag similar to Jupiter’s (k₂Δτ ∼ 2.5 × 10⁻² s from Leconte et al. 2010) and explore a range between 1 and 10² times this value.

Assuming this dissipation for all planets leads to very long evolution time-scales. In particular, the time-scale of circularization for planet b is about ≳10¹⁰ yr and it reaches 10¹³ yr for planet d, which is much higher than the estimated age of ∼7 Gyr. The high eccentricities are therefore not completely surprising and the fact that planet b has not circularized also puts constraints on its dissipation: it cannot be much higher than Jupiter’s. However, the time-scales for the damping of the planetary obliquity (angle between the rotation axis and the perpendicular to the orbital plane) and of synchronization are shorter. Assuming the same dissipation as Jupiter, and even assuming the lower estimate of the age (2.5 Gyr), we find that planets b and c should have a damped obliquity (less than a few degrees) and an evolved rotation. In our model, the evolved rotation period is the pseudo-synchronization period, which depends on eccentricity (Hut 1981). Depending on the age of the system, the obliquity and rotation of planet d might still be evolving: If the system is older than ∼6 Gyr, the obliquity should be very small and the rotation should be very close to the pseudo-synchronization rotation.

Of course, there is a strong uncertainty on the dissipation factor of planets, these planets could dissipate more energy than what is estimated for Jupiter (with processes such as tidal inertial waves in the convective region; Ogilvie & Lin 2004). But unless the age of the system is close to its upper estimate of 11 Gyr, the fact that planet b still has a high eccentricity tends to indicate that dynamical tide processes are not very efficient.

5 INTERNAL STRUCTURE

In order to characterize the internal structure of TOI-125b, TOI-125c, and TOI-125d, we construct models considering a pure-iron core, a silicate mantle, a pure-water layer and an H–He atmosphere. The models follow the basic structure model of Dorn, Hinkel & Venturini (2017), with the equation of state (EOS) for the iron core taken from Hakim et al. (2018), and the EOS of the silicate mantle from Connolly (2009). For water, we use the quotidian EOS of Vazan et al. (2013) for low pressures and the one of Seager et al. (2007) for pressures above 44.3 GPa. The hydrogen–helium (H–He) EOS is SCVH (Saumon, Chabrier & van Horn 1995) assuming a proto-solar composition. We then use a generalized Bayesian inference analysis using a Nested Sampling scheme (e.g. Buchner 2016). We then quantify the degeneracy between interior parameters and produce posterior probability distributions. The interior parameters that are inferred include the masses of the pure-iron core, silicate mantle, water layer, and H–He atmospheres. For this analysis, we use the stellar Fe/Si and Mg/Si ratios from Table 2 as a proxy for the planet abundances.

Table 2.

Spectral parameters derived from the stacked HARPS spectrum with SNR/resolution element ∼500 at 5500 Å, using two different methods. T_eff and [Fe/H] and their uncertainties were used as Gaussian priors on the MCMC joint modelling of the planetary and stellar parameters – we used the average between the two approaches. The errors were inflated to encompass both values at a 1σ level, in order to reflect the model dependence of the atmospheric parameters. V_t denotes micro- and macro-turbulence velocities.

	Equivalent width		SME
Parameter	Value	1σ	Value	1σ
Teff(K)	5295	42	5125	60
log g (cgs)	4.51	0.07	4.4	0.2
Vtmicro (km s−1)	0.72	0.09	0.8	0.1
Vtmacro (km s−1)	–	–	2.5	0.5
v sin i (km s−1)	–	–	1.0	0.5
[Fe/H](dex)	−0.02	0.03	0.00	0.05
Na i/H (dex)	−0.06	0.05	−0.1	0.05
Mg i/H (dex)	0.01	0.05	–	–
Al i/H (dex)	−0.02	0.07	–	–
Si i/H (dex)	−0.04	0.06	−0.1	0.05
Ca i/H (dex)	−0.03	0.07	−0.1	0.05
Sc ii/H (dex)	−0.02	0.04	–	–
Ti i/H (dex)	0.09	0.06	−0.05	0.05
Cr i/H (dex)	0.02	0.05	0.0	0.05
Ni i/H (dex)	−0.06	0.03	−0.05	0.05
Zr/H (dex)	–	–	−0.1	0.05

Table 3.

Median values and 68 per cent confidence intervals for TOI-125b, TOI-125c, and TOI-125d and their host star fitted with EXOFASTv2, while including a quadratic RV drift and orbital eccentricities for all three planets.

Stellar parameters:
M*	Mass (M⊙)	|$0.859^{+0.044}_{-0.038}$|	–	–
R*	Radius (R⊙)	0.848 ± 0.011	–	–
L*	Luminosity (L⊙)	0.519 ± 0.016	–	–
ρ*	Density (cgs)	|$1.99^{+0.13}_{-0.11}$|	–	–
log g	Surface gravity (cgs)	4.516 ± 0.024	–	–
Teff	Effective temperature (K)	5320 ± 39	–	–
[Fe/H]	Metallicity (dex)	−0.02 ± 0.03	–	–
Age	Age (Gyr)	|$6.8^{+4.4}_{-4.1}$|	–	–
AV	V-band extinction (mag)	|$0.032^{+0.032}_{-0.023}$|	–	–
d	Distance (pc)	111.40 ± 0.44	–	–
|$\dot{\gamma }$|	RV slope (m s−1 d−1)	−0.0123 ± 0.0078	–	–
|$\ddot{\gamma }$|	RV quadratic term (m s−1 d−2)	−0.001 83 ± 0.000 25	–	–
Planetary parameters:		b	c	d
P	Period (days)	|$4.653\,82^{+0.000\,33}_{-0.000\,31}$|	|$9.150\,59^{+0.000\,70}_{-0.000\,82}$|	|$19.9800^{+0.0050}_{-0.0056}$|
RP	Radius (RE)	2.726 ± 0.075	2.759 ± 0.10	2.93 ± 0.17
MP	Mass (ME)	9.50 ± 0.88	6.63 ± 0.99	13.6 ± 1.2
ρP	Density (cgs)	2.57 ± 0.33	1.73 ± 0.33	|$2.98^{+0.65}_{-0.52}$|
TC	Time of conjunction (⁠|$\rm {BJD_{TDB}}$|⁠)	58 355.355 29 ± 0.0010	58 361.9085 ± 0.0013	58 342.8516 ± 0.0039
a	Semimajor axis (au)	|$0.051\,86^{+0.000\,86}_{-0.000\,77}$|	0.0814 ± 0.0013	0.1370 ± 0.0022
b	Transit impact parameter	|$0.27^{+0.17}_{-0.18}$|	|$0.522^{+0.086}_{-0.18}$|	|$0.652^{+0.093}_{-0.16}$|
i	Inclination (degrees)	|$88.92^{+0.71}_{-0.60}$|	|$88.54^{+0.41}_{-0.19}$|	|$88.795^{+0.18}_{-0.10}$|
e	Eccentricitya	|$0.194^{+0.041}_{-0.036}$|	|$0.066^{+0.070}_{-0.047}$|	|$0.168^{+0.088}_{-0.062}$|
ω*	Argument of periastron (degrees)	|$-37^{+12}_{-14}$|	|$70^{+100}_{-110}$|	|$46^{+23}_{-44}$|
Teq	Equilibrium temperature (K)	1037 ± 11	827.8 ± 8.6	638.1 ± 6.6
〈F〉	Incident flux (109 erg s−1 cm−2)	0.252 ± 0.012	0.1056 ± 0.0045	0.0363 ± 0.0019
K	RV semi-amplitude (m s−1)	4.11 ± 0.36	2.25 ± 0.33	3.61 ± 0.31
RP/R*	Radius of planet in stellar radii	0.029 50 ± 0.000 70	0.029 85 ± 0.000 99	0.0317 ± 0.0018
a/R*	Semimajor axis in stellar radii	13.16 ± 0.27	20.66 ± 0.42	34.770.70
δ	Transit depth (fraction)	|$0.000\,870^{+0.000\,043}_{-0.000\,040}$|	|$0.000\,891^{+0.000\,060}_{-0.000\,057}$|	|$0.001\,00^{+0.000\,12}_{-0.000\,11}$|
τ	Ingress/egress duration (days)	|$0.003\,80^{+0.000\,61}_{-0.000\,26}$|	|$0.004\,86^{+0.000\,79}_{-0.000\,93}$|	|$0.0068^{+0.0021}_{-0.0017}$|
T14	Total transit duration (days)	0.1234 ± 0.0024	|$0.1231^{+0.0026}_{-0.0030}$|	|$0.1297^{+0.0070}_{-0.0057}$|
TFWHM	FWHM transit duration (days)	|$0.1194^{+0.0023}_{-0.0024}$|	|$0.1182^{+0.0027}_{-0.0031}$|	|$0.1227^{+0.0076}_{-0.0062}$|
TP	Time of periastron (⁠|$\rm {BJD_{TDB}}$|⁠)	|$58\,326.03^{+0.17}_{-0.20}$|	|$58\,334.1^{+2.3}_{-2.8}$|	|$58\,341.1^{+1.1}_{-2.4}$|
TS	Time of eclipse (⁠|$\rm {BJD_{TDB}}$|⁠)	|$58\,325.546^{+0.084}_{-0.081}$|	|$58\,339.06^{+0.30}_{-0.24}$|	|$58\,334.18^{+0.54}_{-0.57}$|
log gP	Surface gravity	3.097 ± 0.047	|$2.931^{+0.068}_{-0.076}$|	3.192 ± 0.064
Θ	Safronov number	|$0.0148^{+0.0014}_{-0.0013}$|	0.0160 ± 0.0024	|$0.0522^{+0.0054}_{-0.0051}$|
Wavelength parameters:		TESS
u1	Linear limb-darkening coeff	0.382 ± 0.035	–	–
u2	Quadratic limb-darkening coeff	|$0.240^{+0.035}_{-0.036}$|	–	–
Telescope parameters:		HARPS	–	–
γrel	Relative RV offset (m s−1)	11 441.90 ± 0.30	–	–
σJ	RV jitter (m s−1)	|$1.63^{+0.24}_{-0.22}$|	–	–
Transit parameters:		TESS Sector 1	TESS Sector 2	–
σ2	Added variance	−0.000 000 023 ± 0.000 000 027	−0.000 000 046 ± 0.000 000 027	–
F0	Baseline flux	1.000 136 ± 0.000 019	1.000 151 ± 0.000 019	–

Note. ^aThe eccentricities presented here are the direct outputs from EXOFASTv2, without any constraints from N-body simulations. Our dynamical analysis in Section 4.1 puts upper limits on the eccentricities for TOI-125b and TOI-125c, but retains the same eccentricities within a 1σ confidence interval.

Fig. 8 shows the mass–radius relation for a pure-water curve and a planet with 95 per cent water and 5 per cent H–He atmosphere subjected to a stellar radiation of F/F_⊕ = 100 (comparable to the case of the TOI-125 planets). All three planets could, in principle, either consist of a rocky core with a massive water envelope (mostly in the form of supercritical steam) or a rocky core with a likely high-metallicity H–He envelope (up to 5 per cent in mass of H–He). The position of the three planets in the insolation radius diagram (Fig. 9), above the evaporation valley (Fulton et al. 2017; Fulton & Petigura 2018; Van Eylen et al. 2018), indicates, however, that the latter scenario (i.e. involving an H2/He envelope) is the most plausible one (Owen & Wu 2017; Ginzburg, Schlichting & Sari 2018). Spectroscopic transit measurements will hopefully help to discriminate between the two aforementioned cases owing to the relative proximity of the TOI-125 system; see Section 6 for a more in-depth discussion. Transit observations of the exoplanet GJ1214b – which lies in a somewhat similar insolation radius–mass parameter space than the TOI-125 planets – have, however, shown that clouds may limit our ability to conclude on the true nature of these objects (Kreidberg et al. 2014).

Figure 8.

Mass–radius diagram of exoplanets with accurate mass and radius determination (Otegi et al. 2019). Also shown are the composition lines of an Earth-like planet, pure water, and 95 per cent |$\mathrm{H_2 O+5{{\ \rm per\ cent}} \text{ }H\!-\!He}$|⁠.

Figure 9.

Insolation flux relative to the Earth plotted against radii for known exoplanets extracted from the NASA Exoplanet Archive, as presented in Fulton et al. (2017) and Fulton & Petigura (2018). The orange contours indicate point density (not occurrence), showing the separate populations of mini-Neptunes and super-Earths. TOI-125b, TOI-125c, and TOI-125d are plotted as three stars in the same colours as in Figs 1, 3, and 8.

Table 4 lists the inferred mass fractions of the core, mantle, water layer, and H–He atmosphere from our structure models. We find median H–He mass fractions of 2.3 per cent for TOI-125b, 2.9 per cent for TOI-125c, and 4.5 per cent for TOI-125d. These estimates are lower bounds since structure models considering H–He envelopes enriched with heavy elements could result in even higher values. This is because enriched H–He atmospheres are more compressed, and can therefore increase the planetary H–He mass fraction. Indeed, formation models of mini-Neptunes suggest that forming such planets without envelope enrichment is very unlikely (Venturini & Helled 2017).

Table 4.

Inferred interior structure properties of TOI-125b, TOI-125c, and TOI-125d.

Interior structure	TOI-125b	TOI-125c	TOI-125d
Mcore/Mtotal	|$0.31 ^{+0.18} _{-0.32}$|	|$0.31 ^{+0.16} _{-0.27}$|	|$0.26 ^{+0.16} _{-0.21}$|
Mmantle/Mtotal	|$0.39 ^{+0.17} _{-0.26}$|	|$0.38 ^{+0.18} _{-0.29}$|	|$0.36 ^{+0.18} _{-0.31}$|
Mwater/Mtotal	|$0.32 ^{+0.20} _{-0.24}$|	|$0.32 ^{+0.17} _{-0.24}$|	|$0.36 ^{+0.16} _{-0.21}$|
MH–He/Mtotal	|$0.020 ^{+0.006} _{-0.008}$|	|$0.027 ^{+0.007} _{-0.010}$|	|$0.041 ^{+0.009} _{-0.012}$|

TOI-125b and TOI-125c are expected to have very similar compositions, with core and water layer mass fractions of |$\text{$\sim$} 30{{\ \rm per\ cent}}$| and a mantle mass fraction of |$\text{$\sim$} 40{{\ \rm per\ cent}}$|⁠. TOI-125d, instead, has a slightly higher water mass fraction of 35 per cent, and a smaller fraction of refractory materials with a core mass fraction of 26 per cent and mantle mass fraction of 35 per cent.

6 POTENTIAL FOR ATMOSPHERIC CHARACTERIZATION

Our analysis of the internal structure (see Table 4), as well as the position of the three planets in the insolation radius diagram (see Fig. 9), indicates that all the three planets might have a water-dominated atmosphere with a small contribution from lighter elements of the order of a few per cent. If these light elements are evaporated over time (especially for TOI-125b, the most irradiated in the system), their observation could be used to study the planets’ exospheres.

Due to the significant distance of the system (111.40 pc), the absorption of the interstellar medium (ISM) puts Lyman-α observations out of reach. However, H-alpha and He i, which do not suffer from ISM absorption, can be used to detect a potential escaping planetary outflow. H-alpha and other Balmer series lines have been detected for several exoplanets, showing deep absorption features observed at high spectral resolution (Jensen et al. 2012; Cauley, Redfield & Jensen 2017; Jensen et al. 2018; Yan & Henning 2018). Likewise, the well-known He i triplet in the infrared (Seager & Sasselov 2000; Oklopčić & Hirata 2018; Oklopčić 2019) has also successfully detected exospheric absorption in other systems (Allart et al. 2018; Nortmann et al. 2018; Salz et al. 2018; Allart et al. 2019).

The possible water-rich composition from Table 4 could be verified via observations in the infrared, and thus provides valuable insights into the water composition in a multiplanet system with three similarly sized planets but different masses and insolations. However, observations from the ground are challenging due to the planets’ sizes and observational windows. We estimated that one transit observations would not be useful to detect water bands for TOI-125b (scale height 38 km) with NIRPS at the ESO 3.6 m telescope (Bouchy et al. 2017). Observing multiple transits would require a dedicated large program spanning several years given the possible observational windows from Chile. It is, however, a prime target for observations with the next generation of ELTs, particularly with the HIRES optical-to-NIR spectrograph at the E-ELT (Marconi et al. 2016) and CRIRES + at the VLT (Follert et al. 2014).

Using the Pandexo Exposure Time Calculator for HST,⁷ we estimate that the precision with which we can measure the transmission spectrum of TOI-125b using the Wide-Field Camera 3 (WFC3) instrument, in five transits, is ∼30 ppm near the 1.4 μm water feature. The expected water signature at 5-scale heights has a depth of approximately 20 ppm; thus, detecting this feature with HST would be challenging for a planet with an atmosphere as compact as TOI-125b. However, all three planets are prime targets for the JWST’s NIRSpec.

7 CONCLUSIONS

We confirm the detection of three mini-Neptunes around TOI-125 found by TESS using HARPS RV measurements. TOI-125b, TOI-125c, and TOI-125d all have similar radii, 2.726 ± 0.075 R_E, 2.759 ± 0.10 R_E, and 2.93 ± 0.17 R_E, respectively. The three planets differ greatly in mass, however, with 9.50 ± 0.88 M_E, 6.63 ± 0.99 M_E, and 13.6 ± 1.2 M_E, yielding a high–low–higher pattern in terms of density when moving outward in the system. For the two marginal planet candidates TOI-125.04 and TOI-125.05, we derive 2σ upper mass limits of 1.6 M_E and 2.7 M_E, respectively. For TOI-125.05, this means that it is unlikely a viable planet candidate.

The system exhibits an intriguing architecture with the two inner planets slightly interior to the 2:1 MMR while the two outer planets are slightly external to the 2:1 MMR. TOI-125b and TOI-125d both show significant orbital eccentricities. We analyse the dynamics of the system using N-body simulations and demonstrate that planetary orbits are stable despite the high eccentricities. Based on N-body simulations coupled with tidal forces and torques, we conclude that the dynamical tide processes cannot be very efficient in order for TOI-125b to retain its high eccentricity of |$e_\mathrm{ b} = 0.194 ^{+0.041}_{-0.036}$|⁠.

Our analysis of the internal compositions of these three planets yields that they all most likely retain H–He atmospheres and a significant water layer, which could be detected through transmission spectroscopy. This is expected for planets sitting on top of the radius gap (see Fig. 9), receiving less than 300 times the stellar insolation than that of the Earth.

ACKNOWLEDGEMENTS

We thank the anonymous referee for providing thoughtful comments that allowed us to improve on this paper. This study is based on observations collected at the European Southern Observatory under ESO programmes 0101.C-0829, 1102.C-0249, 1102.C-0923, 0102.C-0525, and 0102.C-0451. We thank the Swiss National Science Foundation (SNSF) and the Geneva University for their continuous support to our planet search programs. This work has been in particular carried out in the framework of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). This publication makes use of DACE, which is a facility based at the University of Geneva (CH) dedicated to extrasolar planets data visualization, exchange, and analysis. DACE is a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS, federating the Swiss expertise in exoplanet research. The DACE platform is available at https://dace.unige.ch. This paper includes data collected by the TESS mission. Funding for the TESS mission is provided by the NASA Explorer Program. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. 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. DJA acknowledges support from the STFC via an Ernest Rutherford Fellowship (ST/R00384X/1). The IA/Portuguese team was supported by FCT/MCTES through national funds and by FEDER – Fundo Europeu de Desenvolvimento Regional through COMPETE2020 – Programa Operacional Competitividade e Internacionalização by these grants: UID/FIS/04434/2019; PTDC/FIS-AST/32113/2017 and POCI-01-0145-FEDER-032113; PTDC/FIS-AST/28953/2017 and POCI-01-0145-FEDER-028953. VA acknowledges the support from FCT through Investigador FCT contract no. IF/00650/2015/CP1273/CT0001. SH acknowledges support by the fellowships PD/BD/128119/2016 funded by FCT (Portugal). SCCB acknowledges support from FCT through Investigador FCT contracts IF/01312/2014/CP1215/CT0004. ODSD acknowledges the support from FCT (Portugal) through work contract DL 57/2016/CP1364/CT0004. MRD acknowledges support of CONICYT-PFCHA/Doctorado Nacional-21140646 and Proyecto Basal AFB-170002. JSJ acknowledges support from FONDECYT grant 1161218. FM acknowledges support from The Royal Society Dorothy Hodgkin Fellowship. JVS and LAdS are supported by funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project Four Aces; grant agreement no. 724427). DJAB acknowledges support from the UK Space Agency. JNW acknowledges support from the Heising-Simons Foundation. NN is supported by JSPS KAKENHI Grant Numbers JP18H01265 and JP18H05439, and JST PRESTO Grant Number JPMJPR1775. CD acknowledges support from the Swiss National Science Foundation under grant PZ00P2_174028. KWFL acknowledges support by DFG grants RA714/14-1 within the DFG Schwerpunkt SPP 1992, ‘Exploring the Diversity of Extrasolar Planets’. DB and JLB have been funded by the Spanish State Research Agency (AEI) Project No. ESP2017-87676-C5-1-R and No. MDM-2017-0737 Unidad de Excelencia María de Maeztu Centro de Astrobiología (CSIC-INTA). SM acknowledges support from the Spanish Ministry under the Ramon y Cajal fellowship number RYC-2015-17697. RAG acknowledges the support from PLATO and GOLF CNES grants. This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Grant Agreement No. 832738/ESCAPE. MT acknowledges funding from the Gruber Foundation. MF, IG, and CMP gratefully acknowledge the support of the Swedish National Space Agency (DNR 163/16 and 174/18).

Footnotes

1

https://dace.unige.ch/radialVelocities/?pattern = TOI-125

2

See HARPS instrument monitoring pages: https://www.eso.org/sci/facilities/lasilla/instruments/harps/inst/monitoring/thar_history.html.

3

If real, the USP candidate TOI-125.04 is not expected to play a significant dynamical role in the system, due to its large period ratio with TOI-125b.

4

The REBOUND code is freely available at http://github.com/hannorein/rebound.

5

The mutual inclination I_m between two orbits is a quantity that depends on the inclination of each orbit i_k and i_j with respect to the plane of the sky, and on the difference in the longitudes of the nodes ΔΩ = Ω_k − Ω_j (j and k denote the planets). Its expression is cos I_m = cos i_k cos i_j + cos ΔΩ sin i_k sin i_j.

6

The Posidonius code is freely available at https://github.com/marblestation/posidonius.

7

Available at https://exoctk.stsci.edu/pandexo/.

Author notes

NASA Sagan Fellow.