Masses, Revised Radii, and a Third Planet Candidate in the “Inverted” Planetary System Around TOI-1266

Is the population of close-in planets orbiting M dwarfs sculpted by thermally driven escape or is it a direct outcome of the planet formation process? A number of recent empirical results strongly suggest the latter. However, the unique architecture of the TOI-1266 system presents a challenge to models of planet formation and atmospheric escape given its seemingly “inverted” architecture of a large sub-Neptune ( 𝑃 𝑏 = 10 . 9 days, 𝑅 𝑝,𝑏 = 2 . 62 ± 0 . 11 R ⊕ ) orbiting interior to that of the system’s smaller planet ( 𝑃 𝑐 = 18 . 8 days, 𝑅 𝑝,𝑐 = 2 . 13 ± 0 . 12 R ⊕ ). Here we present revised planetary radii based on new TESS and diffuser-assisted ground-based transit observations, and characterize both planetary masses using a set of 145 radial velocity measurements from HARPS-N ( 𝑀 𝑝,𝑏 = 4 . 23 ± 0 . 69 M ⊕ , 𝑀 𝑝,𝑐 = 2 . 88 ± 0 . 80 M ⊕ ). Our analysis also reveals a third planet candidate ( 𝑃 𝑑 = 32 . 3 days, 𝑀 𝑝,𝑑 sin 𝑖 = 4 . 59 + 0 . 96 − 0 . 94 M ⊕ ), which if real, would form a chain of near 5:3 period ratios, although the system is likely not in a mean motion resonance. Our results indicate that TOI-1266 b and c are among the lowest density sub-Neptunes around M dwarfs and likely exhibit distinct bulk compositions of a gas-enveloped terrestrial ( 𝑋 env ,𝑏 = 5 . 5 ± 0 . 7%) and a water-rich world (WMF 𝑐 = 59 ± 14%), which is supported by hydrodynamic escape models. If distinct bulk compositions are confirmed through atmospheric characterization, the system’s unique architecture would represent an interesting test case of inside-out sub-Neptune formation at pebble traps.


INTRODUCTION
The radius valley is a stark feature in the exoplanet population between ∼ 1.6 − 1.9 R ⊕ 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 ★ E-mail: ryan.cloutier@mcmaster.caet al. 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 et al. 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 onto low mass cores (Lee & Chiang 2015;Lee & Connors 2021), or by waterrich formation 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 scenario (Cloutier & Menou 2020).The XUV spectrum of the multi-planet host K2-3 cannot explain the low density compositions of its inner planets using TDML models (Diamond-Lowe et al. 2022).A sub-population 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 et al. 2023).The low bulk density of the intensively studied 1.5 R ⊕ planet Kepler-138 d was shown to be inconsistent with TDML and requires a substantial volatile mass fraction suggestive of waterrich formation 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 (  = 10.9 days,  , = 2.37 R ⊕ ) alongside an outer super-Earth (  = 18.8 days,  , = 1.56 R ⊕ ; D20).Multi-planet system architectures that contain a likely terrestrial planet (i.e.≲ 1.6 R ⊕ ) whose orbit is greater than the system's sub-Neptune are somewhat rare (∼ 35%; 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.

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-minute 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 Figure 1 and the phase-folded transit light curves from the PM and EM are included in Figure 2.
In this paper we use the PDCSAP light curves in addition to custom light curve extractions from the 2-minute TESS Target Pixel Files (TPFs) and the Full Frame Images (FFIs).The latter has cadences of 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).30 and 10 minutes 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-minute light curve from the TPFs using the lightkurve package (Lightkurve Collaboration et al. 2018;Astropy Collaboration et al. 2022;Ginsburg et al. 2019).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 Figure 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 minutes in the PM and 10 minutes 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 Figure 2.
A preliminary inspection of transit light curves in Figure 2 revealed significant transit depth discrepancies between the PM and EM.We postpone our investigation of these anomalies until Section 4.5.

APO/ARCTIC diffuser-assisted photometry
Figure 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.

