Separated Twins or Just Siblings? A Multiplanet System around an M Dwarf Including a Cool Sub-Neptune

We report the discovery of two TESS sub-Neptunes orbiting the early M dwarf TOI-904 (TIC 261257684). Both exoplanets, TOI-904 b and c, were initially observed in TESS Sector 12 with twin sizes of 2.426−0.157+0.163 and 2.167−0.118+0.130 R ⊕, respectively. Through observations in five additional sectors in the TESS primary mission and the first and second extended missions, the orbital periods of the planets were measured to be 10.887 ± 0.001 and 83.999 ± 0.001 days, respectively. Reconnaissance radial velocity measurements (taken with EULER/CORALIE and SMARTS/CHIRON) and high-resolution speckle imaging with adaptive optics (obtained from SOAR/HRCAM and Gemini South/ZORRO) show no evidence of an eclipsing binary or a nearby companion, which, together with the low false-positive probabilities calculated with the statistical validation software TRICERATOPS, establishes the planetary nature of these candidates. The outer planet, TOI-904 c, is the longest-period M dwarf exoplanet found by TESS, with an estimated equilibrium temperature of 217 K. As the three other validated planets with comparable host stars and orbital periods were observed by Kepler around much dimmer stars (J mag > 12), TOI-904 c, orbiting a brighter star (J mag = 9.6), is the coldest M dwarf planet easily accessible for atmospheric follow-up. Future mass measurements and transmission spectroscopy of the similar-sized planets in this system could determine whether they are also similar in density and composition, suggesting a common formation pathway, or whether they have distinct origins.


Introduction
The past decade of exoplanet exploration has revealed that M dwarf stars typically host multiplanet systems of small, rocky planets with orbital periods of P < 20 days (Dressing & Charbonneau 2015;Hardegree-Ullman et al. 2019).As these cool, low-mass stars allow for easier detections and characterizations of exoplanets than larger hosts, whether via deeper transits or more significant radial velocity reflex motion, these short-period planets have been prioritized for in-depth study by the exoplanet community.Yet few planets have been found at greater distances from M dwarf hosts to date.Indeed, the Kepler mission only discovered ∼30 M dwarf planets at distances of >0.15 au (periods P > 25 days; Muirhead et al. 2012; Barclay et al. 2015;Dressing & Charbonneau 2015;Torres et al. 2015Torres et al. , 2017;;Morton et al. 2016;Berger et al. 2018), all of which orbit stars dimmer than V mag = 15 and thus are not easily accessible for ground-based observations.By prioritizing M dwarf stars and focusing on stars in the solar neighborhood, the Transiting Exoplanet Survey Satellite (TESS) now has the potential to populate this underexplored region of parameter space and provide targets for future characterization.
TESS is conducting a full-sky survey of nearby, bright stars, 75% of which are M dwarfs (Ricker et al. 2015).Unlike Kepler, which observed one region of the sky for 4 yr, the TESS primary mission (PM) observed both hemispheres by cycling through 26 sectors observed for ∼27 day intervals, and it is continuing to reobserve the hemispheres (as well as the ecliptic) in the same fashion through its extended missions.While this observation strategy makes TESS most sensitive to planets with P  10 days, there exists the potential to observe planets with longer orbital periods in overlapping sectors near the ecliptic poles.Further, TESS is projected to observe many planets with P  25 days as single-transit events (Villanueva et al. 2019) that could be recovered in later extended missions.Since its launch in 2018, TESS has found 28 confirmed planets with periods of >25 days, five of which (Cañas et al. 2020;Rodriguez et al. 2020;Fukui et al. 2022;Mann et al. 2022;Schanche et al. 2022) orbit low-mass stars.
In this Letter, we introduce and statistically validate TESS object of interest (TOI) 904 c, a sub-Neptune with a period of 84 days that was observed to produce a single-transit event in four different sectors of the TESS PM and its first and second extended missions (EM1 and EM2).We also introduce another planet in the same system (TOI-904 b) with a period of 10.87 days.Orbiting an early M dwarf (TIC 261257684, TOI-904) at T equil ≈ 200 K, TOI-904 c is the coldest M dwarf planet discovered by TESS to date.Only three other transiting planets have been observed around low-mass stars with similar orbital periods: Kepler-186 f (P = 129.94 days; Quintana et al. 2014), Kepler-1229 b (P = 86.83days; Morton et al. 2016), and Kepler-1628 b (P = 76.38 days; Morton et al. 2016).Unlike these planets, TOI-904 c orbits an M dwarf bright enough (J mag = 9.61) for future mass measurements and atmospheric characterization.Through these additional observations, this system embodies a (currently) unique opportunity to test theories of planet formation around low-mass stars.Further, by measuring the densities and constraining the composition of both planets in this system, we can determine whether these similarly sized planets also have twin compositions and formation histories or if they evolved through distinct formation pathways while orbiting the same host star.

