Masses, revised radii, and a third planet candidate in the ‘Inverted’ planetary system around TOI-1266

ABSTRACT

1 INTRODUCTION

The radius valley is a stark feature in the exoplanet population between |$\sim 1.6\,{\text {and}}\,1.9\, \mathrm{R}_{\oplus }$| that separates terrestrial super-Earths from larger sub-Neptunes, which are enveloped in a low-density volatile material such as H/He or water. The radius valley exists around FGK (e.g. Fulton et al. 2017; Fulton & Petigura 2018; Van Eylen et al. 2018; Berger et al. 2020; Hardegree-Ullman et al. 2020; Petigura et al. 2022) and M dwarfs alike (e.g. Cloutier & Menou 2020; Hsu, Ford & Terrien 2020; Van Eylen et al. 2021), albeit with distinct dependencies on planet radius and instellation that are suggestive of different radius valley emergence mechanisms around different spectral types.

Proposed radius valley emergence models include thermally driven mass-loss (TDML) of planetary atmospheres, which may be driven by XUV photons from the host star (e.g. photoevaporation; Owen & Wu 2013; Jin et al. 2014; Lopez & Fortney 2014) or from heating by the planetary core after formation (e.g. core-powered mass-loss; Ginzburg, Schlichting & Sari 2018; Gupta & Schlichting 2019). Alternatively, the radius valley may emerge directly from the planet formation process via either late-stage super-Earth formation in a gas-depleted environment (Lopez & Rice 2018), by limited gas accretion on to low-mass cores (Lee & Chiang 2015; Lee & Connors 2021), or by water-rich formation (WRF) followed by inward migration (Raymond et al. 2018; Venturini et al. 2020; Burn et al. 2021).

A number of lines of recent empirical evidence are suggesting that the M dwarf radius valley likely emerges directly from the planet formation process without the need to invoke atmospheric escape. Planet occurrence rates show evidence that the slope of the M dwarf radius valley in radius-instellation space is inconsistent with TDML models and instead favours a gas-depleted formation (GDF) scenario (Cloutier & Menou 2020). The XUV spectrum of the multiplanet host K2-3 cannot explain the low-density compositions of its inner planets using TDML models (Diamond-Lowe et al. 2022). A subpopulation of water-rich planets around M dwarfs may have been uncovered following a uniform reanalysis of archival radial velocity time series (Luque & Pallé 2022), although this claim has been contended (Rogers, Schlichting & Owen 2023). The low bulk density of the intensively studied |$1.5\, \mathrm{R}_{\oplus }$| planet Kepler-138 d was shown to be inconsistent with TDML and requires a substantial volatile mass fraction suggestive of WRF beyond the snowline (Piaulet et al. 2023). Lastly, the bulk compositions of seven keystone planets (i.e. within the M dwarf radius valley) require substantial volatile mass fractions and thus disfavour a TDML explanation (Cherubim et al. 2023). This collection of results paints the emerging picture that unlike around FGK stars, close-in planets around M dwarfs are unlikely to be sculpted by TDML, and instead, we may be witnessing the emergence of a new channel of small planet formation.

TOI-1266 is a validated planetary system from NASA’s Transiting Exoplanet Survey Satellite mission that features two small planets purportedly sitting on opposing sides of the M dwarf radius valley (Demory et al. 2020; Stefánsson et al. 2020, hereafter D20 and S20, respectively). The system was reported to harbour an inner sub-Neptune (P_b = 10.9 d, |$R_{p,b}=2.37\, \mathrm{R}_{\oplus }$|⁠) alongside an outer super-Earth (P_c = 18.8 d, |$R_{p,c}=1.56\, \mathrm{R}_{\oplus }$|⁠; D20). Multiplanet system architectures that contain a likely terrestrial planet (i.e. |$\lesssim 1.6\, \mathrm{R}_{\oplus }$|⁠) whose orbit is wider than the system’s sub-Neptune are somewhat rare (∼35 per cent; Weiss et al. 2018). Most planet pairs that span the radius valley exhibit the super-Earth on an orbit that is interior to that of the sub-Neptune. Seemingly ‘inverted’ architectures like TOI-1266 are particularly interesting because they can be challenging to reconcile with models of atmospheric escape and population synthesis models based on core accretion theory (Owen & Wu 2017; Gupta & Schlichting 2019; Burn et al. 2021). TOI-1266 is therefore an important testbed for planet formation models, particularly the emergence of the M dwarf radius valley. This motivates our study to measure the planet masses and bulk compositions, and test the consistency of our findings with various radius valley emergence models.

In this paper, we present the results of our campaign to measure the masses of TOI-1266 b and c using precise radial velocity measurements taken with the HARPS-N spectrograph. In Section 2, we describe our observations. In Section 3, we summarize our knowledge of the host star. In Section 4, we present our data analysis. In Section 5, we discuss our results and their implications for sub-Neptune formation around M dwarfs before concluding with a summary of our key findings in Section 6.

2 OBSERVATIONS

2.1 TESS transit photometry

TOI-1266 b and c were originally discovered in the primary mission (PM) of NASA’s Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015), as reported by D20 and S20. Both discovery papers presented data from TESS sectors 14, 15, 21, and 22, which spanned ut 2019 July 18 to 2020 March 17. TOI-1266 has since been reobserved in the TESS’s first extended mission (EM) in sectors 41, 48, and 49 (ut 2021 July 24 to 2022 March 25). The aforementioned discovery papers focused their TESS analyses on the 2-min Presearch Data Conditioning Simple Aperture Photometry (pdcsap) light curves (PDCSAP; Twicken et al. 2010; Morris et al. 2020), which were processed by the NASA Ames Science Processing Operations Center (SPOC; Jenkins et al. 2016). The full PDCSAP light curve is depicted in Fig. 1 and the phase-folded transit light curves from the PM and EM are included in Fig. 2.

Figure 1.

TESS PDCSAP light curves of TOI-1266 from the primary mission (PM) and first extended mission (EM). The mean GP detrending model for each set of consecutive sectors is shown in pink. The bold measurements and vertical ticks highlight the transit times of the planets TOI-1266 b (yellow) and c (blue), respectively.

Open in new tabDownload slide

Figure 2.

Phase-folded transits and transit model fits of both TOI-1266 planets from three TESS data reductions and two ground-based facilities. The transits of TOI-1266 b and c are depicted in the top and bottom rows, respectively. Different TESS extractions from the primary (PM) and extended missions (EM) include PDCSAP and our custom extractions from the TESS TPFs and FFIs. We also include ground-based, diffuser-assisted transit observations from the APO/ARCTIC (TOI-1266 c only) and Palomar/WIRC cameras, which we use to validate the variations in both planets’ transit shapes revealed by TESS between the PM and EM. The central vertical bar separates observations taken during the PM on the left and post-PM on the right (i.e. on and after ut 2021 July 24).

Open in new tabDownload slide

In this paper, we use the PDCSAP light curves in addition to custom light-curve extractions from the 2-min TESS Target Pixel Files (TPFs) and the Full Frame Images (FFIs). The latter has cadences of 30 and 10 min in the TESS PM (i.e. sectors 14, 15, 21, 22) and EM (i.e. sectors 41, 48, 49), respectively. We consider multiple light curve extractions to assess the validity of each planet’s transit shape across all seven TESS sectors (see Section 4.5). We performed a custom extraction of the TESS 2-min light curve from the TPFs using the lightkurve package (Lightkurve Collaboration 2018; Ginsburg et al. 2019; Astropy Collaboration 2022). Similarly to the PDCSAP light-curve extraction, we used the cotrending basis vectors (CBVs) for the appropriate TESS camera and CCD in each sector to correct for systematics. This allows us to have more control over the light-curve extraction to assess the sensitivity of the planetary transit depths on the extraction methodology (see Section 4.5). We tested a variety of photometric apertures including the optimized TESS aperture size and threshold apertures that include pixels within 3, 5, and 10σ of the median target flux. We also corrected for dilution by neighbouring sources based on their known positions and TESS magnitudes in the TESS Input Catalog (TIC; Stassun et al. 2019). We derived dilution values that are consistent with the values reported by the SPOC and used in the PDCSAP analysis. The phase-folded TPF transit light curves for both planets from the PM and EM are included in Fig. 2.

We performed a second custom light-curve extraction from the TESS FFIs following Vanderburg et al. (2019). We retrieved the TPFs using the TESSCut interface integrated into MAST (Brasseur et al. 2019). We extracted each sector’s light curve using 20 photometric apertures that were either circularly or PSF-shaped. In each sector, we selected the aperture size and shape that maximized the photometric precision after correcting for dilution and detrending via linear regression using the spacecraft’s quaternion time series and the PDCSAP CBVs, after binning to the appropriate FFI cadence (i.e. 30 min in the PM and 10 min in the EM). An additional spline component was included to correct for any residual photometric variability that we attribute to the star. The phase-folded FFI transit light curves for both planets from the PM and EM are included in Fig. 2.

A preliminary inspection of transit light curves in Fig. 2 revealed significant transit depth discrepancies between the PM and EM. We postpone our investigation of these anomalies until Section 4.5.

2.2 APO/ARCTIC diffuser-assisted photometry

Fig. 2 includes an archival ground-based transit of TOI-1266 c taken on ut 2020 January 28 with the diffuser-assisted ARCTIC camera on the 3.5-m ARC telescope at Apache Point Observatory (APO). These observations were originally presented by S20, with the details described therein.