Palomar/WIRC diffuser-assisted photometry
We observed one full transit of each of TOI-1226 b and c in the -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 × 8. ′ 7 with a plate scale of 0. ′′ 25 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-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, flatfielded, 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-25 pixels and selected the aperture that minimized the 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 et al. 2019, 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 Figure 2.
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 1 , 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 Figure 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  -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.
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 2 /s 2 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 6.

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 et al. 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 (Figure 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  rot ∼ 45 days 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  ′ 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  BP ,  RP , and 2MASS   -band magnitudes and use the empirical relations for M dwarfs from Mann et al. (2013Mann et al. ( , 2015Mann et al. ( , 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 SEDfitting 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 2, along with all other stellar parameters.the full RVs, the chromatic   RVs, and the following activity indicators: , dLW, FWHM, and BIS.The corresponding false alarm probability (FAP) of each time series is included in Figure 4.

Generalized Lomb-Scargle periodograms
Aside from the aliased signals close to 1-day (Dawson & Fabrycky 2010), the GLS of the full RV time series exhibits one periodicity with FAP < 1% at 10.9 days, 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 days; FAP = 95%), plus a peak at 32.3 days that is unrelated to either known transiting planet (FAP = 4%).We investigated the origin of the latter by assessing its chromatic dependence in the  -bands and its presence in the GLSs of various activity indicators.The chromatic RVs are not informative because the -band RVs are too low S/N, and while the 32.3-day signal does exist in the  and -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 -band to the -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-day signal until Section 4.3.
Turning to the activity indicators, we do not see a low FAP peak close to 32.3 days 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% is exhibited by our  time series at 44.6 days. is a strong indicator of M dwarf chromospheric activity at optical wavelengths (Lafarga et al. 2021) and we interpret the 44.6-day signal as the star's rotation period, or a harmonic thereof.The lack of a low FAP peak at 32.3 days, or 32.3/2 ≈ 16.2 days, 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.

Global RV + Transit Modeling
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 dataset.
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% 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 minutes and 16 minutes for planets b and c, respectively, but these measurements represent only marginal TTV detections at < 3.The corresponding TTV masses are poorly constrained:  , = 13.5 +11.0 −9.0 M ⊕ and  , = 2.2 +2.0 −1.5 M ⊕ (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  0 , scalar jitter , quadratic limb-darkening parameters { 1 ,  2 }, orbital periods , planet-star radius ratios   / ★ , impact parameters , and rescaled orbital eccentricities  and arguments of periastron The RV component includes the additional parameters: the systemic RV  0 , scalar jitter  RV , and each planet's RV semiamplitude .Our adopted priors are summarized in Table 4.
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 days of insofar unknown origin (c.f. Figure 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  time series as it exhibits the periodic signal at ∼ 44.6 days 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  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 et al. 2019, 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 dataset 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 (SHO).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) (Nelson et al. 2020).
We use the PyMC3 Markov Chain Monte Carlo (MCMC) package (Salvatier et al. 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 103 burn-in steps and 10 4 posterior draws.We report maximum a-posteriori (MAP) point estimates from each parameter's marginalized posterior PDF in Table 4. Parameter uncertainties are derived from each marginalized posterior's 16 th and 84 th percentiles unless noted otherwise.
The detrended transit light curves were shown in Figure 2 while the individual RV components are in Figure 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 days 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.Figure 6 depicts the phase-folded RV time series for each planet.

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-day 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 , we estimate its Bayesian model evidence Z  using the Perrakis estimator, for which we use the model's joint posterior from our MCMC calculations as an importance sampler (Perrakis et al. 2014).We then compute Bayes factors, or evidence ratios, to determine the favoured model given our RV data.Our results are reported in Table 3.
We find that in both our 2p and 3p models, the inclusion of a GP is favoured. 3Despite 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, Z 3 +GP /Z 2 +GP = 10 6.4 and Z 3  /Z 2  = 10 5.6 .This suggests that the 32.3-day 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 4.
- * Marginalized posterior distribution from training on our dLW time series.† / ★ priors are based on the planet's period, and the stellar mass and radius.‡ 95% upper limit.§ For the candidate planet TOI-1266 d, we only measure the minimum planet mass  , sin  rather than  , .¶ Equilibrium temperature calculations assume zero albedo and perfect heat redistribution.△ Transmission Spectroscopy Metric (Kempton et al. 2018).

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 modeling 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 fo-  cus on deriving eccentricity constraints by modeling TOI-1266 b's transit times only.We do so by adopting a linear model approach as described in Linial et al. (2018).We are sensitive to the magnitude of 'combined eccentricity' in this model,  4. We conclude that |Z , | < 0.08 (0.11) at 95% (99%) confidence, which is more stringent than our eccentricity constraints derived from the RVs alone.

Transit Depth Discrepancy
As introduced in Section 2.1 and illustrated in Figure 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 artifacts 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.

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.

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.

Residual Artifacts 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.
Figure 7 and Table 5 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.PDC-SAP, 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 (Figure 7).We followed-up with additional transit observations of TOI-1266 b and c in the -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.

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 5), 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 timeseries using an adapted version of SOAP 2.0 (Dumusque et al. 2014;Wurmser et al. in prep.).We explored a grid of  = 43, 890 plage sizes (0 − 0.65  ★ ), 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 et al. 2022) such that we assume TOI-1266's rotation axis is perpendicular to our line-of-sight.We adopt  rot = 30 days as these calculations were performed before the identification of the probable rotation period at 44.6 days.Fortunately, our preliminary tests sampled orbital periods from 30-60 days and revealed that our results are insensitive to the exact choice of  rot .
We compare our grid of plage models to our data in each of our four observing seasons of spectroscopic follow-up (see Figure 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.Figure 8 shows the results of this exercise.We find that small, cool plages at high latitudes ( plage / ★ ≲ 0.17, Δ ≲ 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 Dec 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.
Our finding that plages on TOI-1266 are localized at high latitudes is consistent with the findings of Morales et al. (2010) who analyzed 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.  = 0.77;   = 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   / ★ and unocculted ARs to recover robust planetary radii from transit depth measurements.The basis of the framework is that while the true value of   / ★ 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   / ★ values if the transit light curve has sufficient S/N and limb-darkening is wellconstrained.In light of the possibility of elevated plage coverage in observing season 3, we attempted to use this framework to measure each planet's   / ★ 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 toward the TESS EM epochs, this does not have a demonstrable impact on the observed transit depths.

Stochastic Effects Between the PM and EM TESS Data
We consider whether there are stochastic effects in either the PDC-SAP 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  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 days 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.
Figure 9 compares the resulting distributions of transit depth differences, Δ =  rec −  inj , for the PM and EM and for each set of injected ephemerides.We treat Δ as a measure of the accuracy with which transit depths are recovered.We find that when planets are injected with random ephemerides, the Δ 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 5).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 Δ 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.
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).

Conclusions Regarding the Transit Depth Discrepancy
We conclude that variable flux dilution, orbital precession, and residual artifacts 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.

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-50 days 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  <   = 18.8 days compared to the TESS team, and because the transit probability drops below 1% at ∼ 50 days.TOI-1266 b and c are remarkably coplanar (Δ < 0.1 • ) and a third coplanar planet would no longer exhibit a transiting configuration beyond 27.1 days.We also conducted a focused TLS search between 31.9-32.7 days 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 tify our sensitivity to transiting planets.We sampled orbital periods logarithmically between 19-50 days, uniform times of mid-transit, planet radii logarithmically between 0.5 − 3 R ⊕ , and impact parameters linearly between 0 − 1.We constructed 10 4 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% of the injected period.
The top panel of Figure 10 depicts our detection sensitivity to transiting planets.We do not correct for the geometric probability given the preference of compact multi-planet systems to exhibit small mutual inclinations ≲ 1 − 2 • (Ballard & Johnson 2016).We confidently rule out transiting planets down to 1.4 R ⊕ within 23 days but lack sensitivity to planets within the habitable zone (HZ) bounded by the water-loss and maximum greenhouse limits (36-95 days; Kopparapu et al. 2013).

Radial Velocity Planet Search
We conducted a similar exercise to assess our detection sensitivity to planets in our HARPS-N RV timeseries.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 ±/ ★ , where  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 − 20 M ⊕ .We constructed 10 4 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 Figure 5).We deem successful recoveries of injected planets to be those that exhibit a FAP ≤ 0.1% in the GLS within 5% 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 Figure 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 M ⊕ within 20 days and planets down to 4 M ⊕ 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.

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  , = 2.62 ± 0.11 R ⊕ and  , = 2.13 ± 0.12 R ⊕ .We also measure the masses of both planets using data from the HARPS-N spectrograph ( , = 4.23±0.69M ⊕ and  , = 2.88 ± 0.80 M ⊕ ).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  , sin  = 4.59 +0.96  −0.94 M ⊕ (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   = 1.3 +0.3 −0.2 F ⊕ .The left panel in Figure 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 R ⊕ .If we assume that both planets have Earth-like cores with core mass fractions CMF of 33%, and are enveloped in H 2 gas, then at their present-day equilibrium temperatures of 425 K and 354 K, respectively, we infer envelope mass fractions of  env, = 5.5 ± 0.7% and  env, = 3.5 ± 0.8%.
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;Piaulet et al. 2023;Cherubim 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 multi-phase and non-isothermal equation of state for water, we find that WMF  > 80% at 95% confidence4 and WMF  = 59 ± 14%.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  at  eq = 354 K.
We note that at solar composition, the WMF beyond the water snow line is expected to be at most 50% (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 4.
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 Figure 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 gasenveloped.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 multiplanet systems are not "inverted" (i.e.0.65 ± 0.4% 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.

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.∼ 1.6 R ⊕ ; 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 prep.).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. Figure 11).We therefore do not use the multi-transiting 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.

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  core can accrete is the maximally cooled isothermal mass  iso (Lee & Chiang 2015).This limit leads to the possibility of a primordial radius valley wherein low mass cores (≲ 1 − 2 M ⊕ ; Lee et al. 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 et al. 2022).
We calculate  iso for TOI-1266 b and c following the framework and assumptions outlined in Lee et al. (2022).Assuming in-situ formation, we find that for gas accretion timescales ≥ 1 Myr, both planets are capable of accreting  iso =  iso / core greater than their observed  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.

Hydrodynamic Atmospheric Escape Modeling
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 prep) wherein the XUVdriven 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  eq, .
The star remains in the saturated activity regime for 500 Myr with  XUV / bol = 10 −3 before decaying as a powerlaw, 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 et al. 2009).We do not explicitly model gas accretion onto the planetary cores and instead consider a grid of initial envelope mass fractions  env,0 ∈ [0.01, 0.1].
We evolve the atmospheres under the aforementioned effects and present our results in Figure 12.We find that neither planet is completely stripped after 5 Gyrs for all but the smallest  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 XUVdriven photoevaporation, either solely or alongside core-powered mass loss, because its envelope is efficiently lost, leaving behind  env,b < 5.5%.Meanwhile, the observed radius of TOI-1266 c can be recovered when  env,0,c ≳ 0.03.
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.Because a H/He component is needed to explain the mass and radius of planet b, our atmospheric escape model therefore rules out the possibility that planet c is 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.

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 Myrs.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 2 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).

