TOI-5126: A hot super-Neptune and warm Neptune pair discovered by $\textit{TESS}$ and $\textit{CHEOPS}$

We present the confirmation of a hot super-Neptune with an exterior Neptune companion orbiting a bright (V = 10.1 mag) F-dwarf identified by the $\textit{Transiting Exoplanet Survey Satellite}$ ($\textit{TESS}$). The two planets, observed in sectors 45, 46 and 48 of the $\textit{TESS}$ extended mission, are $4.74^{+0.16}_{-0.14}$ $R_{\oplus}$ and $3.86^{+0.17}_{-0.16}$ $R_{\oplus}$ with $5.4588385^{+0.0000070}_{-0.0000072}$ d and $17.8999^{+0.0018}_{-0.0013}$ d orbital periods, respectively. We also obtained precise space based photometric follow-up of the system with ESAs $\textit{CHaracterising ExOplanets Satellite}$ ($\textit{CHEOPS}$) to constrain the radius and ephemeris of TOI-5126 b. TOI 5126 b is located in the"hot Neptune Desert"and is an ideal candidate for follow-up transmission spectroscopy due to its high predicted equilibrium temperature ($T_{eq} = 1442^{+46}_{-40}$ K) implying a cloud-free atmosphere. TOI-5126 c is a warm Neptune ($T_{eq}= 971^{+31}_{-27}$ K) also suitable for follow-up. Tentative transit timing variations (TTVs) have also been identified in analysis, suggesting the presence of at least one additional planet, however this signal may be caused by spot-crossing events, necessitating further precise photometric follow-up to confirm these signals.


INTRODUCTION
Super-Neptunes are a class of planets potentially bridging the gap between Neptune and Jovian planets 1 (4-8  ⊕ ).Despite being of similar sizes, they exhibit a significant diversity in mass, almost uniformly distributed between 6 and 135  ⊕ (Akeson et al. 2013;Petigura et al. 2017).Given the absence of such planets in the Solar system, studying the population provides extra insights on the versatility of the planet formation process around other stars.
An archetypal example of the super-Neptune class is GJ 436 b (Gillon et al. 2007).It has a Neptune-like radius coupled with a relatively high density.It also has an inferred rocky element core and substantial H/He envelope (Fortney et al. 2007).Its characteristics appear to deviate from the theorised ≈ 5  ⊕ threshold needed for runaway accretion (Lee & Chiang 2016a), highlighting our lack of ★ Email: tyler.fairnington@usq.edu.au† ARC DECRA Fellow ‡ ARC DECRA Fellow § 51 Pegasi b Fellow 1 in the literature, these planets are also often referred as Sub-Saturns understanding about the formation channels from which this class originates.Kepler and the subsequent K2 mission (Borucki et al. 2010;Howell et al. 2014) enabled the first population studies of these planets.Petigura et al. (2017) proposed that super-Neptunes form in gas-depleted disks (Lee et al. 2014) or from planet-planet mergers (Petigura et al. 2017).Both of these formation scenarios lead to heavy elements in the planets atmosphere, making these theories testable with bright systems hosting super-Neptunes.
However, the majority of confirmed super-Neptunes discovered by the Kepler and K2 missions orbit around relatively faint stars, posing challenges for detailed atmosphere characterisation.The allsky Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015) mission aims to survey the entire sky for candidates around some of the brightest nearby stars, providing the best opportunity to further our knowledge of super-Neptunes.TESS has already identified several unique super-Neptunes, including LTT 9779 b (Jenkins et al. 2020a) and TOI-849 b (Armstrong et al. 2020).Both reside in the 'hot Neptune desert' (Mazeh et al. 2016), creating further challenges to their formation and evolution models.Planets in this close-in highly irradiated period space are usually thought to be most susceptable to photoevaporation (Owen & Jackson 2012;Jin et al. 2014;Lopez 2017) and tidal stripping (Matsakos & Königl 2016;Owen & Lai 2018).LTT 9779 b is an ideal target for follow-up atmospheric studies, with a Spitzer phase curve already suggesting a high-metallicity atmosphere (Crossfield et al. 2020).Efforts to understand this population can be most fruitful via multiplanet systems, where comparative planetology is possible.
In this paper, we present the discovery of a multi-transiting system around the F-dwarf ( eff = 6150 +130 −100 K) TOI-5126.The system, first discovered by TESS, hosts a hot super-Neptune located in the 'hot Neptune Desert', and an exterior warm Neptune.The brightness of TOI-5126 and the high predicted equilibrium temperatures of planet b (1442 +46 −40 K) and c (971 +31 −27 K) make these planets ideal for future follow-up, including atmospheric comparative planetology studies.
In section §2, we present the photometric, imaging and spectroscopic observations that led to the discovery and confirmation of the TOI-5126 system.Section §3 discusses the derived stellar parameters of TOI-5126, as well as the global modeling and analysis of the system.The results of our analysis are presented in section §4, including validation of TOI-5126 b and TOI-5126 c, transit timing variations (TTVs), and future follow-up opportunities.