2.3 Palomar/WIRC diffuser-assisted photometry

We observed one full transit of each of TOI-1226 b and c in the J band with the Wide-field Infared Camera (WIRC) on the Hale Telescope at Palomar Observatory, California, USA. Palomar/WIRC is a 5.08-m telescope equipped with a 2048 × 2048 Rockwell Hawaii-II NIR detector, providing a field of view of 8.7 arcmin × 8.7 arcmin with a plate scale of 0.25 arcsec per pixel (Wilson et al. 2003). Our data were taken with a beam-shaping diffuser that increased our observing efficiency and improved the photometric precision and guiding stability (Stefánsson et al. 2017; Vissapragada et al. 2020).

We observed the transit of TOI-1266 c on ut 2022 March 10 between airmasses of 1.65 and 1.19. Our observing window provided close to three full transit durations of coverage, centered on the transit mid-point. We collected 660 science images with total exposure times of 10.5 seconds per image. We observed the transit of TOI-1266 b on ut 2022 June 6 between airmasses of 1.26–2.26. Our observing window provided pre-ingress coverage and full transit coverage, but no post-egress baseline.

We collected 517 science images with the same total exposure time as our first transit. For both observing sequences, we initiated the observations with a nine-point dither near the target to construct a background frame that we used to correct for the effects of a known detector systematic at short exposure times. We dark-subtracted, flat-fielded, and corrected for bad pixels following the methodology of Vissapragada et al. (2020). We scaled and subtracted the aforementioned background frame from each image to remove the full frame background structure. We then performed circular aperture photometry with the photutils package (Bradley et al. 2020) on our target, plus three comparison stars. We tested apertures sizes between 7 and 25 pixels and selected the aperture that minimized the root-mean-square deviation of the final photometry. We also used uncontaminated annuli with inner and outer radii of 25 and 35 pixels, respective, for local background subtraction. We fit the light curves using the exoplanet package (Foreman-Mackey, Barentsen & Barclay 2019; Foreman-Mackey et al. 2021), which included a systematics correction constructed as a linear combination of comparison star light curves, the centroid offset, PSF width, airmass, and local background (Greklek-McKeon et al. 2023). Typically, airmass and sky background are not included simultaneously as detrending parameters, but their inclusion here provided improved fits as our observations showed considerable variations in background flux that did not track with airmass, despite the clear weather on both nights. Our light curves are included in the right most column of Fig. 2.

2.4 Precise radial velocity measurements from HARPS-N

We obtained 145 spectra of TOI-1266 with the HARPS-N spectrograph between ut 2020 May 8 and 2022 July 14. HARPS-N is a high resolution (R = 115, 000) optical échelle spectrograph located at the 3.6m Telescopio Nazionale Galileo (TNG) on La Palma, Canary Islands. HARPS-N boasts exceptional long-term pressure and temperature stability, which enables the instrument to reach sub-m/s stability on bright FGK stars (Cosentino et al. 2012). Despite its operation at optical wavelengths (i.e. 380–691 nm), HARPS-N has consistently demonstrated its ability to reach the m s⁻¹ level of precision on bright (V ≲ 13.5) early-to-mid M dwarfs (LTT 3780; Cloutier et al. 2020a, TOI-1235; Cloutier et al. 2020b, TOI-1634; Cloutier et al. 2021, TOI-1695; Cherubim et al. 2023, HD 79211; DiTomasso et al. 2023, LHS 1903; Wilson et al. in preparation).

Each HARPS-N spectrum was taken as part of the HARPS-N Guaranteed Time Observations (GTO) program. We positioned the non-science fibre on-sky in all of our observations. We fixed all exposure times to 1800 s. We reduced our spectra using the latest version of the HARPS-N Data Reduction Software v2.3.5 (DRS; Dumusque et al. 2021). We obtain a nightly RV measurement from the DRS via the cross-correlation function (CCF) technique using an M2 CCF mask. The DRS also produces time series of the full width at half-maximum (FWHM) and bisector inverse slope (BIS) activity indicators. We obtained a median peak signal-to-noise ratio (S/N) across our spectra of 33.

We performed our RV extraction using the novel line-by-line (LBL) method proposed by Dumusque (2018) and implemented by Artigau et al. (2022). We used version 0.61.0 of the lbl code¹, which performs a simple telluric correction by fitting a TAPAS model (Bertaux et al. 2014). This routine has been demonstrated to significantly improve the RV precision of M dwarfs observed with ESPRESSO (Allart et al. 2022).

The LBL method is conceptually similar to the commonly used template-matching method, which makes LBL similarly well suited to RV extraction from M dwarf spectra (Anglada-Escudé & Butler 2012; Astudillo-Defru et al. 2015). It is a data-driven method wherein the Doppler shifts of many individual spectral lines are calculated with respect to a high S/N template spectrum that we construct by coadding all our HARPS-N spectra. LBL extractions have demonstrated similar performances to template-matching in terms of RV precision and dispersion on M dwarfs observed with optical RV spectrographs (e.g. Artigau et al. 2022). We elect to use an LBL extraction over template-matching because the former is robust to outlying RV measurements that can bias a spectrum’s RV measurement well beyond the typical photon noise uncertainty. Spurious lines can arise from a variety of sources including imperfect telluric contamination, detector defects, and variable line sensitivity to magnetic activity (e.g. Lafarga et al. 2023). Spurious and high S/N signals offset a spectrum’s RV measurement when left unaccounted for, which is the case for template-matching extractions that consider all lines.

By calculating the RV shift of individual lines, we easily identify high S/N outliers and reject them, therefore producing a more robust RV measurement. In each spectrum, we reject lines that are >3σ discrepant from the median RV over all lines. An illustrative example of this is shown in Fig. 3 for a randomly selected HARPS-N spectrum. We perform our outlier rejection on the full HARPS-N wavelength domain, as well as in the BVR bands, which we will ultimately use to assess the chromaticity of RV signals (Section 4.1). The rms values and median uncertainties in each RV time series, including the CCF RVs, are reported in Table 1. Indeed, LBL offers a significant improvement over the CCF technique for TOI-1266.

Figure 3.

An illustrative example of our LBL RV extraction from a representative HARPS-N spectrum. Left panel: The measured RV and uncertainty of 22 577 individual lines (light blue) as a function of wavelength. The dark purple markers depict 48 lines that are rejected as outliers with ΔRV/σ > 3 in this spectrum, where ΔRV represents the RV differences of individual lines from the median value, and σ are the RV uncertainties on individual lines. Right panel: The distribution of normalized LBL RVs compared to the unit normal distribution (dashed purple line).

Open in new tabDownload slide

We also calculated the differential line width activity indicator dLW using the LBL framework (Zechmeister et al. 2018; Artigau et al. 2022). The dLW time series has units of m² s⁻² and scales with the FWHM of the line profile. But similarly to our LBL RVs, the dLWs are constructed in an outlier-resistant way. Because of this, and because of the generally poor performance of M dwarf binary masks in the CCF method, we consider the dLW time series to be a more robust activity indicator than the CCF FWHM and BIS. We provide the dLW, along with all of our spectroscopic LBL and CCF time series, in Table 2.

3 TOI-1266 STELLAR PARAMETERS

TOI-1266 has been studied extensively in previous efforts to validate the planetary system (D20, S20). Spectroscopic reconnaissance, high-resolution imaging, and Gaia astrometry (RUWE = 1.227; Gaia Collaboration 2023) rule out comoving stellar companions and indicate that TOI-1266 is a single star located at a distance of 36 pc.

TOI-1266 is considered magnetically inactive and does not show signs of photometric variability in its TESS light curve (Fig. 1). We are unable to recover the stellar rotation period in the TESS data. We also queried the ground-based photometric monitoring campaigns from the All-Sky Automated Survey for SuperNovae (ASAS-SN; Shappee et al. 2014; Kochanek et al. 2017) and the Zwicky Transient Facility (ZTF; Bellm et al. 2019). We find no signatures of stellar rotation in these photometric light curves but we do find evidence for P_rot ∼45 d in our spectroscopic time series, presumably from bright plages (see Section 4.1). The lack of rotational modulation in the TESS, ASAS-SN, and ZTF light curves indicates that TOI-1266’s filling factor by dark starspots is consistent with zero. We also measure |$\log {R^{\prime }_{\text{HK}}}$||$=-5.50^{+0.35}_{-0.44}$| from the Mount Wilson S-index calculated by the HARPS-N DRS (Astudillo-Defru et al. 2017), which supports the notion that TOI-1266 is chromospherically inactive.

We adopt physical stellar parameters from the TESS Input Catalog v8.2 (TIC; Stassun et al. 2019), which are consistent with the values reported in D20 and S20. Values of the stellar effective temperature, radius, and mass are based on the G_BP, G_RP, and 2MASS K_s-band magnitudes and use the empirical relations for M dwarfs from Mann, Gaidos & Ansdell (2013) and Mann et al. (2015, 2019), respectively. The TIC does not report a metallicity measurement but we note that both discovery papers present multiple [Fe/H] measurements using spectroscopic and SED-fitting techniques. The five literature values of [Fe/H] exhibit a large rms of 0.17 dex and we do not attempt to evaluate the accuracy of these measurements. Instead, we conclude that TOI-1266 appears to be somewhat metal-poor and report a range of metallicity values in Table 3, along with all other stellar parameters.