TESS Observations
TESS observed TOI-904 in Sectors 12 and 13 (2019 May 21-2019 July 17) of its PM; Sectors 27, 38, and 39 (2020 July 5-2020 July 30 and 2021 April 29-2021 June 24) of its EM1; and Sector 61 (2023 January 18-2023 February 12) of its EM2.TOI-904 was included in the TESS Candidate Target List (Stassun et al. 2019) and monitored in both the 2 minute postage stamp and longer-cadence (30 minutes in PM, 10 minutes in EM1, 200 s in EM2) full-frame images during both missions (see Figure 1).The transit signatures of TOI-904 b and c were initially detected by the Science Processing Operations Center (SPOC; Jenkins et al. 2016), located at the NASA Ames Research Center, in a transit search of Sector 12 with a noise-compensating matched filter (Jenkins 2002;Jenkins et al. 2010Jenkins et al. , 2020) ) on 2019 July 1.This filter originally folded a single transit of TOI-904 c onto the second transit of TOI-904 b.Both the SPOC and the Visual Survey Group (Kristiansen et al. 2022), working in conjunction with the TESS Single Transit Planet Candidate (TSTPC) working group, flagged these transit events and attributed both to a single-planet candidate with an orbital period of 18.35 ± 0.005 days.The correct period of TOI-904 b, 10.877 ± 0.001 days, was identified by the SPOC via a transit search of Sector 13 conducted on 2019 July 27.The transit signature passed all of the diagnostic tests presented in the resulting data validation report (Twicken et al. 2018) and was fitted with an initial limbdarkened transit model (Li et al. 2019).The difference image centroiding test located the source of the transit signature to within 1 3 ± 2 9.The TESS Science Office reviewed the diagnostic results and issued an alert for TOI-904.01 on 2019 June 23 (Guerrero et al. 2021), and it was then recognized as a TOI on 2019 July 15 on the TESS data alerts web portal at the Massachusetts Institute of Technology. 20he SPOC conducted a subsequent multisector search of Sectors 12, 13, 27, 38, and 39 on 2021 July 25, where they recovered 12 transits of TOI-904 b and the signature of TOI-904 c at 4× the true period (see Appendix B).The SPOC found that the difference imaging centroiding test located the source of the transit signature to within 3 4 ± 3 4.In 2021 August, a member of the TSTPC working group (Hugh Osborn) found an undetected additional transit of TOI-904 c in SPOC-produced light curves near the end of TESS Sector 39, which, together with the Sector 12 transit, constrained the period of this planet to a discrete set of possible aliases.We proceeded to search all previous TESS observations of this star for additional transits of this outer planet and discovered a third clear transit in TESS Sector 27 from which we were able to derive a unique period of 83.999 ± 0.001 days.This planet was registered as a community TOI on 2021 September 1, and the TESS Science Office issued an alert for it on 2022 April 20.We searched the TESS observations of TOI-904 for any additional transit signals and found no evidence of additional planets in this system from the transit method.
To investigate possible false-positive scenarios of the planets that were observed in this system, we visually inspected the background of each individual light curve that included a transit of the outer planet to rule out false positives due to asteroids or other anomalies.We also used the Gaia DR3 (Gaia Collaboration et al. 2023) to investigate the possibility of blended objects in the aperture and whether the renormalized unit weight error (RUWE), a measure of the normalized χ 2 of the Gaia observations to the astronometric single-star fit corrected for color and magnitude dependencies (Lindegren et al. 2018;Pearce et al. 2020), indicated that the host star existed in a binary that could be responsible for the observed transits.With a RUWE metric of ≈0.999, TOI-904 showed no obvious indication of being a stellar binary (Pearce et al. 2020).As our preliminary inspection gave no indication that either planet candidate's transits were false positives, we looked to follow-up observations (as described in the following sections) to further validate this system.

Ground-based Photometry
The TESS pixel scale is ∼21″ pixel −1 , and photometric apertures typically extend out to roughly 1′, which generally results in multiple stars blending in the TESS aperture.An eclipsing binary in one of the nearby blended stars could thus mimic a transit-like event in the large TESS aperture.With additional ground-based photometric observations, we can attempt to (1) rule out or identify nearby eclipsing binaries (NEBs) as potential sources of the detection in the TESS data, (2) check for the transit-like event on-target using smaller photometric apertures than in the TESS images to confirm that the event is occurring on-target or in a star so close to TOI-904 that it was not detected by Gaia DR3 (which is unlikely, as such a star would be too faint unless it is perfectly aligned with the target star, but which is accounted for in the triceratops transit probability calculation in Section 3.2), and (3) refine the TESS ephemeris.