TESS
TOI-5126 (TIC 27064468) was observed by TESS nearly continuously on both the 2 min cadence target pixel stamps and the 10 min FFIs in sectors 45 (Camera 3 CCD 4) and 46 (Camera 2 CCD 2), with an additional sector (sector 48, Camera 1 CCD 4) a month later.The observations were taken from UT 2021-Nov-06 to UT 2022-Feb-26, for a total of ≈76 days.
In sector 46, systematic errors were observed in the first 13 hours of data due to a pointing adjustment 2 .To mitigate these errors, we utilised quality flags to mask out scattered light and other systematic effects in the light curves .
The light curves of TOI-5126 were processed by both the TESS Science Processing Operations Center (SPOC; Jenkins et al. 2016), located at NASA Ames Research Center.The SPOC conducted a transit search of Sector 45 2-min cadence data on UT 2021-Jan-09 with an adaptive, noise-compensating matched filter (Jenkins et al. 2002(Jenkins et al. , 2010(Jenkins et al. , 2020b)), producing a Threshold Crossing Event (TCE) for which an initial limb-darkened transit model was fitted (Li et al. 2019) and a suite of diagnostic tests were conducted to help make or break the planetary nature of the signal (Twicken et al. 2018).Transit signatures were also detected in a search of Full Frame Image (FFI) 10-min cadence data by the Quick Look Pipeline (QLP) at MIT (Huang et al. 2020).The TESS Science Office (TSO) reviewed the vetting information and issued an alert for TOI-5126 b on UT 2022-Jan-25 based and an alert for TOI-5126 c on UT 2022-Feb-08 (Guerrero et al. 2021) based on QLP data.The signals were repeatedly recovered as additional observations were made including in both Sectors 46 and 48, as well as SPOC multi-sector searches (42)(43)(44)(45)(46).The transit signatures again passed all the diagnostic tests presented in the Data Validation reports.Specifically, difference images from the Sectors 42-46 data located the catalog location of 2 https://archive.stsci.edu/missions/tess/doc/tess_drn/tess_sector_46_drn66_v02.pdf TOI-5126 b and TOI-5126 c within 9.0 ± 5.0 ′′ and 5.2 ± 5.8 ′′ of the transit source.The latest SPOC results captured signals for planets b (Signal-to-Noise (SNR): 26.1;Multiple Event Statistic (MES): 19.9) and c (SNR: 10.4; MES: 9.9). Te transit events of TOI-5126 b were also identified by the Planet Hunters TESS citizen science project (Eisner et al. 2020) and uploaded to exoFOP as communitiy objects of interest (CTOI) on UT 2022-Jan-12.
For our analysis, a modification of the deblended SPOC Simple Aperture Photometry (SAP) light curves (Twicken et al. 2010;Morris et al. 2020) were used in modeling (see Figure 1).Stellar variability was accounted for by simultaneously fitting for a polynomial trend around each transit as part of our global analysis.Two transits were manually excluded due to them occurring near a data gap.Our ephemeris and radius precision did not noticeably suffer from these excluded transits.The details of the global fit will be presented in section §3.2.

CHEOPS
To refine the radius and ephemeris of TOI-5126 b, we obtained spacebased photometric observations with CHEOPS (Benz et al. 2021).Observations were attained through the CHEOPS Guest Observers Programme (AO3-023).A first observation of TOI-5126 b occurred on 2022-02-16 21:44 UT for a duration of 11.72 hours.The preliminary data was processed through the CHEOPS Data Reduction Pipeline (version 13.1;Hoyer et al. 2020) on 2023-Feb-17 18:46 UT.The visit consisted of seven CHEOPS orbits, with a four orbit pre-ingress and one orbit post-egress out-of-transit baseline.An additional visit of TOI-5126 b occurred on 2023-Feb-27 23:55 UT in an attempt to identify an ingress/egress with high efficiency, resulting in an longer 14.14 hour observation, with two orbits pre-ingress and five orbits post-egress.We obtained a third visit of TOI-5126 b on 2023-Mar-27 04:44 where we successfully caught both the ingress and egress over an observation duration of 10.96 hours.The full details of observations are listed in Table 1.
The raw CHEOPS light curves were detrended simultaneously with the transit modeling to best incorporate and propagate the uncertainties associated with the instrument into our final system parameters.Observations with poor data quality flags, and those that fail an iterative 5- outlier rejection, were rejected from subsequent analyses.To correct for the CHEOPS systematics, we followed the trend model in the Pycheops package (Maxted et al. 2022).This model fits for spacecraft and environmental variations, including the roll angle, , the background flux, bg, the contamination in the aperture, contam, smearing correction, smear and the x and y point spread function (PSF) centroid coordinates, dx & dy.The raw CHEOPS light curves were detrended by fitting a least-squares model to the light curve residuals, which provided us with the best parameters for the trend model.These parameters were then applied to the trend function and subtracted from the raw flux to provide us with the detrended CHEOPS light curves , seen in Figure 3.

Ground based photometric follow up
Due to the coarse angular resolution of both TESS and CHEOPS, eclipsing binary stars within their apertures may mimic a transit event on the target star.To confirm the on-target detection of TOI-5126, we obtained ground-based photometric observations of both planets.We used the TESS Transit Finder, which is a customized version of the Tapir software package (Jensen 2013), to schedule our transit observations.TOI-5126 b was observed in the   filter by two    The images were calibrated by the standard LCOGT BANZAI pipeline (McCully et al. 2018) and differential photometric data were extracted using AstroImageJ (Collins et al. 2017).We opted to perform our own detrending of all LCOGT light curves simultaneous to the global fit.The trend function for the light curves were represented by a second-order polynomial, with an additional term counting for the full width half maximum of the PSF.As with the CHEOPS trend model, the coefficients of the detrending vectors were derived from a least squares fit that was performed at each iteration of the global modeling process.

Reconnaissance Spectroscopy
We obtained a series of spectroscopic observations to check for contaminating spectroscopic blend scenarios, and to better constrain the spectroscopic atmospheric parameters of the target star.These reconnaissance spectroscopic observations are detailed in Table 3 and seen in Figure 5.
We obtained three spectra of TOI-5126 using the Tillinghast Reflector Echelle Spectrograph (TRES; Fűrész et al. 2008) on the 1.5 m telescope at the Fred Lawrence Whipple Observatory (Mt.Hopkins, Arizona).TRES has a resolving power of R ≃ 44,000, with a wavelength coverage of 3850 to 9096 Å.The observations were obtained with an average exposure SNR per resolution element of 38 at 5110 Å. Radial velocities are derived from a multi-order cross correlation analysis as per Quinn et al. (2012), and spectroscopic parameters were derived via the Spectroscopic Parameter Classification tool (Buchhave et al. 2010).Observations yielded stellar parameters for effective temperature ( eff = 6217 ± 50 K), stellar surface gravity (log  = 4.3 ± 0.1 cgs), metallicity ([Fe/H] = 0.17 ± 0.08 dex) and projected rotation velocity ( sin  ★ = 14.8 ± 0.5 km s −1 ).
In addition, we also obtained five observations of TOI-5126 with the CHIRON high resolution spectrograph, located on the 1.5 m SMARTS telescope at Cerro Tololo Inter-American Observatory (CTIO), Chile (Tokovinin et al. 2013).Observations were obtained using the fiber fed image slicer, yielding a resolution of  ∼ 80, 000.The spectra were extracted and reduced via the official CHIRON pipeline (Paredes et al. 2021).Radial velocities were determined from the modeling of their least-squares deconvolution profiles, generated against a non rotating synthetic template (Zhou et al. 2021).In addition, we also model the rotational broadening of the target as a combination of the rotational, instrument, and macroturbulent broadening kernels.We derive broadening velocities of  sin  ★ = 13.2 ± 0.5 km s −1 and  macro = 6.8 ± 0.5 km s −1 for TOI-5126.
We find no large radial velocity variations at the > 50 m s −1 level, or line profile deviations, as is expected for a system of small planets.The results indicate no obvious spectroscopic blending are present in the system.

High resolution imaging
As part of our standard process for validating transiting exoplanets to assess the the possible contamination of bound or unbound companions on the derived planetary radii (Ciardi et al. 2015), we observed TOI-5126 with high-resolution near-infrared adaptive optics (AO) imaging at Palomar Observatory and in the optical with speckle imaging at WIYN and SAI.The infrared observations provide the deepest sensitivities to faint companions while the optical speckle observations provide the highest resolution imaging making the two techniques complementary.Combined with Gaia, the high resolution imaging observations find no evidence for additional stellar companions within the system.

Palomar 5m/PHARO
The Palomar Observatory observations were made with the PHARO instrument (Hayward et al. 2001) behind the natural guide star AO system P3K (Dekany et al. 2013) on 2022 Feb 13 UT in a standard 5-point quincunx dither pattern with steps of 5 ′′ in the narrow-band  −  filter (  = 2.1686; Δ = 0.0326 m).Each dither position was observed three times, offset in position from each other by 0.5 ′′ for a total of 15 frames; with an integration time of 5.7 seconds per frame, the total on-source time was 85.5 seconds.PHARO has a pixel scale of 0.025 ′′ per pixel for a total field of view of ∼ 25 ′′ .The science frames were flat-fielded and sky-subtracted.The reduced science frames were combined into a single combined image with a final resolution of 0.97 ′′ FWHM.
To within the limits of the AO observations, no stellar companions were detected.The sensitivities of the final combined AO image were determined by injecting simulated sources azimuthally around the primary target every 20 • at separations of integer multiples of the central source's FWHM (Furlan et al. 2017;Lund & Ciardi 2020).The brightness of each injected source was scaled until standard aperture photometry detected it with 5− significance.The resulting brightness of the injected sources relative to TOI 5126 set the contrast limits at that injection location.The final 5− limit at each separation was determined from the average of all of the determined limits at that separation and the uncertainty on the limit was set by the rms dispersion of the azimuthal slices at a given radial distance.Figure 6 shows no stellar companion with a difference in magnitude (Δmag) of 5.4 at a separation ∼0.3 ′′ , and all companions ruled out in a contrast of Δmag of 9 from ∼2 ′′ onwards.No sources are present in the field-of-view with an upper limit on the angular separation of 4 ′′ .

WIYN 3.5m/NESSI
Further imaging observations were made using the NN-Explore Exoplanet Stellar Speckle Imager (NESSI; Scott et al. 2018) on the WIYN 3.5 m telescope at the Kitt Peak National Observatory (KPNO) in Arizona.The observation and data reduction of NESSI can be found in Howell et al. (2011).Images were taken at 832nm on UT 2022-Apr-18.The 5- sensitivity curve, imaged in Figure 6, outline a Δm of 6 from 0.96 ′′ onwards.The inset image (1.2 ′′ angular separation) also shows no stellar neighbours within the Δm of 6.

SAI 2.5m/Speckle Polarimeter
A third direct image of TOI-5126 was taken using the Speckle Polarimeter on the 2.5 m Sternberg Astronomical Institute (SAI; Safonov et al. 2017).Observations were made in the   filter (≈ 660 − 825).Observations were in agreement with the prior direct imaging observations, reaching a contrast of Δm of 7 from 1 ′′ .

Gaia Assessment
In addition to the high resolution imaging, we have utlized Gaia to identify any wide stellar companions that may be bound members of the system.Typically, these stars are already in the TESS Input Catalog and their flux dilution to the transit has already been accounted for in the transit fits and associated derived parameters.Based upon similar parallaxes and proper motions (e.g., Mugrauer & Michel 2020, 2021;Mugrauer et al. 2022), there are no additional widely separated companions identified by Gaia.
Additionally, the Gaia DR3 astrometry provides additional information on the possibility of inner companions that may have gone undetected by either Gaia or the high resolution imaging.The Gaia Renormalised Unit Weight Error (RUWE) is a metric, similar to a reduced  2 , where values that are ≲ 1.4 indicate that the Gaia astrometric solution is consistent with the star being single whereas  In addition the GAIA DR3 non-single-star flag is no present for the system, implying further that this is a single star system.

Stellar Parameters
Stellar parameters for TOI-5126 were obtained by fitting the star's broadband spectral energy distribution (SED) to stellar atmosphere models using the Python package astroARIADNE (Vines & Jenkins 2022).We used TRES spectroscopic data to define priors on the stellar effective temperature ( eff ) (with its uncertainty inflated to 150 K) and metallicity (Fe/H), while distance constraints were derived from Gaia DR3 parallax (Gaia Collaboration et al. 2022).We used the default prior for Extinction (  ) which constrains an upper limit using the SFD galactic dust maps from Schlegel et al. (1998).
astroARIADNE uses a Bayesian Model Averaging (BMA) approach, fitting multiple models3 accounts for model-specific systematic biases.Each model is weighed individually for a combined weighted and averaged posterior probability.This is then fed into the isochrone fitting package.In SED fitting, we used Tycho-2 (Høg et al. 2000) B and V, Gaia DR3 (Fabricius et al. 2021) G, Bp and Rp, Two-Micron All-Sky Survey (2MASS) (Skrutskie et al. 2006) J, H, K, and Wide-field Infrared Survey Explorer (WISE) (Wright et al. 2010)  1 and  2 bands, with an uncertainties floor applied to all bands based on Eastman et al. (2019).The weighted average results obtained are  eff = 6150 +130 −100 K, log  = 4.31 +0.24 −0.23 cgs,  ★ = 1.249 +0.040 −0.038  ⊙ and [Fe/H] = 0.18 +0.13 −0.13 dex.To identify the rotation period of TOI-5126, we passed the raw deblended SPOC light curve through a Lomb-Scargle periodogram.The highest power of the entire TESS set of observations was constrained to be 4.602 +0.071  −0.067 days.We used the FWHM of the highest  power to derive our uncertainties.We also measured the Lomb-Scargle peaks of each sector separately, finding a scatter in the derived rotation period of  = 0.066 days, consistent with the uncertainty estimated from the peak FWHM.However, the actual uncertainty of the rotation period may be much larger and is difficult to quantify with the short observational baseline available.The periodogram can be seen in Figure 8.Given the system has a measured  sin  ★ of ∼ 14 km s −1 , the inclination of the star is consistent with 90 degrees.

Global Modeling
To characterise the physical parameters of the TOI-5126 system, along with their associated uncertainties, we performed a Markov Chain Monte Carlo (MCMC) fit using the publicily available software package emcee (Foreman-Mackey et al. 2013).We opted to not include the 2022-Apr-01 CTIO observation of TOI-5126 b in our fit because it only captured an ingress.Detrending was performed simultaneously to the transit fits for each photometric observation.During each iteration, we subtracted a transit model created from the batman package (Kreidberg 2015) from the observed data.Next, we fit a trend model to the residuals using least-squares optimization, and compared the best trend model and the residuals to determine the log likelihood for the sampler.Section 2 contains further information on the individual trend models.
The global modeling process involved fitting for the two planets simultaneously.Our fit included the period, epoch, the planet-to-star radius ratio,   / ★ and impact parameter,  of both planets.Eccentricities were initially fixed at 0 as we assumed a circular solution.The limb darkening coefficients for each unique photometric filter, denoted as  1, and  2, , were included in the global fit.These coefficients were parameterised as per Kipping (2013).Additionally, we included a scaling factor for each observation as a free parameter in the fit,   .This accounts for differing levels of background contamination due to the varying apertures of the individual instruments.This parameter is neccessary as the large PSFs of both TESS and CHEOPS allow background stars to contribute to the measured flux, resulting in differing total fluxes to that of the LCOGT light curves .We also opted to fit for  ★ and  ★ to simultaneously constrain them with the planet parameters.
To constrain the period and transit times of both planets, we imposed uniform priors with boundaries at the 5- upper and lower uncertainty provided by the SPOC detection pipeline.For TESS and CHEOPS, we allow the  values to uniformly explore between boundaries of 0.9 to 1.1, this is due to the brightest contaminating star in TESS/CHEOPS aperture account for roughtly 5% of the light of the target star.Since the LCOGT aperture excluded all neighboring stars resolved by the Gaia catalog, we tightly constrain   to only uniformly vary between 0.99 and 1.01.For  ★ and  ★ , we imposed a Gaussian prior derived from the results of the SED fitting obtained in Section 3.1.
We first assumed the orbits of the two planets to be circular.We also ran a set of global modeling allowing non-zero eccentricities for both planets.Our eccentric fit consisted of all previously mentioned free parameters, as well as the addition of

√
cos  and √  sin  as free parameters, where  is the eccentricity and  is the longitude of periastron.These parameters could be derived to find the eccentricity and longitude of periastron of both planets.The results of both fits can be seen in Table 5.
Through visual examination of the detrended light curves against the models, assuming a constant period, we noticed that both planets deviated slightly from their predicted transit timing assuming they are both on circular orbits.We followed the algorithm outlined in Nabbie et al. (in prep) to search for potential transit timing variations (TTVs).The 2022-Apr-01 CTIO observation of TOI-5126 b and Teid observation of TOI-5126 c were excluded due to their insufficient out-of-transit baselines.For the rest of the transits, we fit for their individual transit centers independently while allowing for the trend parameters to vary similarly to the previous model.Instead of using Gaussian priors on the stellar parameters, we added the likelihood term from the SED in our model instead so that all the stellar parameters can be constrained at the same time as the transit parameters.
For the SED likelihood, we make use of the isochrone package (Morton 2015) and the MIST tracks.We fixed stellar parallax, distance and reddening to the best-fit values from 3.1.We added  eff , [Fe/H] and age as our free parameters, and an error scale parameter following Eastman et al. (2019).The priors of  eff and [Fe/H] are constrained by Gaussian priors with uncertainties provided by Section 3.1.We adopt uniform priors for  ★ and  ★ with bounds between 1 and 1.4.The error scale was bounded between 0.01 and 100, and stellar age between 0.5 and 10 Gyr.

Global modeling results
Through the results of our analysis conducted in Section 3.2, we report the discovery and validation of two planets around the bright F star TOI-5126.We have determined TOI-5126 b to be a hot super-Neptune, with a radius of 4.74 +0.16  −0.14  ⊕ in a 5.4588385 +0.0000070 −0.0000072 day orbit.Further, we have discovered and validated TOI-5126 c as a warm Neptune at 3.86 +0.17 −0.16  ⊕ .TOI-5126 c orbits its host star every 17.8999 +0.0018 −0.0013 days.Comparison with the current exoplanet demographic is found in Figure 12.The full parameter results are presented in Table 5.The eccentric fit and the circular fit have predominantly overlapped posterior space.We have also detected a tentative TTV signal which will be discussed in Section §4.3.3.As expected, the best fit parameters retrieved when allowing the transit times to vary in the global model, produced slightly larger transit depths and shorter transit durations for both planets.However, the joint modeling with the SED lead to a slightly smaller stellar radius.The posterior of the final planet parameters largely overlap, with the TTV fit reports TOI-5126 b having a radius of 4.67 +0.14  −0.13  ⊕ and a transit duration of 3.686 +0.032 −0.026 hours.For TOI-5126 c, the TTV fit reports a radius of 3.76 +0.14 −0.13  ⊕ and a transit duration of 5.46 +0.10 −0.16 hours.

Validation of the TOI-5126 system
Both astrophysical phenomena and instrumental artifacts can replicate the signal of a planetary transit.Without careful consideration of false positive scenarios, the legitimacy of the planet can be compromised.In this section, we validate the TOI-5126 system by ruling out false positive scenarios which could have created the observed transit signals.We consider the following scenarios: • Either of the transit signals are a reflection of instrumental artifact: Both transit signals were detected with high significance in the TESS data, CHEOPS and ground-based telescopes as mentioned in section 2.1.3,therefore this possibility is ruled out.
• TOI-5126 is an Eclipsing Binary: This scenario is ruled out by our radial velocity observations from TRES and CHIRON (section §2.2).An analysis on the system shows that a brown dwarf companion (  ≈ 13 Jup ) orbiting around TOI-5126 would produce RV semiamplitudes of   = 1326 m s −1 and   = 892 m s −1 .Neither sets of the radial velocity observations show such deviations.
• Light from a distant eclipsing binary or transiting planet system is blended with TOI-5126: The magnitude difference between TOI-5126 and the faintest companion which could replicate the signals observed follow as: Δ ≲ 2.5 log 10  2 12  2 13   .The faintest star which could replicate the transit signals for TOI-5126 b is Δ = 7, and TOI-5126 c is Δ = 6.Direct imaging observations rule out stars with a Δ of 7 from 0.5 ′′ onwards, and ground-based photometric observations detected the transit events within a 4.7 ′′ aperture.To statistically rule out stars within 0.5 ′′ of TOI-5126, we determined the average stellar density in the nearby field from Gaia data.For stars within a magnitude range of seven magnitudes, the average stellar density is calculated to be 2.4 × 10 −5 arcsec −2 .Narrowing this caclulation to stars six magnitudes fainter than TOI-5126, the average stellar density is 4.0 × 10 −5 arcsec −2 .
An additional investigation into the probability of remaining false positive scenarios was performed using Triceratops (Giacalone & Dressing 2020).We used the detrended TESS light curves of TOI-5126 b and TOI-5126 c, as well as constraints provided by Palomar imaging (see section 2.3) in our analysis.A false-positive probability (FPP) of 0.77% and a nearby false-positive probability (NFPP) of 0% offer statistical validation for TOI-5126 b.Analysis of the TESS light curve for planet c reveals a FPP and NFPP of 5.33% and 4.98%, respectively.This is invoked by a Nearby Eclipsing Binary scenario on 9.53 ′′ away neighbour TIC27064467.This scenario, however, has been ruled out by the on-target event detected by ground-based photometry.As a result, the FPP with this scenario eliminated is 0.35%.

Long-term stability
The masses and eccentricities of the planets are currently not wellconstrained.In this section, we explore the possibility of placing additional constraints on the system through considerations of longterm dynamical stability.We first assess the stability of the observed system using calculations of the Mean Exponential Growth factor of Nearby Orbits (MEGNO; Cincotta & Simó 2000;Cincotta et al. 2003) indicator within the REBOUND gravitational dynamics software package (Rein & Liu 2012a).For a given set of initial conditions, we quantify the stability of the system using the time-averaged MEGNO chaos indicator value, ⟨ ()⟩, which reflects the degree of divergence of initially closely-separated trajectories in phase space.
Stable system configurations will have MEGNO values ≲ 2.
The MEGNO value across a grid of different masses for both planets is displayed in the top panel of Figure 8.We are able to obtain a relatively loose constraint on the planet masses: approximately, if   + 2 3   < 3  , the system will be stable.To estimate the average MEGNO value in each cell, we randomly sample the initial condition 30 times.We assume eccentricities of both planets were randomly drawn from a Rayleigh distribution with scale parameter 0.1.The mean anomaly and argument of periapse angles are sampled uniformly between 0 and 2 .The integration time for each realization was set to 10 6 orbits of the outermost planet with a time-step of 1/20th of the innermost planet's period.
We also found the system will be stable across a wide range of eccentricities for both planets if their masses approximated by the mass-radius relationship following Otegi et al. (2020).Therefore, the stability criterion alone cannot provide additional constraint on the system eccentricity.

Stability of a hypothetical third planet
The two detected planets in the TOI-5126 system are fairly widelyspaced (  /  ∼ 3.3).It is possible that there is a third planet residing between the orbits of the two planets or somewhere else in the system.To explore this possibility, we use the Stability of Planetary Orbital Configurations Klassifier (SPOCK; Tamayo et al. 2020) to compute the probability of stability of a given set of parameters.SPOCK is a machine learning model trained on ∼ 100, 000 orbital configurations of three-planet systems.It can estimate the probability of long-term stability of planetary systems with three or more planets.
In particular, we explored the   −  space of the additional planet, whose mass we set to the average of the nominal masses of planets b and c based on the mass-radius relationship discussed in Section 4.3.1.We assumed circular orbits for planets b and c while varying the eccentricity of planet d from 0−0.3 and semi-major axis from 0−0.3 AU.Setting the orbits of the known planets to circular provides us with the maximum possible stable region for the hypothetical third planet.Each grid cell in the   −   map represents the average of 10 realizations of the system with randomized initial conditions.Specifically, the mean anomalies, longitudes of ascending nodes, and inclinations were randomly drawn according to ,  ∼ Unif [0, 2],  ∼ Rayleigh(1 • ).
The bottom panel of Figure 9 shows the results of the SPOCK stability map for a hypothetical third planet.In addition to stable regions interior to TOI-5126 b and exterior to TOI-5126 c, there is also a stable region between the two planets.It is possible that a planet may be transiting but remain undetected due to its small radius.Alternatively, the planet could be also be either non-transiting or in the trivial case not exist at all.Nevertheless, Figure 9 indicates the Bottom: Stability map generated using SPOCK for a hypothetical third planet within the TOI-5126 system.The vertical white dashed lines denote the locations of TOI-5126 b and TOI-5126 c.Each cell in the map represents the average of 10 realizations with randomized initial conditions.There is a small stable region present between the two planets in which a hypothetical third planet could reside.
ideal orbital proximity in which to search for additional planets in the system.

Transit timing variations signal in the light curves
From our analysis of the TESS light curves , we find that both planets exhibit tentative TTVs, seen in Figure 10.Such signals can potentially arise due to the gravitational perturbations that the planets impart on one another as they reach conjunction, creating non-periodic transit times.Notably, TOI-5126 b and TOI-5126 c are not near a first-order mean-motion resonance (MMR), which would produce the largest TTV amplitudes (Agol et al. 2005).
We utilised REBOUND (Rein & Liu 2012a) to estimate expected TTV signals based on TOI-5126 b and TOI-5126 c alone.We sample the initial orbital elements of both planets (, , , ) from the posterior distribution results from our global models (including the TTV fit) in Section 3. The other orbit angles are sampled uniformly between 0 and 2 .The systems are integrated for more than 500 days and transit times corresponding to the observed transits from TESS, LCO and CHEOPS are recorded.The expected TTV scatter is then computed by calculating the standard deviation of the transit times after subtracting the best fit linear ephemeris.We simulated the system for 10000 realisations and found that TOI-5126 b produced 2- upper limit variations of 0.27 minutes, with TOI-5126 c exhibiting 0.19 minute 2- upper limit TTVs.The current system is therefore unable to explain the observed TTVs and we consider alternate scenarios.
One possible interpretation is that these signals imply the existence of a third, non-transiting planet in between TOI-5126 b and TOI-5126 c (as discussed in the previous section), so that one of the neighboring planet pairs have a period ratio just wide of a 3:2 resonance and the other just wide of a 2:1 resonance, both of these period ratios being commonly found in multi-planet systems (Fabrycky et al. 2014).
However, we caution that it has been shown in previous studies that spot crossing events can mimic TTV signals rather than gravitational interactions due to planets (Sanchis-Ojeda et al. 2011;Fabrycky et al. 2012;Szabó et al. 2013;Mazeh et al. 2013), especially when the individual transit has relatively low signal-to-noise.Ioannidis et al. (2016) demonstrated that a SNR larger than 15 is required to statistically minimizes the probability that a spot crossing event being the cause of the deviations.The SNR of planet b in TESS data is 13, with planet c having a SNR of 11.Given two out of three of our CHEOPS observations does not cover either ingress or egress of the transits, further high signal-to-noise photometry follow up is needed to unlock the full system architecture of TOI-5126.

TOI-5126 b in context with the Super-Neptune Population
TOI-5126 b is the newest addition to the super-Neptune population, the class of planets with radii between 4 and 8  ⊕ , and diverse compositions (with mass ranging between 6-135  ⊕ Akeson et al. 2013).TOI-5126 b's planetary parameters (  ≈ 4.7  ⊕ ,   ≈ 5.46 days) place it in the sparsely populated ≲10 day, 4-8  ⊕ region of the period-radius plane, known as the 'hot Neptune desert' (Mazeh et al. 2016) (see Figure 12) The unusually large radius of TOI-5126 b implies that long-term mass-loss processes are likely still ongoing and may be readily observable.Of all super-Neptunes in the desert, TOI-5126 b is on an exclusive list of three (WASP-166 b, Hellier et al. 2019;LTT 9779b, Jenkins et al. 2020b) with visual magnitudes ≲10 and predicted equilibrium temperatures above the cloud-free threshold (≳1200 K).This makes it one of the most readily characterisable highly irradiated super-Neptune's for transmission spectroscopic observations, where mass-loss mechanisms can be directly probed.Questions surrounding mass-loss are particularly unique for TOI-5126 as planet b is larger than planet c, despite receiving more stellar irradiation.This is in contrast to the expected outcome of photoevaporation, which implies that the inner companion would be smaller if the insolation ratio is larger.An example of the radius and insolation ratio relation, as well as the expected outcome of photoevaporation, can be seen in Figure 11.

Other future follow up opportunities
We compute predicted masses for TOI-5126 b and c based on the mass-radius relationship from Otegi et al. (2020).These predicted masses are adopted for the discussions below.The predicted masses for TOI-5126 b and c are 21 +9 −7  ⊕ and 18 +8 −6  ⊕ , respectively.
The estimated RV semi-amplitudes for planets b and c are   = 6.6 +2.5 −2.5 m s −1 and   = 3.8 +1.5 −1.5 m s −1 , respectively.Measuring the masses of both planets are important to determine if they are formed with large protoplanetary cores and therefore a large gaseous envelope (Lee & Chiang 2016b) or inflated atmosphere due to tidal heating (Millholland et al. 2020).TOI-5126's high rotational velocity makes RV measurements difficult, but with the next generation of instruments such as MAROON-X (Seifahrt et al. 2018), along with innovative methods to correct for stellar variability, mass constraints of both planets are not ruled out.The rapid rotation of TOI-5126 does however enable follow-up studies into the spin-orbit alignment of the system with the The planets in the TOI-5126 system present an ideal opportunity to better understand the atmospheric composition difference between a hot super-Neptune and a warm Neptune.TOI-5126 b still likely retains a significant H/He atmosphere due to its inflated radius.TOI-5126 is relatively bright, with a  band magnitude of 9.04.The inner planet TOI-5126 b has an equilibrium temperature of ( eq = 1442 +46 −40 K), and should host a relatively clear atmosphere.To test its suitability for follow-up, we followed the transmission spectroscopy metric (TSM; Kempton et al. 2018) and find that TOI-5126 b has a TSM of 84 +53 −25 .Given the potential of TOI-5126 c to also be a candidate for transmission spectroscopy (TSM = 40 +27 −12 ), this system enables the potential of atmospheric comparative planetology studies between a hot super-Neptune and warm Neptune.The suitability of the TOI-5126 planets for follow-up studies can be seen in comparison to the super-Neptune population in Figure 13, where despite TOI-5126 c being a Neptune, its high insolation flux (  = 96 +9 −8  ⊕ ) and predicted equilibrium temperature ( eq = 971 +31 −27 K) make it favourable even among super-Neptunes.
To demonstrate the possibility of comparative studies, we conducted preliminary atmospheric simulations using the public software package petitRADTRANS (Mollière et al. 2019).As we are focused on identifying differences in the M/H and C/O ratios of the two planets, we used the NIRSpec Bright Object Time Series (BOTS) template with a G395H grating for our simulations.We first simulated varying M/H, including 1x, 10x and 100x Solar.We also simulate the different C/O ratios at 0.2, 0.53 (Solar) and 1.0 C/O levels, in this case assuming 10x solar M/H.The results can be seen in Figure 14.
The key features of the G395H grating are the CH 4 and CO 2 features in the spectra where both planets can be probed to compare formation theories.In this case, a high C/O ratio and metallicity would imply formation in a gas-depleted disk, while low C/O ratios for the planets would support other formation mechanisms.We show that differing CO 2 abundances can be discerned in both TOI-5126 b and TOI-5126 c in the M/H simulation, while TOI-5126 c appears to also have a small CH 4 feature at lower metallicities due to its lower predicted equilibrium temperature.The varying C/O models show that if TOI-5126 c has a high CH 4 abundance, its chemical signature will be apparent in the spectrum, while lower levels would as a result produce a larger CO 2 bump.The high equilibrium temperature of TOI-5126 b diminishes the signature of the CH 4 feature at high C/O ratios.DRC acknowledges partial support from NASA Grant 18-2XRP18_2-0007.
CHEOPS is an ESA mission in partnership with Switzerland with important contributions to the payload and the ground segment from Austria, Belgium, France, Germany, Hungary, Italy, Portugal, Spain, Sweden and the United Kingdom.We thank support from the CHEOPS GO Programme and Science Operations Centre for help in the preparation and analysis of the CHEOPS observations.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).We acknowledge the use of public TESS data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center.
Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products.
Funding for the TESS mission is provided by NASA's Science Mission Directorate.KAC acknowledges support from the TESS mission via subaward s3449 from MIT.  2020)), astroARIADNE (Vines & Jenkins (2022)) and REBOUND (Rein & Liu (2012b)), as well as their dependencies.

