TOI-712: a system of adolescent mini-Neptunes extending to the habitable zone

As an all-sky survey, NASA's $TESS$ mission is able to detect the brightest and rarest types of transiting planetary systems, including young planets that enable study of the evolutionary processes that occur within the first billion years. Here, we report the discovery of a young, multi-planet system orbiting the bright K4.5V star, TOI-712 ($V = 10.838$, $M_\star = 0.733_{-0.025}^{+0.026} M_\odot$, $R_\star = 0.674\pm0.016 R_\odot$, $T_{\rm eff} = 4622_{-60}^{+61}$ K). From the $TESS$ light curve, we measure a rotation period of 12.48 days, and derive an age between about $500$ Myr and 1.1 Gyr. The photometric observations reveal three transiting mini-Neptunes ($R_b = 2.049^{+0.12}_{-0.080} R_\oplus$, $R_c = 2.701^{+0.092}_{-0.082} R_\oplus$, $R_d = 2.474^{+0.090}_{-0.082} R_\oplus $), with orbital periods of $P_b = 9.531$ days, $P_c = 51.699$ days, and $P_d = 84.839$ days. After modeling the three-planet system, an additional Earth-sized candidate is identified, TOI-712.05 ($P = 4.32$ days, $R_P = 0.81 \pm 0.11 R_\oplus$). We calculate that the habitable zone falls between 0.339 and 0.844 au (82.7 and 325.3 days), placing TOI-712 d near its inner edge. Among planetary systems harboring temperate planets, TOI-712 ($T = 9.9$) stands out as a relatively young star bright enough to motivate further characterization.