Las Cumbres Observatory
We acquired ground-based transit follow-up photometry of TOI-904 b as part of the TESS Follow-up Observing Program Sub Group 1 (Collins 2019). 21We used the TESS Transit We used circular photometric apertures of radius 1 9 and 5 8 to check for possible NEBs that could be contaminating the SPOC photometric apertures, which generally extend ~¢ 1 from the target star.To account for possible contamination from the wings of neighboring star point-spread functions, we searched for NEBs in all known Gaia EDR3 and TESS Input Catalog (TIC) version 8 nearby stars out to ¢ 2. 5 from TOI-904 that are possibly bright enough in the TESS band to produce the TESS detection (assuming a 100% eclipse and 100% contamination of the TESS aperture).To attempt to account for possible delta-magnitude differences between the TESS band and the follow-up Pan-STARRS z-short band, we checked stars that are an extra 0.5 mag fainter in the TESS band than needed.We consider a star cleared of an NEB if the rms of its 10 minute binned light curve is more than a factor of 5 smaller than the adjusted expected NEB depth in the star (adjusted to allow for the potential TESS-band delta-magnitude difference).
We then visually inspect each neighboring star's light curve to ensure that there is no obvious eclipse-like signal.The NEBchecked light-curve data are available at ExoFOP-TESS. 22We rule out an NEB as the source of the TOI-904 b signal in the TESS data.Photometry of TOI-904 was extracted using circular apertures with a radius of 5 8, which exclude flux from the nearest known Gaia EDR3 and TIC neighbor (TIC 724109412), 25 4 southwest.We detect the transit event within the TOI-904 photometric apertures in the two z-shortband light curves and include the data in the analyses of this work (see Figure 2).
We attempted to observe an additional transit of TOI-904 c on UTC 2022 February 27 with the South Africa Astronomical Observatory, resulting in a tentative partial transit detection.This detection, however, lacked a sufficient out-of-transit baseline to be conclusive and is not included in the analyses of this work.

Spectroscopic Follow-up
Brown dwarf or grazing stellar binaries are frequent sources of false-positive transit signals.While few instruments are capable of confirming the exoplanet nature of transits by measuring the masses of small extrasolar planets, many more are equipped to detect the Doppler signal caused by a stellar or brown dwarf mass companion.Reconnaissance radial velocity measurements of the host star are thus able to rule out a false positive due to a bound stellar (or brown dwarf) binary while also providing key observations to refine the spectroscopic parameters of the host star.

SMARTS/CHIRON Spectroscopy
We obtained four observations of TOI-904 with the CHIRON facility on the SMARTS 1.5 m telescope at Cerro Tololo Inter-American Observatory, Chile (Tokovinin et al. 2013).CHIRON is a high-resolution spectrograph with a resolving power of R = 80,000 over a wavelength range of 4100-8700 Å with "slicer" mode observations.The spectra were extracted via the standard pipeline as in Paredes et al. (2021).
We attempted to derive a line broadening velocity from the CHIRON observations.We performed a least-squares deconvolution between each spectrum and a nonrotating synthetic spectrum generated via the ATLAS9 model atmospheres grid .Observations were taken in 36 s intervals (light blue) and binned to 5 minutes (dark blue).The best-fit transit model is overlaid in black, showing that both observations recover full transits of the planet at a comparable depth to the TESS observations.(Castelli & Kurucz 2004) at the stellar effective temperature and surface gravity of our target star.We then modeled the line broadening profile via a combination of the rotational, radialtangential macroturbulent, and instrumental broadening kernels (as in Gray & Corbally 1994).We find the rotational broadening component to have a width of <2 km s −1 (1σ) when the additional broadening terms are considered, consistent with expectations from the photometric modulation-derived rotation period of the target star.

EULER/CORALIE Spectroscopy
TOI-904 was observed by the Swiss 1.2 m Euler telescope with the CORALIE instrument installed at its Nasmyth focus.CORALIE is a fiber-fed high-resolution spectrograph with a spectral resolution of 60,000 (Queloz et al. 2001).CORALIE has a 3 pixel sampling per resolution element.Five spectra of TOI-904 were taken between 2020 February 2 and April 14 (see Figure 2).The observations have an exposure time of 2400 s and reach a signal-to-noise ratio varying between 11 and 16.The observations are targeting the extreme phases of the estimated radial velocity signal (see Figures 2(a) and (b)).The radial velocity is extracted by cross-correlating an M2 stellar mask with each spectrum (Pepe et al. 2002).The radial velocity data have an rms of 15 m s −1 over the time span of the observations.

High Angular Resolution Imaging
Close stellar companions (bound or line of sight) can confound exoplanet discoveries in a number of ways.The detected transit signal might be a false positive due to a background eclipsing binary, and even real planet discoveries will yield incorrect stellar and exoplanet parameters if a close companion exists and is unaccounted for (Furlan & Howell 2017, 2020).Additionally, the presence of a bound companion star leads to the nondetection of small planets residing within the same exoplanetary system (Lester et al. 2021).