DATA AVAILABILITY
All RV measurements used in our joint model analysis have been included in this paper's tables.The TESS data products from which we derived our light curves are publicly available online from the Mikulski Archive for Space Telescopes (MAST).Ground-based observations used in this paper is available from ExoFOP-TESS.CHEOPS data will be shared on reasonable request to the lead author.

Figure 1 .
Figure 1.TESS SPOC SAP light curve of TOI-5126 from all observed sectors.Blue triangles mark each transit of TOI-5126 b, with red arrows for transits of TOI-5126 c. Empty triangles indicate the transits not used in global modeling.

Figure 2 .
Figure 2. Detrended TESS light curve of TOI-5126 from all observed sectors.The light curves are detrended using an iterative B-spline method, and are not used in the global modeling that lead to the final constraint of the planet parameters.Blue triangles mark each transit of TOI-5126 b, with red arrows for transits of TOI-5126 c. Empty triangles indicate the transits not used in global modeling.

Figure 3 .
Figure 3. Upper: Raw normalised CHEOPS light curves indicated by green points with the best-fit trend model plotted in purple.Middle: Detrended CHEOPS light curves shown in green.Binned points are denoted as purple markers and the best-fit transit model is overplotted as a purple line.Lower: Residuals of the detrended light curves in green, with binned averaged points shown in purple.