The Lowest 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 days, RV-derived planetary masses appear to be systematically larger than those measured with TTVs, across all radii < 8 R ⊕ (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 loworder 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 et al. 2016), then TOI-1266 b and c may be tracing such a formation pathway given their low period ratio of   /  = 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 the lowest density planets transiting M dwarfs between 1.7 − 3 R ⊕ (Figure 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 prep.).

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 et al. 2015).With  env, = 5.5 ± 0.7%, 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 ( env = 5.24 +0.30  −0.29 %; 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).

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 disk'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 et al. 2013), which is comparable to the expected location of the water snow line in dusty protoplanetary disks around low mass stars (Mulders et al. 2015).Pebble drift in disks around low mass T-Tauri stars proceeds faster and at smaller grain sizes than in more massive disks around typical T-Tauri disks (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 disk 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 inward, therefore leaving space for c to form at the DZIB.If planet c's formation timescale is sufficiently long such that the disk gas becomes depleted and the luminosity of the pre-mainsequence 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. in prep).
We note that the reliance on the DZIB, as opposed to another location in the disk, may be critical in explaining the architecture of the TOI-1266 system.This is because our measured planet masses (i.e.≲ 4.6 M ⊕ ) are likely too low to induce pressure maxima that are independent of a DZIB5 .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.