SOAR/HRCAM High-resolution Imaging
We searched for stellar companions to TOI-904 with speckle imaging with the 4.1 m Southern Astrophysical Research (SOAR) telescope (Tokovinin 2018) on 2022 January 7 UT, observing in the Cousins I band, a similar visible bandpass as TESS.This observation was sensitive to a 5.4 mag fainter star at an angular distance of 1″ from the target.No nearby stars were detected within 3″ of TOI-904 in the SOAR observations.2.4.2.Gemini South/ZORRO High-resolution Imaging TOI-904 was observed on 2021 October 21 UT and 2022 January 13 UT using the ZORRO speckle instrument on the Gemini South 8 m telescope (Scott et al. 2021;Howell & Furlan 2022).ZORRO provides simultaneous speckle imaging in two bands (562 and 832 nm) with output data products including a reconstructed image with robust contrast limits on companion detections.While both observations had consistent results that TOI-904 is a single star to within the angular and contrast levels achieved, the 2022 January observation had better seeing, which led to deeper contrast levels.Seven sets of 1000 × 0.06 s images were obtained and processed in our standard reduction pipeline (Howell et al. 2011).Figure 2(c) shows our final contrast curves and the 832 nm reconstructed speckle image.We find that TOI-904 is a single star with no companion brighter than 5-8 mag below that of the target star from the 8 m telescope diffraction limit (20 mas) out to 1 2. At the distance of TOI-904 (d = 46 pc), these angular limits correspond to spatial limits of 0.9-55 au.

Host Star Parameters
As an independent determination of the basic stellar parameters from those in the TIC (Stassun et al. 2019), we performed an analysis of the broadband spectral energy distribution (SED) of the star together with the Gaia EDR3 (Gaia Collaboration et al. 2016, 2018) parallax (with no systematic offset applied; see, e.g., Stassun & Torres 2021) in order to determine an empirical measurement of the stellar radius, following the procedures described in Stassun & Torres (2016) and Stassun et al. (2017Stassun et al. ( , 2018)).We pulled the JHK S magnitudes from 2MASS (Skrutskie et al. 2006), the W1-W4 magnitudes from WISE (Wright et al. 2010), the G BP G RP magnitudes from Gaia (Riello et al. 2021), and the near-UV (NUV) flux from GALEX (Martin et al. 2003(Martin et al. , 2005)).Together, the available photometry spans the stellar SED over the wavelength range 0.2-22 μm.
We performed a fit using NextGen stellar atmosphere models, with the effective temperature (T eff ) and metallicity ([Fe/H]) adopted from the spectroscopic analysis using the SpecMatch-Emp software tool (Yee et al. 2017) on the CORALIE spectra (the surface gravity, g log , has very little influence on the broadband SED).We limited the extinction, A V , to the full line-of-sight value from the Galactic dust maps of Schlegel et al. (1998).The resulting fit has a reduced χ 2 of 2.3 (not including the NUV flux, which suggests the presence of chromospheric activity) with best-fit A V = 0.03 ± 0.03.Integrating the model SED gives the bolometric flux at Earth, F bol = 7.62 ± 0.45 × 10 −10 erg s −1 cm −2 .Taking the F bol together with the Gaia parallax gives the stellar radius, R å = 0.527 ± 0.021 R e .The stellar mass can also be estimated via the empirical M K -based relations of Mann et al. (2019), giving M å = 0.557 ± 0.028 M e .
Finally, we can use the star's rotation period, P rot , to estimate its age via empirical gyrochronology relations.We use the full TESS light curve and a Lomb-Scargle periodogram as implemented in astropy (Press & Rybicki 1989) to obtain a rotation period of 8.55 days with a false-alarm probability on the order of 10 −181 .However, similar spot coverage on opposite hemispheres can cause aliasing at half the true rotation period, meaning that the real rotation period could be 17.1 days.If, in fact, the period is 8.55 days, using empirical relations for M dwarfs from Engle & Guinan (2018) suggests that the stellar age is approximately 1 Gyr.This relatively young age would be consistent with the chromospheric activity suggested by the NUV flux in the SED; however, our current constraints on v i sin from radial velocity observations of this target favor the 17.1 day period.Precise characterization of the star's v i sin and other activity indicators, such as log ¢ R HK , is needed to better constrain the age of the host star.
These and additional stellar parameters are shown in Table 1.