Figure 4 .
Figure 4. Phase-folded light curves of TOI-5126 b and TOI-5126 c. TESS SPOC light curves are noted by green points, with binned data in blue and the best-fit transit model line in blue.LCOGT data is represented by gray points with orange binned points and errorsbars in blue.The transit model is also denoted as an orange line.The CHEOPS phase folded light curve is combined with all three visits TOI-5126 b, where the purple, pink and blue points represent each visit respectively.The binned data is in cyan, with the transit model also a cyan line.(a) TESS SPOC TOI-5126 b.(b) McD TOI-5126 b.(c) CTIO 2022 TOI-5126 b.(d) CTIO 2023 TOI-5126 b.(e) CHEOPS combined visits of TOI-5126 b.(f) TESS SPOC TOI-5126 c. (g) Teid TOI-5126 c.

Figure 5 .
Figure 5. Left: Radial velocities against the orbital phase of TOI-5126 b with the notional K semi-amplitude curve based off of predicted mass in Section 4.5.TRES points are represented by purple points, while CHIRON data is in brown.Right: Radial velocities versus orbital phase for TOI-5126 c. with notional K semi-amplitude curve drawn from predicted mass.

Figure 6 .
Figure 6.Left: Palomar 5.0 m companion sensitivity curve.The black points represent the 5- limits.The inset image is of the primary target showing no additional close-in companions.Right: Contrast curve of TOI-5126 with the NESSI instrument on the WIYN 3.5 m telescope, reaching a constrast of Δ 6 mag from 0.96 ′′ .