4 DATA ANALYSES

4.1 Generalized Lomb–Scargle periodograms

Fig. 4 depicts the generalized Lomb–Scargle periodograms (GLS; Zechmeister & Kürster 2009) of our spectroscopic window function, the full RVs, the chromatic BVR RVs, and the following activity indicators: H α, dLW, FWHM, and BIS. The corresponding false alarm probability (FAP) of each time series is included in Fig. 4.

Figure 4.

GLS periodograms of our HARPS-N window function, the full LBL RVs, the chromatic LBL RVs in the BVR bands, and the spectroscopic activity indicators H α, dLW, FWHM, and BIS. Left column: the GLS periodogram of the time series labelled on the y-axis. The vertical dashed lines highlight the known orbital periods of the transiting planets TOI-1266 b and c, along with the prominent peaks at 32.3 and 44.6 d in the full RVs and H α time series that we attribute to a new planetary candidate and to stellar rotation, respectively. The planetary and candidate orbital periods are only depicted in the four RV panels for clarity. The horizontal dotted lines depict each periodogram’s 1 per cent FAP level. Right column: The FAP as a function of GLS power corresponding to each periodogram.

Open in new tabDownload slide

Aside from the aliased signals close to 1 d (Dawson & Fabrycky 2010), the GLS of the full RV time series exhibits one periodicity with FAP <1 per cent at 10.9 d, which corresponds to the transiting planet TOI-1266 b. Other noteworthy peaks are seen at the orbital period of TOI-1266 c (i.e. 18.8 d; FAP = 95 per cent), plus a peak at 32.3 d that is unrelated to either known transiting planet (FAP = 4 per cent). We investigated the origin of the latter by assessing its chromatic dependence in the BVR bands and its presence in the GLSs of various activity indicators. The chromatic RVs are not informative because the B-band RVs are too low S/N, and while the 32.3-d signal does exist in the V and R-band RVs with varying FAP values, so do the known signals from the transiting planets. In particular, the prominent signal from TOI-1266 b decreases from the V-band to the R-band RVs, which is not expected as the relative power of achromatic planetary signals are expected to increase with increasing wavelength as the impact of magnetic activity is lessened (Reiners et al. 2010). Establishing the chromaticity of GLS signals is clearly questionable from our time series. We therefore refrain from commenting on the origin of 32.3-d signal until Section 4.3.

Turning to the activity indicators, we do not see a low FAP peak close to 32.3 d in any of our activity time series. The dLW, FWHM, and BIS time series do not show any signs of stellar activity and are all consistent with noise. The only prominent activity peak with FAP <1 per cent is exhibited by our H α time series at 44.6 d. H α is a strong indicator of M dwarf chromospheric activity at optical wavelengths (Lafarga et al. 2021) and we interpret the 44.6-d signal as the star’s rotation period, or a harmonic thereof. The lack of a low FAP peak at 32.3 d, or 32.3/2 ≈ 16.2 d, in any of the activity indicators leads to the possibility that this periodicity may be due to a third planet around TOI-1266. We reserve the justification of this claim until Section 4.3 wherein we compare the model evidences for two planets versus three planets based on our RV data.

4.2 Global RV + transit modelling

Here, we jointly model our HARPS-N RVs and the TESS PDCSAP light curves. We consider the TESS PM and EM data independently in order to separate the distinct transit depths between each data set.

We set the planetary components of our global model to be pseudo-Keplerian orbits. We consider pseudo-Keplerians because the period ratio of TOI-1266 b and c is within 3.5 per cent of the 5:3 mean motion resonance (MMR) and D20 showed marginal evidence for transit timing variations (TTVs) of both planets in the TESS PM. We note however that the independent analysis of the same TESS data by S20 concluded that the system does not exhibit measurable TTVs. Here, we assume the general case in which TTVs may be present and measurable.

The maximum TTV amplitudes reported by D20 are 4 and 16 min for planets b and c, respectively, but these measurements represent only marginal TTV detections at <3σ. The corresponding TTV masses are poorly constrained: |$M_{p,b}=13.5^{+11.0}_{-9.0}\, \mathrm{M}_{\oplus }$| and |$M_{p,c}=2.2^{+2.0}_{-1.5}\, \mathrm{M}_{\oplus }$| (D20). Using the WHFast integrator in the rebound N-body code (Rein & Liu 2012; Rein & Tamayo 2015), we estimate that non-Keplerian effects result in RV excursions from Keplerian motions at the level of ∼30 cm s⁻¹. These excursions are ∼2−3 × smaller than our forthcoming precision on the planets’ semiamplitudes (i.e. ∼70−90 cm s⁻¹) and do not correlate with our uniformly sampled RV observations such that they are a negligible effect. The planet-induced transit and RV signals that we aim to extract from our data are well approximated as being Keplerian, although our pseudo-Keplerian model retains the flexibility to measure TTVs. This is accomplished by fixing each planets’ orbital elements while allowing individual mid-transit times to vary. This makes us sensitive to TTVs but not to transit duration variations.

The light-curve component of our global model features the following parameters: stellar mass and radius, flux baseline f₀, scalar jitter s, quadratic limb-darkening parameters {u₁, u₂}, orbital periods P, planet–star radius ratios R_p/R_⋆, impact parameters b, and rescaled orbital eccentricities e, and arguments of periastron ω (i.e. |$\lbrace \sqrt{e}\cos {\omega },\sqrt{e}\sin {\omega } \rbrace$|⁠). The RV component includes the additional parameters: the systemic RV γ₀, scalar jitter s_RV, and each planet’s RV semiamplitude K. Our adopted priors are summarized in Table 5.

We consider two sets of models containing two and three planets, respectively. The latter set of models include a third Keplerian component to model the periodic signal at 32.3 d of insofar unknown origin (c.f. Fig. 4). Additionally, we consider separate models that do and do not include a Gaussian process (GP) regression model of temporal correlations from stellar activity in our RV time series. We elect to train the GP on our H α time series as it exhibits the periodic signal at ∼44.6 d that we attribute to stellar rotation, and because GP training on photometry is not informative when photometric signatures of activity are undetected (Tran et al. 2023). The simultaneity of the H α and RV time series also has the desirable effect of making us insensitive to temporal changes in the RV activity signal.

We use the exoplanet software package (Foreman-Mackey, Barentsen & Barclay 2019; Foreman-Mackey et al. 2021) to build our RV + transit model, which uses the celerite2 package (Foreman-Mackey et al. 2017; Foreman-Mackey 2018) to construct the GP and evaluate the likelihood of our data set under the GP model. The starry package is used within exoplanet to construct our model light curve (Luger et al. 2019). For the RV component of our global model, we adopt a covariance kernel whose power spectral density is the superposition of two damped simple harmonic oscillators (SHOs). The SHO’s damping factor produces a covariance kernel that is qualitatively similar to the commonly used quasi-periodic kernel. Similarly, we adopt a GP model of residual systematics in the TESS light curves using a kernel with one damped SHO, whose hyperparameter values are distinct from the RV component. We refrain from describing the details of our covariance kernels, as they have been described extensively throughout the literature (e.g. Cherubim et al. 2023) and in the celerite2 documentation.² Our GP hyperparameter priors are included in Table 5.

We use the pymc3 Markov chain Monte Carlo (MCMC) package (Salvatier, Wieckiâ & Fonnesbeck 2016) to sample the joint posterior probability density function (PDF) of our global model parameters. The MCMC is initialized with two simultaneous chains, each with 10³ burn-in steps and 10⁴ posterior draws. We report maximum a-posteriori (MAP) point estimates from each parameter’s marginalized posterior PDF in Table 5. Parameter uncertainties are derived from each marginalized posterior’s 16th and 84th percentiles unless noted otherwise.

The detrended transit light curves were shown in Fig. 2, while the individual RV components are in Fig. 5 for our fiducial model containing three Keplerian components and a GP. The GLS periodogram of each Keplerian component reveals that the Keplerian’s period is the lowest FAP signal, while the stellar rotation period of 44.6 d is not present in the GLS of the RV activity component. However, in Section 4.3 we will show that our fiducial model, which includes a GP, is favoured. The GLS of the RV residuals shows no evidence for significant residual periodicities. Fig. 6 depicts the phase-folded RV time series for each planet.

Figure 5.

Individual components of our 3-planet RV model of TOI-1266 and corresponding GLS periodograms. From top to bottom: The raw LBL RVs, stellar RV activity, TOI-1266 b, TOI-1266 c, the planet candidate TOI-1266 d, and the residuals. The residuals do not show any evidence for unaccounted periodicities in the GLS. The vertical shaded regions in the left column highlight each TESS sector’s observing window (annotated at the top) and reveals, where we have contemporaneous TESS and spectroscopic observations. Similarly, the vertical dashed lines highlight the epochs of our ground-based transit observations of either TOI-1266 b or c, using the ARCTIC and WIRC cameras. The vertical dashed lines in the right column depict the orbital periods of the transiting TOI-1266 planets, the candidate planet d, and the stellar rotation period inferred from our H α time series.

Open in new tabDownload slide

Figure 6.

Phase-folded RV time series for the transiting planets TOI-1266 b (left) and c (middle), and the planet candidate TOI-1266 d (right).

Open in new tabDownload slide

4.3 Evidence for a third planet around TOI-1266