Planet Validation
In this section, we discuss and rule out false-positive scenarios that could explain the transit signals we have detected in our observations of TOI-904.Our radial velocity observations taken with the CORALIE instrument sampled the expected times of radial velocity maxima and minima for both planets and detected no signal indicating a stellar or brown dwarf mass companion to the host star.Using the radvel package (Fulton et al. 2018) under the assumption of circular orbits, we obtain 3σ mass limits of 238.85 and 350.34 M ⊕ for TOI-904 b and c, respectively, which precludes a stellar mass object at these orbital periods.The high angular resolution images of TOI-904 taken with SOAR and ZORRO show no nearby stars brighter than 5-8 mag less than the host star to distances as close as 0.9 au, from which we can rule out the possibilities that the transits are located on a nearby star or that a nearby star is diluting our observations.Together, the reconnaissance radial velocities and imaging observations lead us to rule out the possibility that a blended eclipsing binary could be causing either of the observed transiting planet candidates.
Our ground-based photometric observations of TOI-904 b further allow us to eliminate the possibility that the TESS observations of TOI-904 b occurred on another star in the TESS pixel.Both observed transits had a strong significance of detection of 11σ (see Figure 2(d)).Although we have yet to recover a complete transit of TOI-904 c from ground-based instruments, statistical analyses from the Kepler survey (Lissauer et al. 2012;Rowe et al. 2014;Morton et al. 2016) and recent calculations based on TESS observations from Guerrero et al. (2021) indicate that multiplanet transiting systems are nearly always true planets, especially those with planets smaller that 6 R ⊕ , adding validity to the planetary nature of TOI-904 c.
Finally, we used the triceratops (Giacalone & Dressing 2020) software library to statistically validate both planet candidates.triceratops begins by using the Mikulski Archive for Space Telescope (MAST) module of astroquery (Ginsburg et al. 2019) to obtain the TIC properties for each star within 10 pixels of the target star.Using the TESS magnitudes, each neighboring star is considered for its flux contribution in the target pixel and as a possible source of the transit signal.The tool then creates models of transiting planet and eclipsing binary light curves that are used to calculate the probability of each scenario using a Bayesian framework.These probabilities are then used to determine whether or not a planet candidate can be classified as "validated" (has a falsepositive probability, FPP, of <1.5% and nearby FPP of <0.1%).Previously, in their work statistically validating several hundred TOIs, Giacalone & Dressing (2020) found the FPP (defined as the summed probability of all scenarios in which a planet-sized object is located on the given host star subtracted from unity) for TOI-904 b to be 3% and designated the planet as a "likely planet."By doing a joint fit using the new Gemini South/ZORRO observations discussed in Section 2.4 and incorporating all TESS transit detections of both planets, we find that for TOI-904 b, the 3σ upper limit on the FPP = 0.049%.Doing the same calculation for TOI-904 c, we calculate the 3σ upper limit to be 0.0011%, allowing us to statistically validate both planets with >99% certainty.With the added validity from the planet multiplicity of this system (Lissauer et al. 2012), we can confidently state that both objects observed around this star are planets.

Fitting Planet Parameters
We conducted our analysis on the raw 2 minute simple aperture photometry light curves produced by the SPOC pipeline (Twicken et al. 2010;Morris et al. 2020).We obtained all available data using the MAST portal.We removed longtimescale correlated noise (caused by stellar variability) by applying the wotan software library's Tukey biweight algorithm (Hippke et al. 2019) to each sector of TESS observations after masking known transit signals.To conduct a transit analysis of both signals in each light curve, we use the juliet software library (Espinoza et al. 2019) built on the batman (Kreidberg 2015) transit modeling software and the dynesty (Speagle 2020) nested sampling algorithm for calculating Bayesian posteriors and evidences.We used the LDTK (Husser et al. 2013;Parviainen & Aigrain 2015) software to calculate the limb-darkening parameters in the TESS passband based on the stellar properties given in Table 1.We assumed linear ephemerides, circular orbits, and quadratic limb-darkening while fitting the following parameters: R p /R * , b, T 0 , P, and stellar density ρ * .We performed both a joint multiplanet fit and individual fits for each planet (in which we masked out all transits of the other planet).We also did an independent check of the juliet individual planet fit analyses using the batman modeling package with the emcee (Foreman-Mackey et al. 2021) python package.The emcee software was used to perform an invariant Markov Chain Monte Carlo sampling, initializing 100 walkers and having each fit take 10,000 steps with 8000 burn-in steps.Using the same limb-darkening values, we fitted the parameters T 0 , P, R p /R * , i cos , and a/R * .Our single-planet analysis was consistent with the results from the juliet analysis (see Appendix A, Table 2).
As an additional cross-check, we repeated this analysis for light curves produced by different pipelines to check whether the transit signals were unaffected by different methods of aperture selection and background subtraction.We repeated our analysis for the 30 minute light curves created by the TESS-SPOC, QLP, and eleanor (Feinstein et al. 2019) pipelines (see Appendix A).As with the 2 minute observations, we used the raw un-detrended data before repeating the analysis detailed above using the emcee and juliet single and multiplanet fits.We found that the transit fits created from the eleanor, QLP, and 2 and 30 minute SPOC data were in agreement within 1σ.
Note.The fits created from the eleanor light curves generally show a deeper transit depth for the outer planet, TOI-904 c.The QLP light curves presented shallower transits but within an uncertainty of the SPOC and TESS-SPOC observations, as did the fits made from eleanor observations of the inner planet.The LCO observations show the inner planet to be both smaller than expected and very comparable with the TESS-found sizes of the outer planet.Ground-based observations of TOI-904 c are needed to further elucidate the difference in the two planets' radii.
We fit the LCO observations described in Section 2.2.1 of TOI-904 b using the juliet software library and our own emcee fit as described above and using the LDTK package to calculate the star's limb-darkening parameters in the SLOAN/ SDSS z passband.We found that the transit depths for TESS were within 1σ of the fit we performed on the 2 minute SPOC data, as shown in Appendix A, Table 2.