Figure 7 .
Figure 7. Spectral energy distribution (SED) of TOI-5126.Circles indicate the photometric bands with horizontal bars reflecting the passband width.

Figure 8 .
Figure 8. Lomb-Scargle periodogram of the raw SPOC TOI-5126 light curve .The highest-power period is denoted by a dashed line.Uncertainties are derived from FWHM of the highest-power resulting in a period of 4.602 +0.071 −0.067

Figure 9 .
Figure 9. Top: MEGNO stability map in the   −   parameter space.Two distinct regions of stability are divided by the solid cyan line, with the region A below the line indicate stable configurations.Each pixel in the 20 × 20 map holds the average of 30 realizations of randomized initial conditions.Bottom: Stability map generated using SPOCK for a hypothetical third planet within the TOI-5126 system.The vertical white dashed lines denote the locations of TOI-5126 b and TOI-5126 c.Each cell in the map represents the average of 10 realizations with randomized initial conditions.There is a small stable region present between the two planets in which a hypothetical third planet could reside.

Figure 10 .
Figure 10.Transit times subtracted from a linear ephemeris for TOI-5126 b and TOI-5126 c. Cyan indicates TESS transits, yellow indicates LCOGT transits and red indicates CHEOPS transits.The red dashed line shows the expected transit time assuming a constant period.