INTRODUCTION
Planetary systems undergo many evolutionary processes within the first billion years, but the mechanisms that shape their physical properties and orbital architectures are not yet fully understood. Early interactions with the protoplanetary disk, with other planets, and with the host star can all influence the properties of the mature systems we observe today, but it can be difficult to constrain early processes by observing mature systems. Observing young (< 1 Gyr) planetary systems gives us a chance to more directly observe evolutionary processes (e.g., tidal circularization, atmospheric escape, cooling) and constrain their timescales.
Most stars form in groups, the densest of which persist for hundreds of millions of years in loosely bound stellar clusters and associations (e.g., Briceño et al. 2007). Young stars rotate rapidly, which enhances magnetic field strength and leads to active surface regions (Hartmann & Noyes 1987). These features produce large stellar activity signals and make it difficult to distinguish planetary signals from those of stellar origin (e.g., Saar & Donahue 1997). Early investigations of young planetary systems were largely composed of radial velocity (RV) surveys targeting known moving groups, stellar clusters, and associations (Paulson et al. 2004;Sato et al. 2007;Lovis & Mayor 2007;Quinn et al. 2012Quinn et al. , 2014 but were limited to detecting giant planets due to the difficulty in detecting low-amplitude signals in the presence of increased stellar activity. While these RV surveys have contributed to our understanding of planetary formation and evolution, transiting planets provide additional constraints on the physical evolution of planets by measuring planet radii and paving the way for subsequent atmospheric observations. Moreover, precise, space-based photometry is capable of identifying small transiting planets even in the presence of the large photometric rotational modulation that is characteristic of young stars. Kepler (Borucki et al. 2010), and especially the ecliptic survey of K2 (Howell et al. 2014) made possible the detection and characterization of small planets orbiting young stars. K2 discoveries include planets in open clusters like the Hyades (Mann et al. 2016a;Livingston et al. 2018;Mann et al. 2018;Vanderburg et al. 2018) and Praesepe (Mann et al. 2016a;Obermeier et al. 2016;Pepper et al. 2017;Rizzuto et al. 2018;Livingston et al. 2019), as well as associations like Upper Scorpius (David et al. 2016;Mann et al. 2016b) and the Taurus-Auriga star forming region (David et al. 2019). These discoveries led to the observation that young planets tend to have larger radii than their mature counterparts. While there may be a bias against detecting small planets orbiting * ESA Research Fellow young stars, the large radii of young planets are broadly consistent with predictions for planets undergoing atmospheric mass loss, though the observed sizes imply more mass loss than the models predict. Further study of the evolution of the planet size distribution can constrain the timescales of the process and therefore the mechanisms driving it. This result demonstrates the promise of young systems, but the K2 sample is relatively small and most of the host stars are too faint or active for mass measurements with radial velocities. To discover new benchmark young systems amenable to detailed characterization, and to build a large enough sample of young planets to derive strong constraints from their ensemble properties, many more young stars must be surveyed.
Fortunately, ESA's Gaia mission (Gaia Collaboration et al. 2016, which has already released extremely precise astrometric and photometric measurements for more than a billion stars, has unveiled previously unknown young stellar groups and identified new members of existing groups (e.g., Kounkel & Covey 2019;Gagné et al. 2021). As an all-sky survey, NASA's TESS mission (Ricker et al. 2016) is searching the brightest of these stars for transiting planets. This has led to the discoveries of many more planets in known young groups (e.g., Newton et al. 2019;Plavchan et al. 2020;Rizzuto et al. 2020;Mann et al. 2020;Tofflemire et al. 2021). Moreover, because Sun-like stars spin down over time (e.g., Barnes 2007;Mamajek & Hillenbrand 2008), the TESS photometry used for transit searches can also identify young stars via their rapid rotation. In some cases, this leads to a revised young age for entire groups of stars (the Pisces-Eridanus stream; Meingast et al. 2019;Curtis et al. 2019b). In other cases, TESS rotation periods reveal a young age for planet-hosting field stars (Carleo et al. 2021;Zhou et al. 2021). Because TESS observes nearly the whole sky, it has begun to find the brightest, nearest examples of young planetary systems, which are amenable to detailed characterization through radial velocity masses and atmospheric measurement.
Though the size distribution of young planets has provided evidence for the physical evolution of small planets, many of them reside in systems with only one transiting planet, so their orbital architectures are not well constrained. In the same way that comparing the radii of young planets to their mature counterparts can illuminate the physical evolution of small planets, comparing the orbital architectures of young planetary systems to their mature counterparts can illuminate their orbital evolution. Kepler revealed that most short-period multi-planet systems are flat, with well-aligned, relatively circular orbits, and a slight preference to reside just outside of low-order mean motion resonances (Lissauer et al. 2011;Fabrycky et al. 2014;Hadden & Lithwick 2014). However, to explain the observed population of systems with only one transiting planet, it has been suggested that there exists a sample of planets that is either characterized by high mutual misalignment or single planet systems (e.g., Lissauer et al. 2011;Fang & Margot 2012;Ballard & Johnson 2016). On the other hand, recent analysis of the mutual inclinations of Kepler systems shows that they can be modeled by a single continuous distribution (Millholland et al. 2021). Planetary system architectures at young ages can help us understand the processes that have shaped the observed characteristics of the more mature Kepler sample. Multi-planet systems also offer an opportunity to constrain the system properties (e.g., masses, eccentricities, inclinations) via dynamical simulations, and to perform comparative studies with planets that have experienced the same stellar environment. For young planets, this may be particularly important in the context of constraining the timescale of thermal mass loss.
In this paper, we report the discovery of three mini-Neptunes transiting the adolescent, K4.5V, field star TOI-712, with periods of 9.53 days (TOI-712 b), 51.69 days (TOI-712 c), and 84.39 days (TOI-712 d). Additionally, we find a small candidate planet with a period of ∼4.3 days, but additional data is needed to confirm its planetary nature. In Section 2, we present our observational data obtained from TESS and ground-based facilities. In Section 3, we describe the stellar characterization, including constraints on the age of TOI-712. In Section 4, we provide our global analysis of the TOI-712 system, and report the derived properties of the planets. In Section 6, we present a dynamical stability analysis and explore the evolution of the TOI-712 system, including the potential survivability of a habitable zone (HZ) planet, given the location of TOI-712 d in the Venus zone at the inner edge of the HZ. We discuss our results in Section 7.
2. OBSERVATIONS 2.1. TESS Photometry NASA's Transiting Exoplanet Survey Satellite (TESS) is an all-sky photometric survey, launched 18 April 2018. Its primary mission is to detect small planets around nearby, bright stars, ideal for follow-up observations for further characterization. Using 4 cameras, each 24 • × 24 • , TESS observes the sky in 24 • × 96 • strips. In the two-year prime mission, these strips generally spanned from one of the ecliptic poles to near the ecliptic plane. Each observing sector lasts 27.4 days before stepping ∼27 • to the next anti-solar sector. After a year (13 sectors) observing the southern ecliptic hemisphere, the spacecraft oriented to observe the north for a year 1 . Overlap of the fields of view of successive sectors at high ecliptic lat-itudes results in longer observation time spans, extending to nearly continuous year-long coverage at the poles.
TIC 150151262 (TOI-712) was pre-selected for 2-minute cadence observations by the TESS mission on the basis of the star's size, brightness, and position on the sky near the southern ecliptic pole, the combination of which was expected to enable detection of small transiting planets. TOI-712 was observed on Camera 4 during the prime and extended missions between 2018 September 20 UTC and 2021 June 24 (sectors 3, 4, 6, 7, 9, 10, 11, 13, 27, 29, 30, 31, 32, 33, 34, 36, 37, and 39). During the extended mission, TOI-712 was selected for 20-second cadence observations in sector 27 as part of two Guest Investigator programs: G03068 (PI D. Kipping) and G03278 (PI A. Mayo).
The raw photometric data from TESS was processed by the TESS Science Processing Operations Center (SPOC) located at NASA Ames Research Center. The SPOC pipeline (Jenkins et al. 2010a(Jenkins et al. , 2016 extracts the light curve using Simple Aperture Photometry (SAP) and also corrects the light curve for instrument systematics using the Presearch Data Conditioning (PDC) algorithm (see Stumpe et al. 2012). The resulting light curves are available through the Mikulski Archive for Space Telescopes (MAST) archive. After filtering for stellar variability and residual noise, the light curve was searched for transit signals using the SPOC Transiting Planet Search (TPS, Jenkins 2002;Jenkins et al. 2010bJenkins et al. , 2020, which revealed two periodic transit events with periods of 9.53 days and 51.69 days. An initial limb-darkened transit model was fitted to the data (Li et al. 2019) and a suite of diagnostic tests were conducted to make or break confidence in the planetary interpretation of the signatures (Twicken et al. 2018). The two candidates passed all the diagnostic tests and the sources of the transits were localized to within 1.9500 ± 4.1489 arc sec and 2.0800 ± 4.5217 arc sec, respectively. The TESS Science Office reviewed the DV Data Validation reports and issued an alert in May 2019 (Guerrero et al. 2021). These were released as TESS Objects of Interest (TOIs), TOI-712.01 and TOI-712.02. An additional candidate with a period of 53.7 days was released as TOI-712.03, but only one transit event looked believable, and a system with 51-and 53-day planets would not be dynamically stable.
We note that the algorithms applied by SPOC are successful for the vast majority of stars, but the stellar variability of young stars often requires a boutique analysis to properly remove stellar signal and identify real transits (see, e.g., Rizzuto et al. 2017;Hedges et al. 2021). Given the apparent presence of this additional transit-like event, we chose to re-extract the light curve while correcting for spacecraft systematics, following Vanderburg et al. (2019). We started with the SPOC simple aperture photometry (SAP) light curve and fit a model to the time series consisting of a basis spline (with breakpoints spaced every 1.0 days to adequately model the stellar variability) and systematics parameters relating to the means and standard deviations of the spacecraft quaterion time series within each exposure (see Vanderburg et al. 2019), the fast timescale (band 3) cotrending basis vectors from the SPOC PDC analysis, and a high-pass-filtered time series of the flux in the SPOC background aperture. We fit for up to quadratic relationships in each systematics parameter. We solved for the best-fit coefficients of our free parameters using a matrix inversion technique, iteratively excluding 3σ outliers from the fit until it reached convergence. We then subtracted the best-fit systematics vectors from the SAP light curve, and corrected for diluting flux in the optimal aperture using the SPOC crowdsap parameter, to yield a final, corrected light curve. While there are 24 stars within two TESS pixels, none of them are brighter than T = 17, and as a result the aperture typically includes only about 2% dilution. The uncertainty in this correction, propagated to the derived parameters, is negligible compared to other sources of uncertainty. We performed the same extraction for the sector 27 twenty-second data, and then binned to 2 minutes. While the overall scatter in our extracted light curve is similar to the PDCSAP scatter, the systematics during periods of increased pointing jitter are reduced. This helped us identify additional transits of the third candidate, for which we determined a period of 84.839 days. We note that SPOC later identified another transit of this candidate in the extended mission. It was assigned a period of 678.7 days (eight times too long) and released as TOI-712.04, though we now know that both TOI-712.03 (53 days) and TOI-712.04 (679 days) correspond to the 84.8-day candidate. The three candidates have transit depths close to 1 mmag, implying sizes between 2 and 3 R ⊕ , and durations ranging from 1.7 to 5.7 hours, as expected for the range of orbital periods. Details of the transit fitting are discussed in Section 4. Our extracted light curve is shown in Figure 1.

Ground-based Time-series Photometry
We conducted ground-based photometric follow-up observations for TOI-712 as part of the TESS Follow-up Observing Program 2 (TFOP; Collins 2019). We used the TESS Transit Finder, which is a customized version of the Tapir software package (Jensen 2013), to schedule our transit observations. The observations are described below and summarized in Table 1.

LCOGT
We observed seven TOI-712 b and three TOI-712 c full or partial transits from the Las Cumbres Observatory Global Telescope (LCOGT; Brown et al. 2013) 1.0 m network. The We include all observations in the global fit (see Section 4) with the exception of the egress observation from 2019-12-25, which suffered from poor precision. It was nevertheless an important observation that confirmed that the transit occurrs at the location of the target star and is not a nearby eclipsing binary unresolved by TESS.

ASTEP
Four transits of TOI-712 c were observed in 2020 and 2021 with the Antarctica Search for Transiting ExoPlanets (ASTEP) program on the East Antarctic plateau (Guillot et al. 2015;Mékarnia et al. 2016 −5 0 5 Figure 1. The upper two panels show the raw TESS light curves from year 1 and year 3. TESS two-minute cadences are plotted in gray and the in-transit points for each planet are colored according to the legend. The spline model is shown in black. The middle two panels show the flattened light curves after removal of the spline, and the best-fit EXOFASTv2 model is plotted for each planet. The lower panels show the phase folded transits and residuals for each planet. Colored circles indicate binned data, and the colored lines are the best-fit EXOFASTv2 models.

Figure 2.
Twelve ground-based follow-up transits with EXOFASTv2 best-fit models of TOI-712 b (purple lines) and TOI-712 c (pink lines). Label indicate the instrument, the filter, and the date of observation (see Table 1). plate splits the beam into a blue wavelength channel for guiding, and a non-filtered red science channel roughly matching an Rc transmission curve. The data are processed on-site using an automated IDL-based pipeline described in Abe et al. (2013).

Hazelwood
We observed an ingress of TOI-712 c from Hazelwood Observatory near Churchill, Victoria, Australia. The 0.32 m telescope is equipped with a 2184 × 1472 SBIG STT3200 camera. The image scale is 0. 55 pixel −1 , resulting in a 20 × 14 field of view. The images were calibrated and the photometric data were extracted using AstroImageJ.
This event was the same as the one observed by LCOGT on 2019-12-25. It shows evidence for systematics in transit, and we therefore do not include it in our fit, but as the first observation of this candidate, we note that it represented a valuable confirmation that the transit occurs on the target star and is not a nearby eclipsing binary.

CHIRON Spectroscopy
We obtained three spectra of TOI-712 with the CHIRON high-resolution echelle spectrograph on the 1.5 m SMARTS telescope at the Cerro Tololo Inter-American Observatory (CTIO), Chile (Tokovinin et al. 2013). Observations were taken on UT 2019-Sep-20, 2019-Oct-31, and 2021-Apr-10 in slicer mode, which employs a fiber-fed image slicer to deliver light to the spectrograph, achieving a resolving power of R∼80, 000 across the wavelength range 4100-8700 Å. The CHIRON spectra can be used to characterize the stellar properties, search for massive outer companions, and rule out false positive scenarios that would induce large radialvelocity (RV) variations or exhibit line profile variations.
We extracted the RVs of TOI-712 by fitting the line profiles of each spectrum, which were measured via least-squares deconvolution (LSD) of the observed spectra against synthetic templates (Donati et al. 1997). The RVs, which reveal no significant velocity variation beyond the uncertainties (∼50 m s −1 ) are listed in Table 2. We also find no evidence for a composite spectrum or line profile variations.

High-resolution Imaging
We searched for close visual companions to TOI-712 in I band using the HRCam speckle imager on the 4.1 m Southern Astrophysical Research (SOAR) telescope (Tokovinin 2018;Ziegler et al. 2018) on UT 2019-Oct-16, and simultaneously at 562 nm and 832 nm using the Zorro 3 speckle imager on the 8 m Gemini South Telescope on UT 2020-Jan-10. We achieved contrast ratios better than 5 magnitudes at a separation of 0. 1 and > 7 magnitudes beyond ∼0. 5 in the Zorro 832 nm data, which have a field of view of 2. 5 × 2. 5 and are therefore sensitive to companions within about 1. 2 of TOI-712. Outside of 1. 5, Gaia can exclude the presence of stellar sources bright enough to produce the observed transit signals (Ziegler et al. 2018). Closer stars unresolved by Gaia will often manifest as elevated astrometric noise (see, e.g., Evans 2018), but TOI-712 has RUWE = 0.994, providing no evidence for excess astrometric motion 4 . The 5σ contrast curves from our speckle imaging are shown in Figure 3.

Spectroscopic Classification
We extract the effective temperature (T eff ), surface gravity (log g), and metallicity ( Fe/H ) of TOI-712 from the CH-IRON spectrum by matching against a library of observed spectra that have been previously classified by SPC (Buchhave et al. 2012), with interpolation performed via a gradientboosting regressor implemented in the scikit-learn python module. We find T eff = 4767 ± 50 K, log g = 4.63 ± 0.10, and Fe/H = −0.15 ± 0.08. The projected rotational velocity, v sin i , is derived following Gray (2005) and Zhou et al. (2018), fitting broadening kernels to the instrumental, macroturbulent, and rotational line profiles. We estimate v sin i = 2.0 ± 1.0 km s −1 and v mac = 1.4 ± 1.0 km s −1 , though the instrumental resolving power is only 3.7 km s −1 ; for such slow rotation, these estimates may suffer systematic offsets. The projected rotational velocity is consistent with the equatorial rotation velocity calculated from the stellar radius and rotation period, but it is not precise enough to confidently state that the stellar spin and orbital axes are well aligned.

Gyrochronology
The TESS light curve of TOI-712 shows rotational modulation with an amplitude of ∼2%, indicative of active regions on the stellar surface like those expected for a young star. We estimate the rotation period with multiple techniques, 3 https://www.gemini.edu/sciops/instruments/alopeke-zorro/ 4 See Lindegren et al. (2018) Figure 3. The 5σ contrast curves of TOI-712, observed by the SOAR speckle imager in the I-band (green), and the Zorro speckle imager at 562nm (blue) and 832nm (pink). Zorro at 832nm achieves contrast ratios better than 5 magnitudes at 0. 1, 7 magnitudes at 0. 5, and nearly 8 magnitudes at 1. 2. An eclipsing binary would have to be within 7 magnitudes to produce an event as deep as TOI-712 b.
first using a Lomb-Scargle periodogram (Lomb 1976;Scargle 1982). We assign the rotation period to the period of maximum power, with the uncertainty taken to be the half width at half maximum. However, the analysis is slightly complicated by the long time span of observations. Visual inspection of the light curve ( Figure 1) reveals changes in the morphology of the rotational phase curve on a relatively short timescale. The first year of TESS observations includes two epochs of a more complex, double-peaked phase curve and two epochs of a single-peaked phase curve. The third year is primarily single-peaked. If the dominant spot latitudes also change between the prime and extended missions, we might expect a slight change in the measured rotation period due to differential rotation, which is both predicted (e.g., Küker & Rüdiger 2011) and measured for some rapidly rotating stars (e.g., Collier Cameron 2007). A periodogram of the full time series would then be sensitive to the phasing of the signals even though we should not expect or require them to be in phase. Indeed, this is what we see in Figure 4, which shows that the periodogram of the combined TESS light curve only allows a few discrete periods, while individual years of data each exhibit broader peaks that are consistent with one another. We therefore measure the rotation period in each year separately. From the TESS year 1 light curves, we measure a stellar rotation period of 12.505 ± 0.441 days, and the TESS year 3 light curves yields a rotation period of 12.452 ± 0.629 days. Additionally, WASP photometry spanning five observing seasons from 2008 to 2012 measures a rotation period of 12.503 The TESS light curve shows a strong rotational modulation. The Lomb-Scargle periodogram of the combined TESS light curves is shown in pink, the periodogram for the extended mission is shown in teal, and the periodogram for the combined TESS light curves is shown in navy.. Each season of the WASP observations is shown in grey. The rotation period of ∼12.5 days is consistently measured in the seventeen sectors observed by TESS and the five seasons observed by WASP. Right: Normalized TESS light curve detrended using the methods described in Vanderburg et al. (2019), phase folded over the measured stellar rotation period extracted from the combined primary and extended TESS observations. The phased points are colored based on the time of observation, showing that while the phase curve morphology changes slightly over time, the measured period is consistent with data from both Year 1 (blue tint) and Year 3 (pink tint).
days. The periodograms from the TESS and WASP observations are shown in Figure 4, and the binned TESS light curve, phase-folded to the measured rotation period, is shown in Figure 4. As a check on the Lomb-Scargle periods, we also fit the full TESS light curve using a Gaussian process with a quasi-periodic kernel. This yields a period consistent with the Lomb-Scargle results, P rot = 12.538 days. We adopt a weighted average of rotation periods from individual seasons, P rot = 12.454 days.
Finally, we caution that we are observing the rotation period at the latitude of the active regions. The unknown strength of differential rotation and current spot latitudes contribute additional (systematic) uncertainties if interpreting the observed rotation period as the equatorial rotation period. If differential rotation in TOI-712 is similar to that of the Sun, then the true equatorial rotation period could be different by 10% or more.
We adopt the activity-rotation-age relation as described in Mamajek & Hillenbrand (2008), to estimate the age of TOI-712. In particular, we use the equation: where P rot is the rotational period of the star, in days, t is the stellar age in Myr, B − V is the color index, a = 0.407 ± 0.021 and b = 0.325 ± 0.024 are gyrochronolgy constants, c = 0.495 ± 0.010 is the 'color singularity', and n = 0.566 ± 0.008 is the time-dependent power law. Using the stellar rotation period above, and B − V = 1.265 ± 0.152, we estimate a gyrochronological age of 492 ± 103 Myr. How-ever, we caution that the formal uncertainty is quite likely underestimated. First, systematic errors may affect the determination of P rot , as described above. More importantly, TOI-712 is outside of the color range (0.5 < B −V < 0.9) in which the Mamajek & Hillenbrand (2008) relation is well calibrated, as their sample consisted of primarily Sun-like stars and TOI-712 is relatively red. More recent investigations of the gyrochronology of K dwarfs, enabled by precise photometry of adolescent stars from K2 and TESS, suggest that in contrast to the G dwarfs, K dwarfs experience an epoch of stalled spin-down. While it is true that K dwarfs are expected to have rotation periods similar to TOI-712 at 500 Myr, they may also have the same rotation period at twice that age (Curtis et al. 2019a). As a result, from its rotation alone, we can only estimate the age of TOI-712 to be between about 0.5 and 1.1 Gyr.

Comparison with stellar populations
Since recent work has identified many new members of known stellar groups, we explore whether TOI-712 might belong to a known association. We use the online tool BANYAN Σ (Gagné et al. 2018), to assess possible membership in known groups, but the result suggests that it is not part of any previously known groups.
Similarly, there are likely many sparse groups that are currently unrecognized, but the large volume of public data (from Gaia, TESS, broadband photometric surveys, and spectroscopic surveys) makes possible a targeted search for comoving, coeval companions to a star of interest. We searched for evidence of a previously unknown comoving population of stars via Friend Finder 5 . Within 25 pc of TOI-712, we find only two stars with similar kinematics, but this is likely no more than would be expected by chance. TOI-712 has either escaped from its birth environment or formed in such a small environment that almost none of its stellar siblings remain within < 25 parsecs. While no convincing kinematic matches were found, the age can also be constrained by the precise UV photometry provided by GALEX (e.g., Shkolnik et al. 2011;Findeisen & Hillenbrand 2010). Friend Finder illustrates that the near-UV (NUV) flux of TOI-712 is consistent with the NUV sequence of the stars populating the Hyades cluster (see Figure 5), which has a well constrained age of ∼625 Myr (see, e.g., Perryman et al. 1998;Martín et al. 2018) . Older stars should exhibit lower levels of NUV flux. This provides supporting evidence that the age of TOI-712 is consistent with the age of the Hyades, though if NUV flux and rotation evolve together, it is possible that NUV flux would plateau while spin-down stalls. If that were the case, then the NUV flux of Hyades K dwarfs could be similar to that of somewhat older stars. Richey-Yowell et al. (2019) do show that the NUV flux is inversely related to stellar age in K dwarfs, but that study included no stars between Hyades age and field stars. NUV sequences in slightly older clusters like NGC 752 (1.3 Gyr; Agüeros et al. 2018) and Ruprecht 147 (2.7 Gyr; e.g., Torres et al. 2020;Curtis et al. 2020) might help resolve the age of TOI-712, but because of their dis-5 https://github.com/adamkraus/Comove tance, GALEX provides only non-informative upper limits for most of the K dwarf members.
Given the age ambiguity associated with the stalled spindown of K dwarfs, we ultimately adopt an age of 830 +250 −280 Myr and bounded between 0.3 and 1.3 Gyr, as determined by the global fit described in Section 4.
3.2.3. The stellar flare rate of TOI-712 TOI-712 was observed with 20-second cadence in sector 27, and during that time a flare was observed at time BTJD 2054.672 with a maximum flux excursion of ∼10%. Though it lasted only three 20-second cadences, two independent flare-finding algorithms (Medina et al. 2020;Feinstein et al. 2020a,b) identify it as a high-probability flare. We calculate an energy in the TESS bandpass of 8.75 × 10 32 ergs, and we derive a flare rate of approximately 10 −3 flares per day above an energy of 2.32 × 10 32 ergs. While there are no well calibrated flare rates for adolescent K dwarfs, this flare rate falls slightly below the observed rates for younger K dwarfs (Feinstein et al. 2020b). This does not resolve the gyrochronology age degeneracy, but the presence of any significant flares further supports that this is an adolescent star.
Finally, we note that this flare lasted only 1 minute, which implies that there may be flares in the two-minute data that last only one cadence and appear as outliers of a few percent or less. This should not be surprising, as many previous studies with TESS and Kepler data have focused on active stars with flares lasting many minutes (Günther et al. 2020;Medina et al. 2020;Feinstein et al. 2020b) or longer (e.g., Davenport 2016), but there is an indication that flare energies in the optical correlate with flare duration (Maehara et al. 2015). While it is beyond the scope of this work, TESS data could offer a means to calibrate flare rates of adolescent stars, which may be prone to short duration and low energy flares that appear as a single modest outlier in two-minute cadence. Distinguishing single-cadence flares from instrumental effects is difficult, but a statistical approach may be feasible, particularly given the large number of TESS light curves of both young and old stars.

EXOFASTV2 GLOBAL FIT
To derive the properties of TOI-712 and its planets, we fit the system with the exoplanet modeling software EXO-FASTv2 (Eastman et al. 2019). We used all available 2minute photometric data from TESS and the ground-based light curves described in Section 2.2 in conjunction with broadband flux measurements (Table 3), the Gaia DR2 parallax, spectroscopic stellar parameters, and the MESA Isochrones and Stellar Tracks (MIST) stellar isochrones (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015Choi et al. 2016;Dotter 2016a). In our Markov Chain Monte Carlo (MCMC) fit,  we set priors on the stellar metallicity derived from the CH-IRON spectra, the parallax measured by Gaia, and the upper limit on extinction from the dust maps of Schlegel et al. (1998) and Schlafly & Finkbeiner (2011). At each step in the MCMC, the quadratic limb-darkening coefficients are interpolated from the values tabulated by Claret (2017). Following Ford (2006), we require the Gelman-Rubin statistic of all parameters to be < 1.01 and the number of independent draws to exceed 1000 to indicate convergence. The results of our global fit are reported in Table 4, and the best-fit models are plotted with the phase-folded light curves in Figure 1. Notably, the radii (2.049 R ⊕ , 2.701 R ⊕ , 2.474 R ⊕ ) suggest that all three planets are mini-Neptunes, likely to possess significant volatile envelopes. The eccentricities of the outer planets are relatively low, and while the eccentricity of TOI-712 b is not inconsistent with zero, it is poorly constrained. The medians of the marginalized posteriors are reported in Table 4, but the eccentricity posterior for planet b is broad and non-Gaussian. This is important to note, because while we present further analysis based on the broad EXOFASTv2 posteriors, our expectation based on Kepler multi-planet systems is that all three planets will have eccentricities much closer to zero than the medians of the posteriors reported in Table 4 (see, e.g., Moorhead et al. 2011;Van Eylen & Albrecht 2015;He et al. 2020). While multiplanet systems can exhibit transit time variations due to mutual dynamical interactions, the EXOFASTv2 transit times are consistent with a linear ephemeris.

Validation of the TOI-712 planets
In the previous section, we presented a model of the planetary system orbiting TOI-712, but have not formally validated the signals as bona fide planets. In particular, we note that while we detect the transits of TOI-712 b and c from the ground in small apertures that localize the transit signal to within a few arcseconds of TOI-712, there is no groundbased transit detection of TOI-712 d. Given the large TESS pixels, the signal could potentially be produced by a faint, nearby eclipsing binary (NEB). Gaia detects 25 stars within 2 pixels (42. 0), though the centroid analysis made available in the SPOC Data Validation report rules out NEBs more distant than 18. 0. The brightest of the remaining stars is 8.1 magnitudes fainter in the TESS bandpass, and therefore contributes a flux fraction of only 0.0006 (0.6 mmag). The transit of TOI-712 d is deeper than 1 mmag, so none of the nearby stars can be responsible for the transit signal. The high resolution imaging and Gaia astrometry discussed in Section 2.4 rule out nearly all stellar companions within 1 arcsecond, and the CHIRON RVs exclude many close, bound eclipsing binaries that would have accelerated TOI-712 during the 600-day observing timespan.
Previous analysis has shown that the overwhelming majority of Kepler candidates in multi-planet systems are real planets, and moreover that systems with more than 2 candidates are even more likely to be real (e.g., Latham et al. 2011;Lissauer et al. 2012). The same cannot be assumed for TESS, mainly because of the increased potential for NEBs due to its large pixels and, in some cases, crowded fields. However, our follow-up observations have excluded the presence of contaminating eclipsing binaries even more confidently than is the case for Kepler candidates. As a result, we consider TOI-712 b, c, and d to be bona fide planets.

TOI-712.05
Upon subtracting the best-fit EXOFASTv2 three planet model from the TESS light curve, we perform a Transit Least Squares (TLS, Hippke & Heller 2019) analysis on the residuals to search for the signature of any additional transiting planets. TLS reveals a ∼4.32-day candidate signal. We refer to this signal as TOI-712.05. A new EXOFASTv2 fit including all four candidates converged, with no change to the stellar parameters or properties of the other three candidates.   The derived properties of this candidate are shown in Table  Table 5, and the phase folded TESS light curve is shown in Figure 6. If real, TOI-712.05 would be slightly smaller than the Earth, R P = 0.81 R ⊕ . We present the results for TOI-712.05 separately because we have less confidence that it is real. Notably, the signal is not clearly present in the full frame images (FFIs). While this could be due to smearing of the short-duration transit in the 30-and 10-minute cadence data, it may also be due to a difference in the apertures used to produce the FFI and 2minute light curves. If this were the case, it could suggest that the signal originates on a nearby star at the edge of the 2-minute aperture. To test this possibility, we ran two experiments. First, we extracted FFI light curves for each of the stars within 2.5 and folded them at the ephemeris of TOI-712.05. None showed signs of an eclipse that would cause the observed signal. Because the sensitivity of the FFIs to a shallow, short duration event may be compromised we also extracted 2-minute light curves using different apertures.
Planetary candidates that are released as TOIs also undergo ephemeris matching against all other signals detected by the pipeline, in order to test for a false positive introduced by a nearby transit-like signal. TOI-712.05 was not released as a TOI, so we performed this test using all signals detected by MIT's Quick Look Pipeline (QLP; Huang et al. 2020a,b), and we did not find any ephemeris matches.
It is also possible that SPOC introduced the transit signal while processing the light curve since TOI-712 is relatively active. As a check, we ran TLS on the unprocessed SAP light curve to see if the transit is still present. After flattening, and removing the transits from the three planets, the phase folded SAP light curve does show the transit signal for TOI-712.05.
A final concern is that the signal could be introduced by stellar activity. If that were the case, we might expect that the signal would not appear consistent in time due to evolution of the stellar signal, which is apparent by eye in the unflattened TESS light curve. We tested this by comparing the signal in Year 1 to the signal in Year 3. To within the uncertainties, the two are consistent, as would be expected for a transiting object.
While TOI-712.05 passed the above tests, it remains a low signal-to-noise event, so we present it here as a candidate planet. Its planetary nature may be solidified as the signal strength increases in future TESS extended missions, or perhaps by additional analysis of existing data.

HABITABLE ZONE AND VENUS ZONE
As described above, the three detected planets in the system are all mini-Neptune in terms of size. While they are too large to expect Earth-like compositions, their locations with respect to the Habitable Zone (HZ) and the Venus Zone (VZ) can motivate further analysis and future observations. The HZ boundaries have been previously calculated from Earthbased climate models that determine the radiative balance conditions for retaining surface liquid water Kasting et al. (1993); Kopparapu et al. (2013Kopparapu et al. ( , 2014. The conservative HZ (CHZ) is defined by the runaway greenhouse limit at the inner boundary, and the maximum CO 2 greenhouse limit at the outer boundary (Kane et al. 2016). The optimistic HZ (OHZ) is an extension to the conservative HZ based on empirical evidence regarding retention of surface liquid water on Venus (inner boundary) and Mars (outer boundary) during their evolutionary histories (Kane et al. 2016). The location of the HZ boundaries are sensitively related to the stellar parameters (Kane 2014), provided in Table 4. Adopting the methodology of Kopparapu et al. (2014), we calculate the CHZ boundaries for the TOI-712 system as 0.436-0.802 au, and the OHZ boundaries as 0.344-0.846 au. The HZ regions are represented in Figure 7, which shows a top-down view of the TOI-712 system. Because the eccentricities are poorly constrained, we plot circular orbits, though note that modest Though none of the TOI-712 planets lie unambiguously within the HZ, they do lie within the VZ, defined as the region interior to the runaway greenhouse limit and exterior to the atmospheric mass loss limit . Such planets provide exceptionally useful targets for atmospheric characterization, both individually and in the context of describing the processes that occur within the VZ (Kempton et al. 2018;Ostberg & Kane 2019). The TOI-712 planets are unlikely to have rocky surfaces, but their atmospheric response to the local radiation environment can provide useful information regarding the processes driving atmospheric erosion. Moreover, because the planetary system extends (at least) to the inner edge of the HZ, their presence can place dynamical constraints on the existence of exterior planets that may lie within the HZ.

DYNAMICAL ANALYSIS
The orbital parameters of the TOI-712 planets (see Table 4) provide an opportunity to examine the orbital dynamics that exists within the system, and the dependence on their proximity to each other. Such dynamical interactions are particularly important for assessing eccentricity constraints for the known planets, and for placing dynamical limits on the pres-ence of other potential planets in the system (Gladman 1993;Chambers et al. 1996;Hadden & Lithwick 2018). Sampling from the posteriors of the EXOFASTv2 fit, we performed N-body integrations using the Mercury Integrator Package (Chambers 1999) and adopted a methodology similar to that described by Kane & Raymond (2014); Kane (2016Kane ( , 2019. We used a hybrid symplectic/Bulirsch-Stoer integrator with a Jacobi coordinate system (Wisdom & Holman 1991;Wisdom 2006), and adopted a time resolution of 0.2 days to adequately sample the orbit of the inner planet. Each of the simulations were executed for a duration of 10 7 simulation years.
Shown in the top panel of Figure 8 is a 300,000 year segment of the results from a simulation initialized with eccentricities near the median of the posteriors. While the true eccentricities may be closer to zero, the system remains stable even for relatively high eccentricities, and this simulation illustrates some of the mutual interactions experienced by the planets. The outer planets (c and d), which lie near, but not within, the 8:5 mean motion resonance (MMR), exhibit high frequency angular momentum exchanges that cause frequent fluctuations in their eccentricities. In turn, planets c and d interact with the innermost known planet (b), resulting in high amplitude oscillations in their respective eccentricities with a period of ∼50,000 years. These complex dynamical interactions raise the question of whether the known planets may exclude potential terrestrial planets from residing within the HZ. To test for this, we conducted an exhaustive suite of 10 7 year dynamical simulations using the same methodology as described above-with initial conditions drawn randomly from the EXOFASTv2 posteriors-but with the insertion of a hypothetical Earth-mass planet located in the semi-major axis range of 0.3-0.9 au, in steps of 0.05 au. The simulations were allowed to run until either the simulation completed, or planet-planet interactions caused the Earth-mass planet to be swallowed by the star or ejected from the system. The bottom panel of Figure 8 shows the results from this simulation in terms of the percentage of the simulations that survived the full 10 7 years as a function of the semi-major axis for the test locations of the Earth-mass planet, indicated by the solid line. Orbits beyond ∼0.7 au are robustly stable, and some planets closer than 0.7 au survive near MMRs. Since there are three interior planets, the locations of MMR are numerous and complicated, but we highlight one example of a resonance stable location corresponding to the 3:2 MMR of an HZ planet with planet d. Simulation Survival (%) Figure 8. Results of the dynamical simulations, described in Section 6. Top: The eccentricity variations of the known planets for a period of 300,000 years, showing the stability of the system and the high frequency eccentricity oscillations of planets c and d. Bottom: The dynamical simulation results for an Earth-mass planet located in the semi-major axis range 0.3-0.9 au, shown as percentage survival of the simulation as a function of semi-major axis (solid line). As for Figure 7, the CHZ and OHZ are shown in light green and dark green, respectively. The vertical dashed line indicates an example stable resonance location (3:2). teractions to constrain planetary masses and orbital architectures. This avenue for characterization is especially important when other methods like radial velocity monitoring are difficult-for example, due to stellar activity, faint host stars, or small or long-period planets that induce low-amplitude orbital motion. Because short-period planets are more likely to transit, relatively few systems include long-period planets, and most that do were discovered by Kepler orbiting faint stars. As a result, TOI-712 stands out as one of the brightest long-period planetary systems currently known (see Figure 9). While many mini-Neptunes transit bright stars, most reside close to their star with high equilibrium temperatures, so TOI-712 represents an opportunity to study temperate mini-Neptunes, which are likely to have experienced different formation environments and evolutionary histories than their short-period counterparts.
Previous investigations of young planets have suggested that they tend have larger radii than their older counterparts (e.g., Mann et al. 2017;David et al. 2021), which is qualitatively consistent with atmospheric mass loss, but the details of the process and the timescale over which it occurs are not precisely known. Moreover, most known young planets have relatively short orbital periods for which photoevaporation is likely an important process (e.g., Owen & Wu 2013;Jin & Mordasini 2018), so whether a mechanism like corepowered mass loss (Ginzburg et al. 2018) drives significant atmospheric evolution in temperate planets is an open question. We illustrate in Figure 10 that the radii of the TOI-712 planets are smaller than the youngest planets, but there are  Figure 9. Known validated small planets (RP ≤ 6.0R⊕) in multiplanet systems (as of 29 October 2021, via NASA's Exoplanet Archive). The host stars' K-magnitudes are plotted against the orbital periods. The points are colored based on the planet radius. TOI-712 d is one of the longest period, small planets orbiting a bright star. The brightness of the host allows for the planetary radii to be well constrained through photometric observations with relatively few observed transits, and also enables future characterization with other facilities. no young planets with insolation fluxes similar to the outer TOI-712 planets.
For TOI-712, it is unclear if the planets are small because they have experienced significant atmospheric loss, or because they simply formed smaller than the young planets discovered to date. The first possibility was discussed in Section 3.2; if they are older than 1 Gyr, their small sizes may not be a powerful constraint. On the other hand, the outer planets should not experience significant mass loss from photoevaporation, and David et al. (2021) show evidence for evolution in the radius gap as late as 2 to 3 Gyr, so even if the planets are ∼1 Gyr old, core-powered mass loss may operate on long enough timescales that they can provide constraints on the process, particularly at long periods where fewer planets are known. Even at short periods, where the radius valley is well characterized, TOI-712 b is of note. At R p = 2.049 R ⊕ , it sits near the edge of the radius valley for its size and incident flux. If the timescale for core-powered mass loss is as long as 2 Gyr or more, it could be losing atmosphere now and ultimately be destined to become a rocky world. Future advances in understanding K dwarf gyrochronology and activity may provide a more precise age and a more definitive conclusion about the evolution of these planets.
If TOI-712 is as young as gyrochronology suggests, then the small planet sizes are more likely to reflect the formation process rather than (or in addition to) the mass loss history. In particular, TOI-712 c and d orbit at greater distances from their host star than other known young planets, and may differ in composition. Under the assumption that mini Nep- tunes form in situ, local conditions in the protoplanetary disk during formation would differ for the long-period TOI-712 planets, the outer of which resides at the edge of the habitable zone. More data, such as transmission spectroscopy, could shed light on the compositional differences between hot and temperate planets. While the brightest examples of such planets are statistically likely to be old, adolescent examples like TOI-712 c and d, for which core-powered mass loss may not have had time to strip significant atmosphere, could provide complementary information about formation by probing primordial compositions.

Future characterization of TOI-712
As discussed above, young planets-including those orbiting TOI-712-open avenues for observational constraints on the dynamical and physical evolution of planetary systems. While young stars can be challenging to observe, future observations of the TOI-712 system remain promising.
Our modest precision RVs suggest some variation beyond the internal uncertainty estimates, which is consistent with expectations for an adolescent star with rotating surface features. Given the rotation (v sin i = 2 km s −1 ) and the amplitude of photometric modulation (∼1%), we estimate there should be tens of meters per second RV variation from stellar surface features (see, e.g., Saar & Donahue 1997). Since the predicted RV amplitudes induced by the planets are only a few meters per second (see Table 4), attempts to measure the planet masses with stabilized spectrographs would be difficult without significant advances in stellar activity modeling.
Dynamical interactions between the planets, if detected via transit timing variations or transit duration variations, could help determine their masses and eccentricities. There is currently no evidence for transit timing variations, but TOI-712 . The log of the TSM is plotted as a function of equilibrium temperature, Teq. Each planet is colored based on the effective temperature of the host star, Teff. The dashed lines denote the suggested JWST TSM cuttoff for mini-Neptunes as described in Kempton et al. (2018). We show the calculated TSM for planets ranging between 2.0 ≤ RP ≤ 4.0 R⊕. TOI-712 b, TOI-712 c, and TOI-712 d are shown by the starred points, with a black border. There are very few temperate planets with a high TSM, so while TOI-712 c and TOI-712 d fall below the suggested threshold, they are promising targets for studying atmospheric evolution in the absence of strong irradiation.
will be observed by TESS in future extended missions, providing longer baselines over which to measure TTVs. Moreover, the signal-to-noise ratio (SNR) of the TOI-712.05 transit should increase if it is a real companion, and we can therefore more confidently confirm or refute its existence. Similarly, additional transiting planets in the system may become apparent with the inclusion of more data. In particular, we showed in Section 5 that there are dynamically stable orbits for small planets in the habitable zone of TOI-712. We also know that systems of small planets tend to be packed, so it would not be surprising to find an additional planet between TOI-712 b and c (10 P 51 days). While planets larger than TOI-712 b (2 R ⊕ ) would have likely been detected already, smaller planets may be revealed in future TESS data. Better knowledge of the system architecture and improved measurements of orbital properties will make dynamical modeling more powerful, which also feeds back into improved planetary properties. As discussed in Section 7.1, the atmospheric composition of the outer planets could inform our understanding of the formation of temperate mini-Neptunes. We note that TOI-712 is located 1.4 • from the Southern Ecliptic Pole, in the JWST continuous viewing zone, and it stands out as one of the brightest stars hosting long-period transiting planets (Figure 9). The star is relatively small, and as a result, transmission spectroscopy of the planets may be possible. Following Kempton et al. (2018), we calculate the Transmission Spectroscopy Metric (TSM) of the planets to assess the prospect for transmission spectroscopy with JWST. We find TSM b = 53.22, TSM c = 41.35, and TSM d = 31.88. Based on expected yields from TESS, Kempton et al. (2018) suggested that a survey with JWST should adopt a TSM threshold of 92 for planets of this size. However, this threshold is purely based on the SNR and ignores the outsized scientific value of some individual planets that happen to have lower expected SNRs. TOI-712 c and d are examples of planets that would not meet the cutoff for such a survey, but may still be worth observing. Figure 11 shows the TSM values of transiting mini-Neptunes (2.0 ≤ R P ≤ 4.0R ⊕ ) as a function of equilibrium temperature and stellar effective temperature. It is clear that although the planets lie below the TSM threshold for the predicted TESS yield, they stand out among the real sample as two of the best temperate mini-Neptunes. When also considering the age of the system, these planets become more scientifically valuable and may form the basis for a comparative study of temperate planets as a function of age.

CONCLUSION
Young planetary systems offer a glimpse into planetary evolution, most of which occurs within the first billion years after formation. In this paper, we presented the discovery of three planets transiting the adolescent K dwarf TOI-712. We determine the stellar age from the TESS rotation period, but because adolescent K dwarfs experience stalled spin-down, can only determine its age to be in the range of 500 Myr to 1.1 Gyr. There is no evidence that it is a part of a moving group, so we conclude that TOI-712 is an adolescent field star.
Our global analysis determined precise planetary radii (R b = 2.049 +0.12 −0.080 R ⊕ , R c = 2.701 +0.092 −0.082 R ⊕ , R d = 2.474 +0.090 −0.082 R ⊕ ), indicating that all three planets likely possess significant volatile envelopes. We identify an additional 4.32-day rocky planet candidate (TOI-712.05; R P = 0.81R ⊕ ), but do not claim it as a validated planet due to its low SNR. TOI-712 d orbits at the inner edge of the optimistic habitable zone, making it the first such planet orbiting a young star. As a result of their cool temperatures and relatively bright host star, TOI-712 c and d represent promising targets to study the evolution of temperate small planets, for which thermal mass loss cannot shape planetary atmospheres.
Dynamical simulations of the system reveal it to be stable and find that small planets in the conservative habitable zone could also survive. Future TESS observations may be sensitive to such objects, if they exist, and can confirm or refute the small inner rocky planet candidate.