Discussion
In this Letter, we have announced the discovery and validation of TOI-904 c, a cool sub-Neptune (T eq ≈ 217.26 ± 10.22 K; R p = -+ Å R 2.167 0.118 0.130 ; P = 83.999± 0.001 days) that is the coldest planet orbiting an M dwarf star discovered by TESS to date.We also report on TOI-904 b, another sub-Neptune (R p = -+ Å R 2.426 0.157 0.163 ; P = 10.8772 ± 0.0003 days) located much closer to the early M dwarf host star.
M dwarfs are known to host a number of Earth-to Neptunesized planets (see Figure 3).These planetary systems are generally very compact and tend to exist within ∼0.2 au of their host star.By contrast, TOI-904 is an extended system with ∼0.23 au between the two planets.While distant planets have been found around other M dwarfs, all of these planets were found either on the periphery of a compact system comprised of several other small planets (e.g., Kepler-186) or as the lone planet in their systems (e.g., Kepler-1628 and Kepler-1229).As this is already a multiplanet system, we searched for indications of any other planets orbiting TOI-904 to determine if this system belongs in the former category.One way in which we looked for nontransiting planets was by searching for transittiming variations (TTVs) in the TESS and LCO observations of TOI-904 b and the TESS observations of TOI-904 c.We found the maximum TTV amplitude of a sine wave fit with periods from 0.1 to 100 days to be 5-10 minutes for the inner planet and 25-30 minutes for the outer planet (see transit times presented in Appendix B, Table 3).We found no current evidence of a nontransiting planet in the TTVs of the observed transits for TOI-904 b or c.
We also used the dynamite software library created by Dietrich & Apai (2020) to statistically calculate the potential periods and probabilities of unseen planets in multiplanet systems.dynamite calculates these probabilities by implementing a triple integral over the probability density function (PDF) of planet inclination, period, and radius based on the occurrence rates calculated by Mulders et al. (2018), assuming each variable is independent (Dietrich & Apai 2020).This PDF is then sampled using the Monte Carlo method before producing all dynamically stable results based on calculations using Dietrich & Apai's (2020) Equation (6) (see Figure 4).Our implementation of dynamite predicted four potential planets in the TOI-904 system with orbital periods of 6.74, 18.5, 50.7, and 139.05 days (see Figure 4).These additional planets may be detectable using extreme-precision radial velocity observations or potential TTVs in future observations of TOI-904 b and c, both of which will be reobserved by TESS in Sectors 65, 66, and 67.
Especially in the absence of another planet, the planets orbiting TOI-904 represent an interesting case study for exoplanet formation scenarios around low-mass stars.TOI-904 c is one of three known transiting M dwarf planets that orbit beyond 0.2 au and have R p > 1.8 R Earth , implying that they could have a gaseous envelope.The other two (Morton et al. 2016;Berger et al. 2018), however, are Kepler systems with J mag > 12, while TOI-904 has J mag = 9.6.This system thus provides an unprecedented opportunity to constrain planet formation at larger distances from M dwarfs by probing its atmosphere.Assuming a predicted mass of 6.1 M ⊕ (based on the mass-radius relations of Chen & Kipping 2017), the transmission spectroscopy metric (TSM; Kempton et al. 2018) of TOI-904 c is 23.Thanks mainly to its higher equilibrium temperature, the TSM of TOI-904 b (assuming a mass of 6.9 M ⊕ ) is 50.By comparison, the TSM values of Kepler-1628 b and Kepler-1229 b are 13 and 4, respectively.
Additional study of this system could resolve the currently ambiguous composition of the two planets.Simulations conducted by Burn et al. (2021) and Pan et al. (2022) of planet formation via core accretion around low-mass stars predict that planets with radii of ∼2.3 R ⊕ could have a range of densities spanning the ultradense sub-Neptune, water-world, and puffy sub-Neptune paradigms described by Luque & Pallé (2022).Both studies expect that in any case, planets of this size are likely to have some atmosphere, so it is not realistic to expect these planets to be rocky in composition.The mass limits we obtained using radvel on our CORALIE measurements rule out none of these scenarios, with 2σ upper limits of 174 and 244 M ⊕ on the mass of TOI-904 b and c, respectively.With follow-up mass and atmospheric measurements, we can resolve this degeneracy for the outer planet and gain insight into formation of cold planets around low-mass stars.With the predicted masses mentioned above, we expect radial velocity semiamplitudes of 2.9 and 1.3 m s −1 for TOI-904 b and c, respectively.Based on Luque & Palléʼs (2022) analysis of M dwarf small planets, the three different regimes of planet density differ in scale by a factor of 2 (0.25, 0.5, and 1 ρ ⊕ for puffy, water, and rocky planets, respectively).The masses for planets in each regime would similarly differ, resulting in detectable differences in the expected semiamplitudes of radial velocity measurements.We find that the three compositions would result in masses (semiamplitudes) of 3.3 M ⊕ (1.4 m s −1 ), 6.5 M ⊕ (2.9 m s −1 ), and 13.1 M ⊕ (5.8 m s −1 ), respectively, for the inner planet; for the outer planet, the masses (semiamplitudes) would be 2.6 M ⊕ (1.1 m s −1 ), 5.2 M ⊕ (2.3 m s −1 ), and 10.4 M ⊕ (4.6 m s −1 ), respectively.While a 5σ mass measurement (Batalha et al. 2019) could be achieved with HARPS (Mayor et al. 2003) or Magellan II/PFS (Crane et al. 2010;Teske et al. 2016), particularly for TOI-904 b, the mass precision needed to distinguish between different bulk compositions is realistically only within reach of VLT-ESPRESSO (Pepe et al. 2021).Indeed, ESPRESSO was used to measure a semiamplitude of 2.2 m s −1 for the candidate planet LHS 1140d, which has a period of 79 days and orbits a star five times fainter than TOI-904 (Lillo-Box et al. 2020).
If TOI-904 c is an ultradense Neptune and thus most easily detectable for mass measurements, both planets may have formed in situ (Kennedy et al. 2006;Hansen & Murray 2012;Hansen 2015).This case would have resulted in TOI-904 c forming as a rocky, ultradense Neptune (similar to K2-110 b; Osborn et al. 2017) if the snow line of the M dwarf receded toward the star on planet formation timescales causing a late stage of mass accretion for planets located near the snow line (Kennedy et al. 2006).In situ formation seems unlikely for this system, however, as the larger inner planet, TOI-904 b, could not have reached its size at its current location due to the high The yellow circles indicate the periods of the two known planets in the system, while the dark blue lines show the PDF calculated for this system, and the light blue histograms represent the stable Monte Carlo iterations sampled from the PDF.These scenarios were calculated using the dynamite Exoplanets Systems Simulator model (syssim; He et al. 2019).By incorporating the location of known planets and the stellar type, the dynamite software finds that additional nontransiting planets are most likely to have periods of 6.74, 18.5, 50.7, and 139.05 days.temperature (>2000 K; Ali-Dib et al. 2020) of a forming protoplanet at this distance from the host star, precluding the accretion of any gaseous envelope.In situ formation also often predicts 4-10 rocky planets on coplanar orbits (Hansen & Murray 2012;Pan et al. 2022), for which there is currently no evidence in this system.
It is more probable that TOI-904 c migrated to its current location from a greater distance via type I migration (Burn et al. 2021;Luque & Pallé 2022;Pan et al. 2022).If the planet did form further in the disk, its core was likely formed of both rock and ice materials before the planet migrated inward.In this case, this planet's density could allow us to differentiate between a "water world" with ∼50% ice/rock ratios and little to no envelope (Burn et al. 2021;Luque & Pallé 2022) and a planet with a lower ice/rock ratio and a larger atmosphere (Pan et al. 2022).In this model, the inner planet's core could be made up of rocky materials or an ice/rock mixture.If the latter, Bitsch et al. (2019) stated that the hot inner water world would be considerable evidence for the planet migration model, and TOI-904 b's overall water content could inform the exact migration history of that planet.If, instead, planet b has a large gaseous envelope, its water content could still inform which formation path this system followed and where it acquired its envelope, either beyond the snow line (Burn et al. 2021) or after being trapped in the inner regions of the protoplanetary disk, where little water is available to accrete (Pan et al. 2022).
Type I migration is a very efficient process, however, and models of this theory predict very compact multiplanet systems, unlike the widely separated planets orbiting TOI-904.One possibility to resolve this discrepancy is that TOI-904 c formed at a significantly greater distance from the host star than TOI-904 b, which would explain why the inner planet migrated in so much further and why the mutual separation of the planets has persisted.This could also result from a lower disk density at greater distances, where type I migration could be triggered for smaller planets like TOI-904 c (Burn et al. 2021).This theory could be tested by determining if TOI-904 c has a larger C/O ratio than TOI-904 b, which would indicate that TOI-904 c may have formed beyond the methane ice line and that TOI-904 b did not.
Given that the planets are so similar in size, a final alternate theory could be that they followed the same formation path by forming sequentially, following the inside-out formation theory suggested by Chatterjee & Tan (2014).In this theory, pebbles form in outer regions of this disk before migrating inward and stalling at a high-pressure region (potentially at the water-ice line) before forming a planet at that location.The planet would then migrate inward and cause the high-pressure zone to move outward, allowing the process to repeat and form another planet.Pan et al. (2022) conducted the first simulations of this type of formation around M dwarfs and found that it is likely to result in systems of two to four sub-Neptune-sized planets located at greater distances and separations than type I migration would produce and thus more closely resembling the TOI-904 system.If both planets have a high water content, TOI-904 b's larger size could be explained by an inflated hydrosphere caused by its closer proximity to the star (Luque et al. 2019).If the compositions of these planets are similar, the TOI-904 system could strongly support this relatively new theory of planet formation.Whatever theory of formation dominates in these systems, through investigating the mass and atmospheric composition of both planets, we can determine whether they have twin compositions or are only sibling planets orbiting the same star.