In Section 4.2, we considered four RV + transit models with either two or three planets, which we denote as 2p and 3p, respectively. The third planetary component in our 3p models is used to model the 32.3-d signal seen in the GLS of our RV time series. We also considered two versions of our 2p and 3p models that do and do not include a stellar activity component in the form of a trained GP. For each model M, we estimate its Bayesian model evidence |$\mathcal {Z}_M$| using the Perrakis estimator, for which we use the model’s joint posterior from our MCMC calculations as an importance sampler (Perrakis, Ntzoufras & Tsionas 2014). We then compute Bayes factors, or evidence ratios, to determine the favoured model given our RV data. Our results are reported in Table 4.

We find that in both our 2p and 3p models, the inclusion of a GP is favoured.³ Despite TOI-1266 being an inactive star, the inclusion of a GP activity component produces a more complete model. We also find that our 3p models are heavily favoured over our 2p models, regardless of whether a GP is included. Specifically, |$\mathcal {Z}_{3p+\mathrm{GP}}/\mathcal {Z}_{2p+\mathrm{GP}} = 10^{6.4}$| and |$\mathcal {Z}_{3p}/\mathcal {Z}_{2p} = 10^{5.6}$|⁠. This suggests that the 32.3-d signal in the GLS of the RVs is well described by a Keplerian orbital solution. We consider this to be tentative evidence for a third planet candidate in the system whose planetary parameters are included in Table 5.

4.4 TTV analysis

We do not find strong evidence for TTVs in the TESS photometry. This result is wholly consistent with our measured RV masses given the relatively large photometric rms of the individual TESS transit events, which inhibit the precise measurement of individual transit times. This result is consistent with previous findings based on TESS PM data (S20). However, with the planets’ masses constrained by our RV measurements, here we supplement our pseudo-Keplerian analysis by modelling the lack of TESS TTVs to place stringent upper limits on the transiting planets’ eccentricities following the formalism of Hadden et al. (2019).

Since the transit times of TOI-1266 c may be influenced by dynamical interactions with the candidate planet d, we focus on deriving eccentricity constraints by modelling TOI-1266 b’s transit times only. We do so by adopting a linear model approach as described in Linial, Gilbaum & Sari (2018). We are sensitive to the magnitude of ‘combined eccentricity’ in this model, |$|\mathcal {Z}_{{\rm b},{\rm c}}|\approx 1/\sqrt{2}\sqrt{e_{\rm b}^2 + e_{\rm c}^2-2e_{\rm b}e_{\rm c}\cos (\varpi _{\rm b}-\varpi _{\rm c})}$|⁠, which depends on the planet eccentricities e_i and arguments of periastron ϖ_i. We use the emcee package (Foreman-Mackey et al. 2013) to sample the likelihood of TOI-1266 b’s transit times under this model, which is conditioned on the RV masses reported in Table 5. We conclude that |$|\mathcal {Z}_{{\rm b},{\rm c}}| \lt 0.08$| (0.11) at 95 per cent (99 per cent) confidence, which is more stringent than our eccentricity constraints derived from the RVs alone.

4.5 Transit depth discrepancy

As introduced in Section 2.1 and illustrated in Fig. 2, our initial inspection of the PDCSAP transit light curves revealed transit depth discrepancies for both transiting planets between TESS’s PM and EM. Phase-folding the PM and EM light curves based on the linear ephemerides from our preliminary analysis revealed that TOI-1266 b’s transit depth increases by 1.9σ from |$2.86^{+0.12}_{-0.14}$| ppt to |$3.09^{+0.14}_{-0.14}$| ppt. More significantly, the transit depth of TOI-1266 c increases by 3.5σ from |$1.45^{+0.18}_{-0.17}$| ppt to |$2.05^{+0.18}_{-0.17}$| ppt. We checked that these distinct transit depths persist across multiple transit events in the PM and EM and are not dominated by any single event that may result from a transient event (e.g. a solar system object occultation).

We seek to resolve these discrepancies by considering five possible explanations: variable flux dilution, orbital precession, residual artefacts in the PDCSAP light curve, evolution in stellar activity, and stochastic effects between the PM and EM TESS data. In the following subsections, we rule out the first three scenarios and ultimately attribute the observed transit depth discrepancy to a marginal increase in stellar activity during the TESS EM and to stochastic effects in the TESS PM data that limit the accuracy of recovered low S/N signals with the ephemeris of TOI-1266 c.

4.5.1 Variable flux dilution?

As discussed in Section 2.1, we accounted for dilution in our TPF and FFI-extracted light curves using the known positions and TESS magnitudes of neighbouring sources from the TICv8 (Stassun et al. 2019). Yet variable dilution between the PM and EM can be immediately ruled out as the culprit for the transit depth discrepancies because its impact on both planets would be equal. The transit depth ratios between the PM and EM are 0.88 ± 0.02 and 0.71 ± 0.04 for TOI-1266 b and c, respectively. From these differing levels of depth increases, we conclude that variable flux dilution is not responsible for the transit depth discrepancies.

4.5.2 Orbital precession?

Gravitational interactions between planets in compact systems can cause precession of the planets’ orbital planes, which can produce an observable signature by altering a transiting planet’s impact parameter (Miralda-Escudé 2002). Changing the impact parameter can result in the occultation of differentially limb-darkened transit chords, which may be able to explain the transit depth discrepancy. This effect has been observed in the K2-146 system containing two sub-Neptunes in a 3:2 resonance (Hamann et al. 2019). However, in comparison, the TOI-1266 planets are near a higher order period ratio (i.e. 5:3) that produces weaker gravitational interactions such that we do not expect to see significant orbital precession. We test this by revisiting our transit light curve fits from the PM and EM (see Section 4.2) and noting that each planet’s transit duration are consistent within 1σ between the PM and EM. We conclude that the TOI-1266 system does not show evidence for orbital precession over the 16-month interval between the PM and the EM such that orbital precession cannot explain the observed transit depth discrepancy.

4.5.3 Residual artefacts in the PDCSAP light curve?

We investigate whether the transit depth discrepancies are unique to the SPOC’s PDCSAP light-curve construction by considering alternative TESS data reductions and precise transit light curves that were not obtained with TESS. As outlined in Section 2.1, we used our custom light curve extractions from each TESS sector based on the TPFs and FFIs. We corrected each light curve for dilution and detrended using a GP following our methodology in Section 4.2.

Fig. 7 and Table 6 report each planet’s transit depth and duration measured by TESS, or from the ground, in the TESS PM, EM, and post-EM. Across the three TESS extractions considered (i.e. PDCSAP, TPF, and FFI), both planet’s transit depths increase from the PM to the EM epochs while their transit durations remain constant. These findings are consistent across the different TESS extractions.

To verify that these depth discrepancies are not unique to TESS, we sought to obtain high-precision transit observations from the ground using diffuser-assisted imagers. While no such observations are known to exist for TOI-1266 b around the PM epochs, we retrieved the full transit observation of TOI-1266 c from S20 taken with the ARCTIC camera (Section 2.2). These data confirm the PM transit depth and duration of TOI-1266 c observed by TESS (Fig. 7). We followed-up with additional transit observations of TOI-1266 b and c in the J band from the WIRC camera (Section 2.3). These ground-based data are consistent with the deeper transits seen by TESS in the EM versus the PM, and confirm the invariability of both planets’ transit durations over time. We conclude that the observed depth discrepancies are not unique to TESS and that they are robust to the data reduction steps and across different observing facilities.

Figure 7.

Observed transit depths and durations for TOI-1266 b (yellow) and c (blue) as a function of time and based on different observing facilities or different TESS data reductions (different markers). The horizontal bands depict the median and standard deviations of each planet’s depth and duration values across all facilities but in the PM and EM epochs separately. One exception is the transit duration of TOI-1266 c from the TESS PM PDCSAP light curve, which we omit due to its anomalous nature when compared to the other extraction methods and facilities.

Open in new tabDownload slide

4.5.4 Evolution in stellar activity?

We investigate whether active region (AR) evolution could plausibly explain the observed transit depth discrepancies. The TESS photometric rms does not vary significantly between the PM and EM (Table 6), such that we can rule out dark starspots as the culprit. Here, we consider bright plages by simulating their impact on the TESS photometry, HARPS-N LBL RVs (less the known planets), CCF FWHM, and CCF BIS time series using an adapted version of soap 2.0 (Dumusque, Boisse & Santos 2014; Wurmser et al. in preparation). We explored a grid of N = 43 890 plage sizes (⁠|$0{\!-\!}0.65\, R_\star$|⁠), temperature contrasts (0−1000 K), latitudes (⁠|$0{\!-\!}90\deg$|⁠), and longitudes (⁠|$0{\!-\!}360\deg$|⁠), and determined the activity-induced rms values on each of the aforementioned time series. Planet population studies suggest that small planets are preferentially well-aligned with their host stars’ rotation axes (Albrecht, Dawson & Winn 2022) such that we assume TOI-1266’s rotation axis is perpendicular to our line of sight. We adopt P_rot = 30 d as these calculations were performed before the identification of the probable rotation period at 44.6 d. Fortunately, our preliminary tests sampled orbital periods from 30 to 60 d and revealed that our results are insensitive to the exact choice of P_rot.