Three-body Dynamical Considerations
In Figure 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,   /  , places it somewhat close to the second-order 5:3 MMR with Δ  = 3(  /  )/5 − 1 = 0.035.However, using the resonance width formula given in Equation 5of 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 Δ  = 0.032 from the 5:3 MMR.

SUMMARY AND CONCLUSIONS
We have presented an RV plus transit modeling 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  , = 2.62 ± 0.11 R ⊕ and  , = 2.13 ± 0.12 R ⊕ .The planet pair does not span the radius valley.• We measure planetary masses of  , = 4.23 ± 0.69 M ⊕ ,  , = 2.88 ± 0.80 M ⊕ , and find evidence for an outer planet candidate in the system with  , sin  = 4.59 +0.96  −0.94 M ⊕ .• Despite both period ratios   /  and   /  being close to 5:3, we demonstrate that the system likely does not form a resonant chain.
• TOI-1266 b and c are the lowest density sub-Neptunes known around M dwarfs (1.7 − 3 R ⊕ ).
• TOI-1266 b requires a H/He envelope to explain its mass and radius ( env, = 5.5 ± 0.7%) 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  = 59 ± 14%).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 toward 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.

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

Figure 3 .
Figure 3.An illustrative example of our line-by-line (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 Δ 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).

Figure 4 .
Figure 4. GLS periodograms of our HARPS-N window function, the full LBL RVs, the chromatic LBL RVs in the  -bands, and the spectroscopic activity indicators  , dLW, FWHM, and BIS.Left column: the GLS periodogram of the time series labeled 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 days in the full RVs and   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% FAP level.Right column: the FAP as a function of GLS power corresponding to each periodogram.

Figure 5 .
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   time series.

Figure 6 .
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).
which depends on the planet eccentricities   and arguments of periastron   .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

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

Figure 8 .
Figure8.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.

Figure 9 .
Figure 9. Cumulative distributions of transit depth measurement accuracies Δ 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.

Figure 10 .
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 (  = 18.8 days) and out to 50 days where our sensitivity becomes consistent with 0%.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 days.

Figure 11 .
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%, 33% (i.e.Earth-like), and 100%, and water-rich models from Aguichine et al. 2021 with water mass fractions (WMF) of 10%, 20%, and 50% 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 thermally driven mass loss (TDML) and gas-depleted formation (GDF) (Cloutier & Menou 2020).Multi-planet systems are connected by the dashed lines in both panels.

Figure 12 .
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  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 9 years in our model of TOI-1266 b with  env,0 = 0.01 occurs when complete atmospheric loss is achieved, leaving behind a bare core.
Figure 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

Table 1 .
Summary of RV time series

Table 4 .
Priors and posterior point estimates of the TOI-1266 planetary system model parameters.

Table 5 .
Comparative transit parameters across different observing facilities and epochs.
* WIRC transit times are excluded here and are reserved for a forthcoming detailed TTV analysis(Greklek-McKeon et al. in prep.)

Table 6 .
HARPS-N Spectroscopic Time Series of TOI-1266.For conciseness, only a subset of rows are depicted here to illustrate the table's contents.The entirety of this table is provided in the arXiv source code.