Finder,
which is a customized version of the Tapir software package (Jensen 2013), to schedule our transit observations and AstroImageJ (Collins et al. 2017) to extract differential photometry.We observed the predicted transit windows of TOI-904 b in the Pan-STARRS z-short band from the Las Cumbres Observatory Global Telescope (LCOGT; Brown et al. 2013) 1.0 m network nodes at Siding Spring Observatory and South Africa Astronomical Observatory on UTC 2020 October 15 and 2020 December 19, respectively.The 1 m telescopes are equipped with 4096 × 4096 SINISTRO cameras having an image scale of 0 389 pixel -1 , resulting in a ¢ ´¢ 26 26 field of view.The images were calibrated by the standard LCOGT BANZAI pipeline (McCully et al. 2018).

Figure 1 .
Figure 1.TESS observations of TOI-904.Panel (a): raw observations of each TESS sector, with 2 minute data shown in light gray and 10 minute binned data shown in dark gray.We show the Tukey biweight trend calculated with the wotan (Hippke et al. 2019) library of the variability in dark blue.Each planet transit is denoted by an arrow at the bottom of the graph, with blue denoting transits of TOI-904 b and gold denoting transits of TOI-904 c. Panel (b): TESS data detrended by the same fit, with each transit highlighted in the same color as the arrows in the previous figure.Panel (c): phase-folded light curves of transits of TOI-904 b, with the transit fit including a shaded region denoting 1σ uncertainty.Panel (d): phase-folded light curves of TOI-904 c, with the transit fit including a shaded region denoting 1σ uncertainty.