We compare our grid of plage models to our data in each of our four observing seasons of spectroscopic follow-up (see Fig. 5). In this way, we are sensitive to changes in the plage properties that may produce the observed transit depth discrepancies. In each model, we consider plages to be observationally ‘allowed’ if each of the model rms values in photometry, RVs, FWHMs, and BISs are less than the observed rms of the corresponding time series. Fig. 8 shows the results of this exercise. We find that small, cool plages at high latitudes (R_plage/R_⋆ ≲ 0.17, ΔT ≲ 300 K, ϕ ≳ 80°) are consistent with the data in all four observing seasons. However, a larger fraction of models is consistent with the data in season 3 wherein plages, of a fixed size, extend down to lower latitudes in season 3 (ut 2021 December 18 to 2022 July 14). We conclude that season 3 exhibits the highest level of activity out of the four seasons. Season 3 corresponds to epochs that overlap with sectors 48 and 49 of the TESS EM wherein we witness the transit depth increases compared to the PM.

Figure 8.

The results of our active region simulations using soap 2.0. The panels depict the fraction of plage configurations that are consistent with our observations in each of our spectroscopic follow-up observing seasons and as a function of latitudinal position and temperature contrast with the photosphere (i.e. marginalized over longitude and plage size). Cool plages at high latitudes are allowed in each observing season with season 3 exhibiting the highest number of allowed models, indicative of a slight increase in TOI-1266’s magnetic activity level.

Open in new tabDownload slide

Our finding that plages on TOI-1266 are localized at high latitudes is consistent with the findings of Morales et al. (2010) who analysed a set of eclipsing binaries and concluded that polar magnetic ARs on low-mass stars are a preferred solution to explain differences between stellar model predictions and EB measurements of fundamental stellar parameters. The existence of polar plages on TOI-1266 also provides a plausible explanation for the more severe transit depth discrepancy exhibited by TOI-1266 c than b due to its larger impact parameter (i.e. b_c = 0.77; b_b = 0.55).

Because TOI-1266’s magnetic activity is dominated by bright plages rather than by dark spots, the impact of ARs on the observed transit depths will be to increase the transit depth relative to an unspotted photosphere. Morris et al. (2018) suggest a framework to disentangle the effects of R_p/R_⋆ and unocculted ARs to recover robust planetary radii from transit depth measurements. The basis of the framework is that while the true value of R_p/R_⋆ affects both the transit depth and the ingress/egress durations, only the depth is impacted by unocculted ARs. The ingress/egress durations may therefore be used to infer accurate R_p/R_⋆ values if the transit light curve has sufficient S/N and limb-darkening is well-constrained. In light of the possibility of elevated plage coverage in observing season 3, we attempted to use this framework to measure each planet’s R_p/R_⋆ value in the full TESS EM, the individual EM sectors, and in the diffuser-assisted WIRC transits. We find that in all light curves, the photometric S/N values are too low to precisely recover the ingress/egress times such that we measure planetary radii that are consistent between spotted and unspotted photospheres. That is, while there exists some evidence that TOI-1266’s activity level increases slightly towards the TESS EM epochs, this does not have a demonstrable impact on the observed transit depths.

4.5.5 Stochastic effects between the PM and EM TESS data

We consider whether there are stochastic effects in either the PDCSAP PM or EM data that drive the disagreement in transit depths between their respective light curves. Here, we focus on TOI-1266 c, which shows a greater depth discrepancy than b. We proceed by injecting transit signals into our cleaned TESS data, recovering them within our model fitting framework (Section 4.2), and comparing the recovered transit depths to those injected. We sample transit depths from a log uniform grid spanning Z_inj ∈ [0.5, 3] ppt. We separately consider two sets of ephemerides for our injected planets: (i) randomized mid-transit times and orbital periods between 15–25 d and (ii) injected planets with the same ephemeris as TOI-1266 c. We inject these signals into the cleaned PDCSAP light curves after removing our MAP GP and transit models from the PM and EM, respectively.

Fig. 9 compares the resulting distributions of transit depth differences, ΔZ = Z_rec − Z_inj, for the PM and EM and for each set of injected ephemerides. We treat ΔZ as a measure of the accuracy with which transit depths are recovered. We find that when planets are injected with random ephemerides, the ΔZ distributions are Gaussian-distributed in both the PM and EM, and with similar rms values of 0.21 and 0.25 ppt, respectively. The values are comparable to the photometric rms in the PM and EM light curves (Table 6). This implies that transiting planets with random transit times and orbital periods can be accurately recovered in both the PM and EM PDCSAP light curve. Conversely, we find that injected planets with the same ephemeris as TOI-1266 c exhibit depth measurement accuracies that differ between the PM and EM. This effect is not due to the smaller number of transits observed because our randomized injected ephermides have comparable periods to TOI-1266 c. We note that the depth of a c-like planet is more accurately recovered in the EM than in the PM. This is evidenced by the distribution of ΔZ in the EM being narrowly distributed around 0.03 ppt, whereas the typical accuracy of transit depths measured in the PM is considerably larger at 0.14 ppt.

Figure 9.

Cumulative distributions of transit depth measurement accuracies ΔZ from a set of injection/recovery tests. Four sets of models are shown here comparing results from the TESS PM and EM, and for two different sets of injected planet ephemerides: (i) with random orbital periods and mid-transit times and (ii) with the linear ephemeris of TOI-1266 c.

Open in new tabDownload slide

Recall that the observed transit depth discrepancy for TOI-1266 c is 0.6 ppt, which is greater than the level of inaccuracy that our injection/recovery tests suggest in the PM. However, our results do indicate that the measured transit depth of planet c is more reliable in the EM than in the PM. Our results also indicate that this is a stochastic effect related to the exact ephemeris of TOI-1266 c because signals injected with random ephemerides do not show the same systematic overestimate in transit depth. This may explain why the transit depth discrepancy is significant for planet c (3.5σ) while it is marginal for b (1.9σ).

4.5.6 Conclusions regarding the transit depth discrepancy

We conclude that variable flux dilution, orbital precession, and residual artefacts in the PDCSAP light curve cannot explain the observed transit depth discrepancies between the TESS PM and EM. While there is some evidence that the star’s magnetic activity level increased during our third observing season of spectroscopic follow-up, which aligns with the EM epochs that show increased depths, the impact is small and is likely an incomplete explanation of the apparent transit depth increase. There is also evidence that the transit depth of a planet with TOI-1266 c’s ephemeris is recovered more accurately in the EM versus the PM. Although the level of improvement is similarly insufficient to explain the observed transit depth discrepancy. We conclude that part of the transit depth discrepancy can be explained by a small increase in the stellar activity and to stochastic effects that limit the accuracy with which TOI-1266 c’s transit depth is measured in the PM. We proceed with adopting our fit results from the EM but caution that there is still room for further clarification on the transit depth of TOI-1266 c, requiring additional follow-up with ultra-precise photometers.

4.6 Transiting planet search

We conducted an independent search for transiting planets using the Transit Least Squares algorithm (tls; Hippke & Heller 2019). We searched for transit-like signals between 19 and 50 d in the combined PDCSAP light curve after detrending for residual systematics and the two known transiting planets. We chose our period bounds because we do not expect to have enhanced sensitivity to transiting planets with P < P_c = 18.8 d compared to the TESS team, and because the transit probability drops below 1 per cent at ∼50 d. TOI-1266 b and c are remarkably coplanar (Δi < 0.1°) and a third coplanar planet would no longer exhibit a transiting configuration beyond 27.1 d. We also conducted a focused TLS search between 31.9 and 32.7 d to check if the RV planet candidate revealed in Section 4.1 is transiting. We also conducted a by-eye search for this planet candidate based on its assumed linear ephemeris from our RV analysis.

Our search yielded no new transit-like signals. We proceeded with injecting synthetic planetary signals into our TESS data to quantify our sensitivity to transiting planets. We sampled orbital periods logarithmically between 19 and 50 d, uniform times of mid-transit, planet radii logarithmically between |$0.5\,{\text{and}}\,3\, \mathrm{R}_{\oplus }$|⁠, and impact parameters linearly between 0 and 1. We constructed 10⁴ synthetic light-curve realizations using the batman package (Kreidberg 2015) to compute the transit models. We then processed each light curve using the TLS and recovered planets with signal detection efficiencies >6 (Hippke & Heller 2019) within 5 per cent of the injected period.

The top panel of Fig. 10 depicts our detection sensitivity to transiting planets. We do not correct for the geometric probability given the preference of compact multiplanet systems to exhibit small mutual inclinations ≲1°−2° (Ballard & Johnson 2016). We confidently rule out transiting planets down to |$1.4\, \mathrm{R}_{\oplus }$| within 23 d but lack sensitivity to planets within the habitable zone (HZ) bounded by the water-loss and maximum greenhouse limits (36–95 d; Kopparapu et al. 2013).

Figure 10.

Our detection sensitivity to additional planets in the TOI-1266 system based on our PDCSAP TESS photometry (top panel) and HARPS-N RVs (bottom panel). We only search for transiting planets beyond the orbit of TOI-1266 c (P_c = 18.8 d) and out to 50 d where our sensitivity becomes consistent with 0 per cent. The dashed lines depict lines of constant transit S/N (upper panel) and constant RV semiamplitude (lower panel). The diamond markers highlight TOI-1266 b (yellow), c (dark blue), and the non-transiting planet candidate d (light blue). The vertical dotted line depicts the maximum orbital period that a hypothetical third planet could have and still exhibit a transiting configuration if it was coplanar with TOI-1266 b and c. The shaded region highlights the habitable zone between 36 and 95 d.

Open in new tabDownload slide

4.7 Radial velocity planet search

We conducted a similar exercise to assess our detection sensitivity to planets in our HARPS-N RV time series. We sampled the injected planetary parameters identically as in Section 4.6 with the exception that because we do not require injected planets to be transiting, we revise our impact parameter sampling to be linearly sampled from ±a/R_⋆, where a is determined by the sampled value of the orbital period. These sensitivity calculations therefore represent conservative estimates, as we are not preferentially sampling coplanar systems. Additionally, we logarithmically sample planet masses between |$0.5 \,{\rm and}\, 20\, \mathrm{M}_\oplus$|⁠. We constructed 10⁴ synthetic RV time series by computing the Keplerian orbit for each synthetic planet and injected individual planets into our HARPS-N residuals (i.e. bottom panel in Fig. 5). We deem successful recoveries of injected planets to be those that exhibit a FAP ≤0.1 per cent in the GLS within 5 per cent of the injected period and have ΔBIC <−10, where ΔBIC is the Bayesian Information Criterion comparing the injected Keplerian model to the null hypothesis (i.e. a flat line). The bottom panel of Fig. 10 depicts our detection sensitivity to planets in our RV time series as a function of orbital period and true mass. We confidently rule out planets down to |$3\, \mathrm{M}_\oplus$| within 20 d and planets down to |$4\, \mathrm{M}_{\oplus }$| within the HZ, noting however that we are only sensitive to the planet’s minimum mass if it is only detected in the RV data.

5 RESULTS AND IMPLICATIONS FOR THE SYSTEM’S FORMATION

5.1 Fundamental planetary parameters

We measure orbital parameters for the transiting planets TOI-1266 b and c that are consistent with values reported in their discovery papers (D20, S20). Conversely, our transit analysis suggests increased planet radii that arise from increased transit depths exhibited by the TESS EM data compared to the PM. In Section 4.5 we deemed the PM to suffer from stochastic effects that weaken the accuracy of TOI-1266 c’s measured transit depth. We consequently adopt the results from our EM model in the forthcoming discussion.

We recover planetary radii of |$R_{p,{\rm b}}= 2.62\pm 0.11\, \mathrm{R}_{\oplus }$| and |$R_{p,{\rm c}} = 2.13\pm 0.12\, \mathrm{R}_{\oplus }$|⁠. We also measure the masses of both planets using data from the HARPS-N spectrograph (⁠|$M_{p,{\rm b}} =4.23\pm 0.69\, \mathrm{M}_{\oplus }$| and |$M_{p,{\rm c}}=2.88\pm 0.80\, \mathrm{M}_{\oplus }$|⁠). Our results correspond to 24σ and 19σ radius measurements and 6.1σ and 3.6σ mass measurements for TOI-1266 b and c, respectively. We also uncover a third planet candidate in the system and measure its minimum planet mass |$M_{p,{\rm d}}\sin {i}=4.59^{+0.96}_{-0.94}\, \mathrm{M}_{\oplus }$| (4.8σ). If real, TOI-1266 d would be a temperate super-Earth or sub-Neptune located just inside the inner edge of the water-loss HZ at an instellation of |$F_{\rm d} = 1.3^{+0.3}_{-0.2}\, \mathrm{F}_{\oplus }$|⁠.

The left panel in Fig. 11 compares the masses and radii of TOI-1266 b and c to the current population of small planets around M dwarfs with well-characterized masses (i.e. ≥3σ). The TOI-1266 planets occupy a unique region of the mass–radius parameter space as they represent the lowest density known sub-Neptunes around M dwarfs in the size range of |$1.7{\!-\!}3\, \mathrm{R}_{\oplus }$|⁠. If we assume that both planets have Earth-like cores with core mass fractions CMF of 33 per cent, and are enveloped in H₂ gas, then at their present-day equilibrium temperatures of 425 and 354 K, respectively, we infer envelope mass fractions of X_{env, b} = 5.5 ± 0.7 per cent and X_env,c = 3.5 ± 0.8 per cent.

Figure 11.

Left panel: the mass-radius-instellation diagram for small transiting planets around M dwarfs. The marker symbols depict planets with different bulk compositions based on their masses and radii with the TOI-1266 planets, highlighted with bold markers. The solid curves illustrate interior structure models from Zeng & Sasselov 2013 with iron core mass fractions (CMF) of 0 per cent, 33 per cent (i.e. Earth-like), and 100 per cent, and water-rich models from Aguichine et al. 2021 with water mass fractions (WMF) of 10 per cent, 20 per cent, and 50 per cent at the equilibrium temperatures of TOI-1266 b (yellow) and c (blue). The shaded region highlights the densest small planets predicted by models of maximum collisional stripping (Marcus et al. 2010). Right panel: the same planet sample in period-radius-instellation space. The shaded region is bounded by model predictions of the location of the M dwarf radius valley from TDML and GDF (Cloutier & Menou 2020). Multi-planet systems are connected by the dashed lines in both panels.

Open in new tabDownload slide

However, multiple lines of recent empirical evidence have suggested that unlike around FGK-stars, sub-Neptunes around M dwarfs are inconsistent with being gas-enveloped terrestrials (Cloutier & Menou 2020; Diamond-Lowe et al. 2022; Luque & Pallé 2022; Cherubim et al. 2023; Piaulet et al. 2023). Instead, the small, close-in planets around M dwarfs that are too low density to be super-Earths are likely irradiated ocean worlds with substantial water mass fractions (WMF) and negligible envelope mass fractions. Adopting the water-rich mass–radius relations from Aguichine et al. (2021), which are state-of-the-art by their use of a multiphase and non-isothermal equation of state for water, we find that WMF_b > 80 per cent at 95 per cent confidence⁴ and WMF_c = 59 ± 14 per cent. The mass–radius relationships from Aguichine et al. (2021) are only computed down to equilibrium temperatures of 400 K, which required us to extrapolate to compute WMF_c at T_eq = 354 K.

We note that at solar composition, the WMF beyond the water snow line is expected to be at most 50 per cent (Lodders 2003) such that TOI-1266 c may be a bona-fide water-world. We highlight that while we cannot definitively classify TOI-1266 c as a water-world with solely the data herein, we present atmospheric escape models in Section 5.2 that provide additional evidence that planet c likely has a water-rich bulk composition. However, we note that breaking this degeneracy and establishing the true bulk composition of TOI-1266 c will require detailed atmospheric characterization to measure its atmospheric mean molecular weight. We report the TSM values of both planets in Table 5.

We note that if TOI-1266 c turns out to be a confirmed water-world, the system’s architecture would be one-of-a-kind as the only M dwarf system to feature a gas-enveloped sub-Neptune orbiting interior to a rocky planet or a water-world (c.f. right panel of Fig. 11). The only comparable ‘inverted’ planet pairs are TOI-270 c and d, which are both consistent with being water worlds (Van Eylen et al. 2021), and Kepler-26 b and c, which are both consistent with being gas-enveloped. While the planet formation process is a stochastic one, ‘inverted’ architectures with planets of distinct compositions are not predicted by population synthesis models as gas accretion should be enhanced at larger orbital separations due to the increased Hill radius and the cooler local gas temperatures that facilitate gas accretion (e.g. Raymond et al. 2018; Burn et al. 2021). Indeed, most Kepler multi-planet systems are not ‘inverted’ (i.e. 0.65 ± 0.4 per cent of Kepler multis; Weiss et al. 2018). Comparative studies of the atmospheric chemistry of both planets in inverted planet pairs will allow for direct tests of sub-Neptune formation models and may help shed light on their potentially unique formation pathway.

5.2 Implications of the TOI-1266 system for the emergence of the M dwarf radius valley

Recall that both TOI-1266 discovery papers reported a radius for TOI-1266 c that was suggestive of a terrestrial composition (i.e. |$\sim 1.6\, \mathrm{R}_{\oplus }$|⁠; D20, S20). If true, the TOI-1266 planets would span the radius valley and serve as an important testbed for radius valley emergence models (e.g. Owen & Campos Estrada 2020; VanWyngarden et al. in preparation). However, our mass and revised radius measurements indicate that TOI-1266 c is not consistent with a terrestrial composition such that the planet pair does not span the radius valley (c.f. Fig. 11). We therefore do not use the multitransiting TOI-1266 system to test for consistency with atmospheric escape and primordial radius valley formation models, although we do test their plausibility in the following subsections.

5.2.1 H/He gas accretion

Here we confirm that the low masses of TOI-1266 b and c are able to accrete enough H/He to explain their masses and radii. The maximum envelope mass that a core of mass M_core can accrete is the maximally cooled isothermal mass M_iso (Lee & Chiang 2015). This limit leads to the possibility of a primordial radius valley wherein low-mass cores (⁠|$\lesssim 1-2\, \mathrm{M}_{\oplus }$|⁠; Lee, Karalis & Thorngren 2022) cannot accrete massive envelopes and consequently form the population of super-Earths. Conversely, more massive cores do accrete H/He envelopes to become sub-Neptunes. This primordial radius valley model has been shown to reproduce the empirical radius valley around FGK stars (Lee, Karalis & Thorngren 2022).

We calculate M_iso for TOI-1266 b and c following the framework and assumptions outlined in Lee, Karalis & Thorngren (2022). Assuming in situ formation, we find that for gas accretion timescales ≥1 Myr, both planets are capable of accreting X_iso = M_iso/M_core greater than their observed X_env. If we relax the assumption of in situ formation, as is suggested by the system’s proximity to a low-order MMR via convergent disk migration (e.g. Cresswell & Nelson 2006; Tamayo et al. 2017), then formation in the outer planetary system only serves to increase both planet’s maximum envelope masses as accretion is enhanced at lower gas temperatures. We conclude that both planets are capable of accreting massive envelopes that are in excess of what their present-day masses and radii suggest.

Given that H/He accretion onto the TOI-1266 planet cores is plausible, here we estimate their expected initial envelope mass fractions |$X_{\mathrm{env,0}} \approx 0.02\, (M_{\mathrm{core}}/\mathrm{M}_{\oplus })^{4/5}\, (T_{\mathrm{eq}}/10^3\, \mathrm{K})^{-1/4}$| (Ginzburg, Schlichting & Sari 2016). We calculate that both planets formed with X_{env, 0} ∼ 6−8 per cent, which is plausible given that these values exceed their present-day envelope mass fractions of 5.5 per cent and 3.5 per cent, respectively.

5.2.2 Hydrodynamic atmospheric escape modelling

Given that both planets are capable of accreting substantial H/He envelopes, here we test the plausibility that those primordial envelopes can be retained. We consider thermally driven escape by XUV-driven photoevaporation and core-powered mass loss, both separately and simultaneously. We model these processes using the numerical model isofate (Cherubim & Wordsworth in preparation) wherein the XUV-driven component follows the prescription detailed in Cherubim et al. (2023), while the CPML model follows from Gupta & Schlichting (2020). Here, we provide a brief summary.

We adopt models of planetary structure following conventional assumptions of a clear H/He envelope at solar metallicity with a deep convective zone beneath an outer, isothermal radiative zone at T_eq,i. The star remains in the saturated activity regime for 500 Myr with L_XUV/L_bol = 10⁻³ before decaying as a power law, consistent with expectations from semi-empirical MUSCLES spectra (Peacock et al. 2020). Our model also includes the thermal contraction (Lopez & Fortney 2014) and Ly-α cooling following the recombination-limited escape regime (Murray-Clay, Chiang & Murray 2009). We do not explicitly model gas accretion onto the planetary cores and instead consider a grid of initial envelope mass fractions X_env,0 ∈ [0.01, 0.1].

We evolve the atmospheres under the aforementioned effects and present our results in Fig. 12. We find that neither planet is completely stripped after 5 Gyr for all but the smallest X_env,0 of 0.01. Yet despite retaining some H/He, we find that it is difficult to reconcile the observed radius of planet b with any model that includes XUV-driven photoevaporation, either solely or alongside core-powered mass loss, because its envelope is efficiently lost, leaving behind X_env,b < 5.5 per cent. Meanwhile, the observed radius of TOI-1266 c can be recovered when X_env,0,c ≳ 0.03.

Figure 12.

Results from hydrodynamic escape models that combine XUV photoevaporation, core-powered mass-loss, and thermal contraction for TOI-1266 b (top) and c (bottom). We present models for different initial envelope mass fractions X_env,0 ∈ [0.01, 0.1]. The shaded regions in each panel depict the planet’s radius posterior. We note that the apparent discontinuity beyond 10⁹ yr in our model of TOI-1266 b with X_env,0 = 0.01 occurs when complete atmospheric loss is achieved, leaving behind a bare core.

Open in new tabDownload slide

We conclude that TOI-1266 b is more susceptible to atmospheric escape compared to c. While there are instances of our escape model that are consistent with the mass and radius of TOI-1266 c, the same mass loss effect on planet b always results in too much atmospheric loss. Therefore, because a H/He layer is needed to explain the mass and radius of planet b, our atmospheric escape model strongly suggests that planet c is inconsistent with a purely H/He-enveloped terrestrial planet without a substantial volatile mass fraction. Our escape model therefore supports the hypothesis that TOI-1266 c is water rich.

5.2.3 Reconciling with expectations from hydrodynamic escape

Our hydrodynamic models are based on what may be referred to as a ‘standard picture’ of atmospheric escape of H/He envelopes atop terrestrial cores. Such models have successfully reproduced the observed radius valley around FGK stars and ruled out core compositions that are either water or iron-rich (Owen & Wu 2017; Wu 2019; Gupta & Schlichting 2020). Our results indicate that the TOI-1266 planets are inconsistent with this standard picture.

For example, our models suggest that TOI-1266 b is particularly susceptible to XUV-driven escape. This planet may therefore require an alternative explanation that prevented its H/He envelope from being completely stripped during the first ∼500 Myr. Speculative explanations include the possibility of its late arrival at its current orbit after the star’s saturated XUV phase or by the inhibition of hydrodynamic escape by either (i) atmospheric shielding by high albedo hazes or (ii) the sequestration of H₂ in the planet’s surface magma that is subsequently outgassed at late times as the core cools (Chachan & Stevenson 2018). Other amendments to the so-called standard picture are the inclusion of additional chemical components that impact atmospheric thermal structure. This includes silicate vapour, which introduces a near-surface radiative layer and alters a planet’s size compared to H/He gas alone (Misener & Schlichting 2022). It may also include a substantial volatile mass fraction if the planet forms water-rich from beyond the snow line (e.g. Venturini et al. 2020).

5.3 Lines of evidence of a unique formation pathway

5.3.1 Low-density Sub-Neptunes around M dwarfs

There is an outstanding tension between the mass–radius relations of sub-Neptune-sized planets based on RV-derived versus TTV-derived masses (Jontof-Hutter et al. 2014; Steffen 2016). Specifically, beyond orbital periods of 11 d, RV-derived planetary masses appear to be systematically larger than those measured with TTVs, across all radii |$\lt 8\, \mathrm{R}_{\oplus }$| (Mills & Mazeh 2017). This tension is exhibited at the planet population level rather than in individual systems that have high quality RV and TTV data (Malavolta et al. 2017; Borsato et al. 2019). A plausible explanation for the population level discrepancies is that TTV signal detections are increasingly likely in compact systems with small period ratios (i.e. ≲ 2), particularly those close to low-order MMRs. This introduces a bias in which TTVs preferentially probe compact planet pairs with small period ratios compared to RV measurements. If it is true that sub-Neptune pairs preferentially form with tight spacings due to gas damping (Dawson, Lee & Chiang 2016), then TOI-1266 b and c may be tracing such a formation pathway given their low period ratio of P_c/P_b = 1.72. However, the crux of this argument remains that these planets appear to have distinct compositions of gas enveloped and water rich, which is challenging to reconcile with this picture of formation.

It is worth highlighting that TOI-1266 is somewhat unique in that it is an M dwarf system that exhibits a low period ratio and has precisely measured RV masses. While the comparison of our RV-derived versus TTV-derived masses in the system is premature due to the lack of a robust TTV detection from TESS, we highlight that these planets are among the lowest density planets transiting M dwarfs between |$1.7{\,\text{and}\,}3\, \mathrm{R}_{\oplus }$| (Fig. 11). This is consistent with multi-planet systems exhibiting lower planet densities (Rodríguez Martínez et al. 2023) and suggests that the TOI-1266 planets are more consistent with the aforementioned population of systematically low-density TTV planets that may form via a distinct pathway than the bulk of other small planets around M dwarfs. High-precision ground-based measurements capable of detecting smaller TTVs than TESS can enable a more complete comparison of the RV versus TTV-derived masses in the TOI-1266 system (Greklek-McKeon et al. in preparation).

5.3.2 Impact of metallicity on planet composition

Although uncertain, TOI-1266 appears to have a sub-solar metallicity ([Fe/H] ∈ [−0.31, −0.08] dex), which puts the planets’ low bulk densities further at odds with predictions from planet formation models. This is because high metallicity environments are where we expect massive cores to form rapidly enough to trigger the accretion of a substantial gas envelope (Dawson, Chiang & Lee 2015). With X_env,b = 5.5 ± 0.7 per cent, TOI-1266 b’s envelope mass fraction is one of the highest among planets around M dwarfs and is comparable to the puffy sub-Neptune GJ 1214 b (⁠|$X_{\mathrm{env}} = 5.24^{+0.30}_{-0.29}$| per cent; Cloutier et al. 2021). However, the large envelope mass of GJ 1214 b is easily reconciled with the host star’s high metallicity (0.29 ± 0.12; Newton et al. 2014), which likely facilitated rapid core formation and makes GJ 1214 b a rare gas-enveloped planet around mid-to-late M dwarfs (Ment & Charbonneau 2023). The TOI-1266 planets do not benefit from forming in a similarly metal-rich environment, which adds to the system’s complexity. Although the impact of photoevaporative mass-loss may be less severe around metal-poor stars, which have been argued to exhibit lower levels of stellar activity (See et al. 2021).

5.3.3 Inside–out planet formation

The concept of inside–out planet formation offers a plausible explanation for the unique architecture of the TOI-1266 system, including the existence of its ‘inverted’ inner planet pair and a third outer planet. Inside–out formation features planetary core formation at the pressure maximum of the inner boundary of the disc’s magnetorotational instability dead zone (Chatterjee & Tan 2014). Although the exact location of a dead zone’s inner boundary (DZIB) is sensitive to the disk parameters, it is typically located within ∼1 au (Matsumura & Pudritz 2005; Mohanty, Ercolano & Turner 2013), which is comparable to the expected location of the water snow line in dusty protoplanetary discs around low-mass stars (Mulders et al. 2015). Pebble drift in discs around low-mass T-Tauri stars proceeds faster and at smaller grain sizes than in more massive discs around typical T-Tauri discs (Pinilla et al. 2013). This allows massive planetary cores to form rapidly out of the local pebble isolation mass that is enhanced at the DZIB (Bitsch et al. 2018). In the inside–out formation picture, after its formation at the DZIB, this first planetary core may undergo inward type I migration or, if sufficiently massive, can open a gap in the disc that results in the DZIB retreating to larger orbital separations. In either case, the planetary core and DZIB become spatially decoupled, which then allows for additional pebbles to pile up at the now-cleared DZIB, producing a second planetary core from a depleted pebble reservoir. This process can proceed to form additional planets sequentially.

Many of the features of the TOI-1266 system appear to fit this picture if tuned somewhat finely. While neither planet b nor c are massive enough to induce DZIB retreat, they are both sufficiently massive to undergo type I migration (Burn et al. 2021). The distinct bulk compositions of TOI-1266 b and c may be naturally explained if the inner planet b formed during the earliest stages of the system’s lifetime when the DZIB was located inside of the snow line and the gas surface density was sufficient for substantial H/He accretion. After approaching the local pebble isolation mass, TOI-1266 b could have migrated inwards, therefore leaving space for c to form at the DZIB. If planet c’s formation timescale is sufficiently long such that the disc gas becomes depleted and the luminosity of the pre-main-sequence TOI-1266 drops, then it is conceivable that while planet c may have formed at the same location as b, the environment in which it formed may have been gas-poor and beyond the snow line (Martin & Livio 2012). The existence of the massive outer planet candidate d, if already formed by an independent process, would have also contributed to depleting the reservoir of pebbles available to form planet c. This could naturally explain why the sandwiched planet c is the lowest mass planet in the system (Pritchard et al. 2023).

We note that the reliance on the DZIB, as opposed to another location in the disc, may be critical in explaining the architecture of the TOI-1266 system. This is because our measured planet masses (i.e. |$\lesssim 4.6\, \mathrm{M}_{\oplus }$|⁠) are likely too low to induce pressure maxima that are independent of a DZIB.⁵ The TOI-1266 planets are also too low mass to open gaps and induce DZIB retreat. Consequently, type 1 migration would likely be responsible for moving each planet away from the DZIB, which is required to enable sequential planet formation.

5.4 Three-body dynamical considerations

In Fig. 13, we show that the orbital configuration of the TOI-1266 does not lie close to any low-order two or three-body MMRs. The nominal period ratio of the inner planet pair, P_c/P_b, places it somewhat close to the second-order 5:3 MMR with Δ_bc = 3(P_c/P_b)/5 − 1 = 0.035. However, using the resonance width formula given in equation 5 of Hadden & Lithwick (2018), we determine that, for the RV-measured planets masses and eccentricity constraints from the lack of TTVs in the TESS data, the planet pair is outside of the 5:3 MMR. We reach the same conclusion for the c–d pair, which sits at a similar distance Δ_cd = 0.032 from the 5:3 MMR. Fig. 13 also shows that the system does not sit at any low-order three-body MMR, which are shown as the loci of points where

for the combinations (j, k, j′, k′) = {(3, 1, 3, 1), (2, 1, 2, 1), (3, 1, 5, 2), (2, 1, 5, 2), (5, 2, 3, 1), (5, 2, 2, 1)}. The strength of such three-body MMRs scale, to leading order in eccentricity, as |$e^{k+k^{\prime }-2}$|⁠, and so these combinations represent the strongest nearby three-body MMRs up to order e² in the vicinity of the system’s current orbital configuration. We conclude based on the current data that the system is likely not in a resonant chain, but more precise transit timing measurements may be able to directly confirm whether the two- and three-body resonant angles are librating and therefore definitively rule out a resonant configuration.

Figure 13.

Period ratios of the TOI-1266 planets b, c, and candidate d, in relation to nearby low-order two- and three-body MMRs. Locations of 3:2 and 2:1 first-order two-body MMRs are shown by the solid horizontal and vertical lined. Dashed horizontal and vertical lines indicate the location of the second-order 5:3 MMRs. Grey lines indicate the location of strongest nearby three-body MMRs, occurring where j′n_d + (1′ − j′)n_c = jn_c + (1 − j)n_b, 5n_d − 3n_c = jn_c + (1 − j)n_b, or j′n_d + (1′ − j′)n_c = 5n_c − 3n_b for j, j′ = 2, 3.

Open in new tabDownload slide

6 SUMMARY AND CONCLUSIONS

We have presented an RV plus transit modelling analysis of the TOI-1266 planetary system using multiple years of data from TESS and the HARPS-N spectrograph. Our main conclusions are summarized below:

• Both planets, but TOI-1266 c in particular, exhibit increased transit depths from the TESS primary mission to the first extended mission. Their revised radii are |$R_{p,{\rm b}}=2.62\pm 0.11\, \mathrm{R}_{\oplus }$| and |$R_{p,{\rm c}}=2.13\pm 0.12\, \mathrm{R}_{\oplus }$|⁠. The planet pair does not span the radius valley.

• We measure planetary masses of |$M_{p,{\rm b}}=4.23\pm 0.69\, \mathrm{M}_{\oplus }$|⁠, |$M_{p,{\rm c}}=2.88\pm 0.80\, \mathrm{M}_{\oplus }$|⁠, and find evidence for an outer planet candidate in the system with |$M_{p,{\rm d}}\sin {i}=4.59^{+0.96}_{-0.94}\, \mathrm{M}_{\oplus }$|⁠.

• Despite both period ratios P_c/P_b and P_d/P_c being close to 5:3, we demonstrate that the system likely does not form a resonant chain.

• TOI-1266 b and c are among the lowest density sub-Neptunes known around M dwarfs (⁠|$1.7-3\, \mathrm{R}_{\oplus }$|⁠).

• TOI-1266 b requires a H/He envelope to explain its mass and radius (X_env,b = 5.5 ± 0.7 per cent), while TOI-1266 c is consistent with having a water mass fraction whose value is expected from models of formation beyond the snow line (WMF_c = 59 ± 14 per cent). The water-rich interpretation of planet c is supported by hydrodynamic escape models, although a H/He envelope cannot be ruled out without atmospheric characterization.

• If the aforementioned compositions are true, then TOI-1266 would be the only known M dwarf system whose larger planet has a smaller orbital separation (i.e. an ‘inverted’ system) and that features planets with distinct bulk compositions.

• The system’s unique architecture and distinct planet compositions may be explained by expectations from inside–out sequential planet formation.

The unique architecture of the TOI-1266 planetary system presents unique challenges to planet formation models that have been largely successful at explaining properties of the small planet population around M dwarfs. Progress towards validating the inside-out picture of planet formation, or the development of alternative theories, will rely on more stringent constraints on the planets’ bulk compositions. Given the degeneracies in inferring accurate bulk compositions from small planets from masses and radii alone, comparative atmospheric planetology will be needed to firmly distinguish between water versus H/He-rich compositions and to shed light on the formation pathway of this system.

ACKNOWLEDGEMENTS

We thank Jon Jenkins and the Transiting Exoplanet Survey Satellite (TESS) Science Processing Operations Center (SPOC) team for their investigation of the TESS light curve upon the discovery of the transit depth discrepancy. In particular, the team’s careful reassessment of the effects of crowding on the light curve. We thank Gudmundur Stefánsson for sharing the ARCTIC/Apache Point Observatory light-curve data. We also thank Farzana Meru and Madison VanWyngarden for insightful discussions.

RC acknowledges support from the Natural Sciences and Engineering Council of Canada (NSERC) and the Banting Postdoctoral Fellowship Program administered by the Government of Canada.

MP acknowledges the financial support from the Agenzia Spaziale Italiana-Istituto Nazionale di Astrofisica (ASI-INAF) Addendum no. 2018-24-HH.1-2022 ‘Partecipazione italiana al Gaia DPAC – Operazioni e attività di analisi dati’.

FR is funded by the University of Exeter’s College of Engineering, Maths and Physical Sciences, UK.

This paper includes data collected with the TESS mission, obtained from the Mikulski Archive for Space Telescopes (MAST) data archive at the Space Telescope Science Institute (STScI). Funding for the TESS mission is provided by the NASA Explorer Program. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555.

Based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated by the Fundación Galileo Galilei (FGG) of the Instituto Nazionale di Astrofisica (INAF) at the Observatorio del Roque de los Muchachos (La Palma, Canary Islands, Spain).

The HARPS-N project has been funded by the Prodex Program of the Swiss Space Office (SSO), the Harvard University Origins of Life Initiative (HUOLI), the Scottish Universities Physics Alliance (SUPA), the University of Geneva, the Smithsonian Astrophysical Observatory (SAO), the Italian National Astrophysical Institute (INAF), the University of St Andrews, Queens University Belfast, and the University of Edinburgh.

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.

DATA AVAILABILITY

This paper includes data collected by the TESS mission, which is publicly available from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute (STScI) https://mast.stsci.edu. The HARPS-N data underlying this article are available in the article and in its online supplementary material.

Footnotes

References