Figure 11 .
Figure 11.Radius ratios of all known inner/outer planet pairs as a function of their insolation ratios.The colours represent the inner planets' insolation flux in Earth units.The horizontal dashed line at one represents a same sized planet pair, with the arrow outlining the decrease in radius ratio expected by photoevaporation of the inner hotter planet.This plot makes use of data from the NASA exoplanet archive (Akeson et al. 2013) downloaded on 2023-02-07 UT.The data extracted only planets with measured periods, radii and host star masses, radii and effective temperatures known.

Figure 12 .
Figure 12.Radius relative to insolation flux for the population of confirmed planets.Colours represent the host star effective temperature.Plots were made with data downloaded from NASA exoplanet archive on 2023-02-07 UT.

Figure 13 .
Figure 13.TOI-5126 planets among all super-Neptunes in terms of orbital period and Transmission Spectroscopy Metric (TSM).TOI-5126 c is notably not a super-Neptune and is outlined for comparative purposes as a dashed line.The color of each point describes the planets' insolation flux, with point size representing their scaled radius.This plot makes use of data from the NASA exoplanet archive downloaded on 2023-02-07 UT.
B.S.S. and A.A.B acknowledge the support of M.V. Lomonosov Moscow State University Program of Development.Some of the observations in the paper made use of the NN-EXPLORE Exoplanet and Stellar Speckle Imager (NESSI).NESSI was funded by the NASA Exoplanet Exploration Program and the NASA Ames Research Center.NESSI was built at the Ames Research Center by Steve B. Howell, Nic Scott, Elliott P. Horch, and Emmett Quigley.Data were reduced using a software pipeline originally written by Elliott Horch and Mark Everett.This work makes use of observations from the LCOGT network.Part of the LCOGT telescope time was granted by NOIRLab through the Mid-Scale Innovations Program (MSIP).MSIP is funded by NSF.This research has used data from the CTIO/SMARTS 1.5m telescope, which is operated as part of the SMARTS Consortium by RECONS This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.This research has made use of the Exoplanet Follow-up Observation Program (ExoFOP; DOI: 10.26134/ExoFOP5) website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

Table 1 .
List of CHEOPS Observations

Table 2 .
SG1 Follow-Up Observations a Transit observed simultaneously by two distinct telescopes at Teid. b Transit observed simultaneously by two distinct telescopes at McD.

Table 3 .
Ziegler et al. 2020)pectroscopyZiegler et al. 2020).TOI-5126 has a Gaia DR3 RUWE value of 1.16 indicating that the astrometric fits are consistent with the single star model.
a Relative velocities from TRES self cross correlation analysis RUWE values ≳ 1.4 may indicate an astrometric excess noise, possibily caused the presence of an unseen companion (e.g.,