Figure 2 .
Figure 2. Follow-up observations of the planetary system around TOI-904.Left: CORALIE observations phase-folded on the periods of both TOI-904 b (panel (a), dark blue) and TOI-904 c (panel (b), gold).The best-fit radial velocity trends for both phase-folded data sets are shown in light blue.Neither figure shows evidence of a stellar or brown dwarf mass companion.Panel (c): high-resolution imaging taken at 562 and 832 nm (blue and gold, respectively).There is no clear flux from another object between 0 1 and 1 2 down to 4 and 8 mag differences, respectively, and thus no evidence of a blending or binary star.Panel (d): ground-based photometric observations of TOI-904 b taken with the LCOGT Siding Spring Observatory (panel (d), left) and the South Africa Astronomical Observatory (panel (d), right).Observations were taken in 36 s intervals (light blue) and binned to 5 minutes (dark blue).The best-fit transit model is overlaid in black, showing that both observations recover full transits of the planet at a comparable depth to the TESS observations.

Figure 3 .
Figure3.Plots of all known transiting multiplanet systems (left) and cold planets (right; T eq < 300 K) around M dwarf stars.Each system is shown as a function of the system name (left y-axis), stellar mass (right y-axis), semimajor axis (x-axis) planet equilibrium temperature (color), and planet radius (point size).TOI-904 is shown at the bottom of each panel with a gold horizontal line.

Figure 4 .
Figure 4. Created using dynamite(Dietrich & Apai 2020), the figure shows the relative likelihood of the most probable locations of unseen planets orbiting TOI-904 in log-period space.The yellow circles indicate the periods of the two known planets in the system, while the dark blue lines show the PDF calculated for this system, and the light blue histograms represent the stable Monte Carlo iterations sampled from the PDF.These scenarios were calculated using the dynamite Exoplanets Systems Simulator model (syssim;He et al. 2019).By incorporating the location of known planets and the stellar type, the dynamite software finds that additional nontransiting planets are most likely to have periods of 6.74, 18.5, 50.7, and 139.05 days.

Figure 6 .
Figure 6.Same as Figure 5 but fitting for TOI-904 c.

Table 2
Parameters Derived from the juliet, emcee, and Joint juliet Planet Fits for Each Pipeline's Data Products Used in This Study and the LCO Observations of TOI-904 b

Table 3
The Transit Times for All TESS Transits of TOI-904 b and c as They Were Observed inSectors 12, 13, 27, 38, 39, and 61Two transits of TOI-904 b were excluded due to poor quality data near 2333.999and 2975.741. Note: