TESS Hunt for Young and Maturing Exoplanets (THYME) VI: an 11 Myr giant planet transiting a very low-mass star in Lower Centaurus Crux

Mature super-Earths and sub-Neptunes are predicted to be $\simeq$Jovian radius when younger than 10 Myr. Thus, we expect to find 5-15$R_\oplus$ planets around young stars even if their older counterparts harbor none. We report the discovery and validation of TOI 1227 b, a $0.85\pm0.05R_J$ (9.5$R_\oplus$) planet transiting a very low-mass star ($0.170\pm0.015M_\odot$) every 27.4 days. TOI~1227's kinematics and strong lithium absorption confirm it is a member of a previously discovered sub-group in the Lower Centaurus Crux OB association, which we designate the Musca group. We derive an age of 11$\pm$2 Myr for Musca, based on lithium, rotation, and the color-magnitude diagram of Musca members. The TESS data and ground-based follow-up show a deep (2.5\%) transit. We use multiwavelength transit observations and radial velocities from the IGRINS spectrograph to validate the signal as planetary in nature, and we obtain an upper limit on the planet mass of $\simeq0.5 M_J$. Because such large planets are exceptionally rare around mature low-mass stars, we suggest that TOI 1227 b is still contracting and will eventually turn into one of the more common $<5R_\oplus$ planets.


INTRODUCTION
Young planets offer a window into the early stages of planet formation and evolution. Planets younger than 100 Myr are particularly useful for this work, as planetary systems likely evolve most rapidly in the first few hundred million years after formation (Lopez & Fortney 2013;Owen & Wu 2013). Populations of such planets are critical to understanding planetary migration (Nelson et al. 2017), photoevaporation (Raymond et al. 2008;Owen & Lai 2018), and atmospheric chemistry (e.g., Segura et al. 2005;Gao & Zhang 2020). Given the timescale for planet formation (1-10 Myr; Yin et al. 2002;Alibert et al. 2005), planets aged 30 Myr can even tell us about the conditions of planets right after their formation.
The number of known young transiting planets has grown significantly in recent years (e.g., Obermeier et al. 2016;David et al. 2018;Benatti et al. 2019;Newton et al. 2021). This growth was primarily driven by a combination of the K2 and TESS missions surveying nearby young clusters and star-forming regions (Ricker et al. 2014;Van Cleve et al. 2016), improvements in filtering variability in young stars (e.g., Aigrain et al. 2016;Rizzuto et al. 2017), and more complete identification of young stellar associations from Gaia kinematics (e.g., Cantat-Gaudin et al. 2018;Kerr et al. 2021). Despite this progress, there are still only a handful of transiting planets at the youngest ages of 30 Myr (David et al. 2016;Mann et al. 2016a;David et al. 2019a;Plavchan et al. 2020;Rizzuto et al. 2020) and a few candidate * UNC Chancellor's Science Scholar † 51 Pegasi b Fellow ‡ NSF GRFP Fellow nontransiting planets from radial velocity surveys (e.g., Johns-Krull et al. 2016;Donati et al. 2017). Models predict that gas giant planets younger than < 50 Myr will be larger and brighter than their older counterparts (Linder et al. 2019). At 10-20 Myr, progenitors of mature Jovian-mass planets are expected to be 1.2-1.6 R J and sub-Neptunes ∼ 1R J (Linder et al. 2019;Owen 2020). Completeness curves from the best search pipelines (Rizzuto et al. 2017), and the discovery of much smaller planets (e.g., Mann et al. 2018;Zhou et al. 2020) demonstrate that giant planets are readily detectable even in the presence of complex stellar variability common to young stars. Thus far, transit surveys of young stars identified only a few Jovian-radius planets in the youngest associations (David et al. 2019b;Bouma et al. 2020;Rizzuto et al. 2020).
As part of the TESS Hunt for Young and Maturing Exoplanets survey (THYME; Newton et al. 2019), our team searches TESS data using a specialized pipeline to identify young planets missed by standard searches (e.g., Rizzuto et al. 2020), and checks previously-identified planet candidates for signs of membership in a young association (e.g., Mann et al. 2020;Newton et al. 2021). We identified TOI 1227 (2MASS J12270432-7227064, TIC 360156606) as a member of a young association, with a planet candidate (TOI1227.01) identified by the TESS mission, and astrometry consistent with membership in the Lower Centaurus Crux (LCC) region of the Scorpius-Centaurus (Sco-Cen) OB association (Goldman et al. 2018;Damiani et al. 2019). The TESS -identified transit signal is 2% deep, suggesting a Jovian-sized planet orbiting a pre-main-sequence low-mass star. As giant planets with periods <100 days are rare around older stars of similar mass (< 1% Bonfils et al. 2013), validation of this planet would provide strong evidence of radius evolution.
In this work, we present validation, characterization, and age estimates for TOI 1227 b, a 0.85R Jup planet orbiting a 11 Myr pre-main-sequence M5V star (0.17M ) every 27.26 days. In Section 2, we detail the photometric and spectroscopic follow-up of the planet and host star. We demonstrate that the star is a member of a recently identified substructure of LCC and derive an updated age of 11 ± 2 Myr for the parent population in Section 3. We estimate parameters for the host star in Section 4, parameters of the planet in Section 5, and combine these with our ground-based follow-up to statistically validate the planet in Section 6. We place the large size of TOI 1227 b in context with its age and host star mass and explore its likely evolution in Section 7. We conclude in Section 8 with a summary and discussion of follow-up and implications for future searches for young exoplanets. In all three sectors, the target fell on Camera 3. The first two sectors had 2 min cadence data as part of a search for planets around M dwarfs (G011180; PI Dressing). Sector 38 had 20 s cadence data, as it was known to be a young TOI at that phase (G03141; PI Newton). We initially used the 2 min cadence for all analysis for computational efficiency and TESS cosmic-ray mitigation 1 (not included in 20 s data), but also ran a separate analysis using the Sector 38 20 s data (see Section 5. In both cases, we used the Pre-Search Data Conditioning Simple Aperture Photometry (PDCSAP; Stumpe et al. 2012;Smith et al. 2012;Stumpe et al. 2014) TESS light curve produced by the Science Process Operations center (SPOC; Jenkins et al. 2016) and available through the Mikulski Archive for Space Telescopes (MAST) 2 .

Identification of the transit signal
The planet signal, TOI1227.01, was first detected in a joint transit search of sectors 11 and 12 (one transit in each) with an adaptive wavelet-based detector (Jenkins 2002;Jenkins et al. 2010). The candidate was fitted with a limb-darkened light curve (Li et al. 2019) and passed all performed diagnostic tests (Twicken et al. 2018). Although the difference image centroiding failed to converge, the difference images indicate that the transit source location was consistent with the location of the host star, TOI 1227. A search of the residual light curve failed to identify additional transiting planet signatures. The TESS Science Office reviewed the diagnostic test results and issued an alert for this planet candidate as a TESS object of interest (TOI) on 2019 August 26 (Guerrero et al. 2021). A third transit was detected in the Sector 38 TESS data, consistent with the expected period and depth.
We searched for additional planets using the Notch and LoCoR pipelines, as described in Rizzuto et al. (2017) 3 . This included using the significance of adding a trapezoidal Notch to the light curve detrending, as characterized by the Bayesian information criterion (BIC). The method was more effective than periodic methods for finding planets with 3 transits, as was the case for HIP 67522 b (Rizzuto et al. 2020). The transits were quite clear from the BIC test. However, no additional significant signals were detected.
TOI 1227 is a relatively faint star (T = 13.8) in a crowded region (see Figure 1). Within 2 ( 6 TESS pixels), there are four sources brighter than TOI 1227, one of which is T = 8.4 (TIC 360156594). Two sources brighter than TOI 1227 are within 1 . While the TESS aperture was small (4-6 pixels over the 3 sectors), background contamination was likely to be a problem for the TESS data. Correcting for contamination and confirming the planet were the major motivations for our ground-based follow-up.

Ground-Based Photometry
The photometry is summarized in Table 1. We provide details by instrument below.

Goodman/SOAR
We observed two transits of TOI 1227 b using the SOAR 4.1 m telescope at Cerro Pachón, Chile, with the Goodman High Throughput Spectrograph in imaging mode (Clemens et al. 2004). For both observations, we used the red camera. In this mode, Goodman has a default 7.2 diameter circular field of view with a pixel scale of 0.15 pixel −1 . The first transit was observed on 2020 Jan 15 with the SDSS i' filter, and the second transit was observed on 2021 Mar 28 with the SDSS g' filter. Both nights had photometric conditions for the full duration of the transit observations. Although TOI 1227 is much fainter in g' than in i', we used a 120 s exposure  The TESS PSF in this region has a full-width half-max of approximately 1.9 TESS pixels. The large TESS PSF combined with a faint source meant that there were high levels of contamination around TOI 1227.
time for both filters to resolve any flares and sample the ingress. We selected a smaller region-of-interest in the readout direction (1598 pixels) to decrease readout time to 1.7 s. Due to the small pixel scale, no defocusing was required and counts were well below half-well depth for both the target and all but one target in the field of view (HD 108342). About 1 hour into the 2021 Mar 28 transit observations, the SOAR guider began to fail, causing the target to shift 5 pixels per exposure and forcing us to shift it back manually every ∼10 exposures. During egress, the guider completely failed and had to be reset before re-acquiring the field. We positioned the target back to its starting location and continued the observations as normal. This resulted in poorer photometric precision than normally achievable with Goodman/SOAR, and a 15 min gap in the data near the end of the egress. We applied bias and flat-field corrections before extracting photometry for TOI 1227 using a 10-pixel radius aperture and used an annulus of 30-60 pixels to determine the local sky background. The aperture center for each exposure was the stellar centroid, calculated within a 10-pixel radius of the nominal location. We repeated this on eight nearby stars that were close in brightness to the target and showed little or no photometric variation when compared to other stars in the field. We corrected the target light curve using the weighted mean of the comparison star curves.

LCO Photometry
We observed a total of seven transits with 1 m telescopes in the Las Cumbres Observatory network (LCO ; Brown et al. 2013). These were all observed with Sinistro cameras, with a pixel scale of 0.389 pixel −1 . Two transits were observed using the SDSS g filter, one with SDSS r , two with i , and two with z s . Exposure times are given in Table 1. Because of the difficulty of scheduling observations of a long-period planet (∼27 d) and weather interruptions, all LCO observations covered only part of the transit.
The images were initially calibrated by the standard LCOGT BANZAI pipeline (McCully et al. 2018). We then performed aperture photometry on all datasets using the AstroImageJ package (AIJ; Collins et al. 2017). The aperture varied based on the seeing conditions at the observatory, but we generally used a 6-10 pixel radius circular aperture for the source and an annulus with a 15-20 pixel inner radius and a 25-30 pixel outer radius for the sky background. For all observations, we centered the apertures on the source and weighted pixels within the aperture equally. Because the event was detected on the source, the usual check of nearby sources for evidence of an eclipsing binary was not necessary. Light curves of nearby sources are available with the extracted light curves and further details on the follow-up at ExoFOP-TESS 4 (ExoFOP 2019).
Except for a 2021 Apr 24 g' transit, all photometry showed the expected transit behavior. Two transits had poor coverage (no out-of-transit baseline), but still showed a shape consistent with ingress. We discuss the 2021 Apr 24 transit in more detail in Section 6.4.

ASTEP Photometry
On UTC 2021 Apr 24, one additional transit was observed with the Antarctica Search for Transiting ExoPlanets (ASTEP) program on the East Antarctic plateau (Guillot et al. 2015;Mékarnia et al. 2016). The 0.4 m telescope is equipped with an FLI Proline science camera with a KAF-16801E, 4096 × 4096 frontilluminated CCD. The camera has an image scale of 0.93 pixel −1 resulting in a 1 o ×1 o corrected field of view. The focal instrument dichroic plate splits the beam into a blue wavelength channel for guiding, and a non-filtered red science channel roughly matches an R c transmission curve.
Due to the extremely low data transmission rate at the Concordia Station, the data are processed on-site using an automated IDL-based pipeline described in Abe et al. (2013). The calibrated light curve is reported via email and the raw light curves of about 1,000 stars are transferred to Europe on a server in Rome, Italy, and are then available for deeper analysis. These data files contain each star's flux computed through 10 fixed circular 4 https://exofop.ipac.caltech.edu/tess/target.php?id=360156606 aperture radii so that optimal light curves can be extracted. For TOI 1227, a 9.3-pixel radius aperture gave the best result.
TOI 1227 was observed under mixed conditions, with a windy clear sky, and air temperatures between −62 o C and −66 o C. The Moon was ∼ 90% full and present during the observation. A strong wind (∼ 6 ms −1 ) led to telescope guiding issues at the beginning of the observation and prevented us from observing the ingress. The data points corresponding to these issues were removed from the analysis. However, the resulting light curve showed the expected egress.

Corrections for second-order extinction
Since atmospheric extinction is strongly colordependent, changes in airmass produce color-dependent flux losses that depend on the spectral energy distribution (SED) of the target. These color terms are weaker in redder and narrower passbands (Young et al. 1991). The effect is often small, as comparison stars typically span a range of colors, and observations are timed for favorable airmass. However, TOI 1227 was unusual both because of its long orbital period (27.3 d) and red color; of stars <5 from the target and differing by <1 magnitude in G, the reddest has a Gaia color of BP −RP = 2.1 while TOI 1227 has BP − RP = 3.3.
Because we had few options for mitigating secondorder extinction (color terms), we corrected for the effect following Mann et al. (2011). To summarize, we estimated T eff for all comparison stars based on their Gaia colors and the tables from Pecaut & Mamajek (2013) 5 . We then combined the relevant BT-SETTL models (assuming Solar metallicity) with an airmassdependent model of the atmosphere above the observing site (Rothman et al. 2009) and the appropriate filter profile. The output was a predicted change in flux from second-order extinction alone. The trend was negligible (< 1 mmag) for i' and z s observations, so we did not apply a correction. The effect was small but nonnegligible for r (>1 mmag), and significant for g'-band observations (as large as 5 mmag), so we applied it to those datasets.

Goodman/SOAR
We observed TOI 1227 with the Goodman spectrograph (Clemens et al. 2004) Kesseli et al. (2017). Because of extinction around TOI 1227 and imperfect flux calibration from Goodman, we fit out a 2nd-order polynomial on the Goodman spectrum to better match the template. The Goodman spectrum is also higher resolution than the template (R of 5880 vs 2000). The inset highlights the Li 6708Å line, which is absent in the spectrum of older stars like the template, but strong in TOI 1227's spectrum.
the 1200 l/mm grating in the M5 setup, and the 0.46 slit rotated to the parallactic angle. This setup gave a resolution of R 5880 spanning 6150-7500Å. We obtained the first spectrum on 2019 Dec 12 (UT) shortly after the target was alerted as a TOI, to check for signs of youth. We obtained three back-to-back exposures of TOI 1227, each with an exposure time of 800 s. The resulting spectrum showed the Li and Hα expected for a young star and TiO features consistent with an M5 dwarf, as we show in Figure 2. Once we confirmed the planetary signal (based on ground-based transits), we obtained an additional spectrum on 2021 Feb 23 (UT) under clear conditions. Our goal was a spectrum with a signal-to-noise ratio (SNR) 100 per pixel across the full wavelength range, to use for stellar characterization (Section 4). To this end, we used an exposure time of 1800 s for five back-to-back exposures.
To better determine the age of TOI 1227, we also observed 13 nearby stars that are likely part of the same grouping in LCC as TOI 1227 (see Section 3). The 13 targets observed with Goodman were selected from the parent sample to map out lithium levels from M0-M5, which can constrain the age of the population. These targets were observed between 23 Feb 2021 and 24 Apr 2021 using an identical setup to the observations for TOI 1227. Exposure times were set to ensure a SNR greater than 50 (per pixel) around the Li line and varied from 90 s to 420 s per exposure with at least 5 exposures per star for outlier (cosmic ray) removal.
Using custom scripts, we performed bias subtraction, flat fielding, optimal extraction of the spectra, and mapping pixels to wavelengths using a 5th-order polynomial derived from the Ne lamp spectra obtained right before or after each spectrum. Where possible, we applied a small linear correction to the wavelength solution based on the sky emission or absorption lines. We stacked the individual extracted spectra using the robust weighted mean. For flux calibration, we used an archival correction based on spectrophotometric standards taken over a year. Due to variations on nightly or hourly scales, this correction is only good to 10%.

HRS/SALT
We observed TOI 1227 with the Southern African Large Telescope (Buckley et al. 2006) High Resolution Spectrograph (Crause et al. 2014) on six nights (2020 May 08 to 2020 June 17). On each of the six visits, we obtained three back-to-back exposures. We used the high-resolution mode with an integration time of 800s per exposure. The resulting spectral resolution was R ∼ 46,000. The HRS data were reduced using the MIDAS pipeline (Kniazev et al. 2016), which performs flat fielding, bias subtraction, extraction of each spectrum, and wavelength calibration with arclamp exposures.
We measured radial velocities from the SALT spectra by computing spectral line broadening functions (BFs; Rucinski 1992) with the saphires python package 6 (Tofflemire et al. 2019). The BF results were computed from the linear inversion of an observed spectrum with a narrow-lined template. We computed the BF for 16 individual spectral orders in the SALT red arm that range from 6400 to 8900Å; the remaining orders had strong telluric contamination. We used a 3100 K, log g= 4.5 PHOENIX model as our narrow-lined template (Husser et al. 2013). The BFs only contained one profile; we did not detect a secondary star. The BFs from each order were then combined, weighted by the SNR. We fit the combined BF with a rotationallybroadened profile to determine the stellar radial velocity and v sin i * . Radial velocities from each epoch, corrected for barycentric motion, are presented in Table  2. The mean and standard deviation of the v sin i * measurements were 17.5±0.3 km s −1 , although this does not include corrections for activity/magnetic broadening or macroturbulence (∼1-3km s −1 , Sokal et al. 2020).
All observations were done following commonly used strategies for point-source observations with IGRINS 7 . Each target was placed at two positions along the slit (A and B), taking an exposure at each position in an ABBA pattern. Individual exposure times were between 120s and 425s, and the (total) times per epoch were between 720s and 2280s to achieve peak SNR 100 per resolution element in the K-band. To help remove telluric lines, we observed A0V standards within 1 hour and 0.1 airmasses of the observation of TOI 1227.
We reduced IGRINS spectra using version 2.2 of the publicly available IGRINS pipeline package 8 (Lee et al. 2017), performing flat fielding, background removal, order extraction, distortion correction, wavelength calibration, and telluric correction using the A0V standards and an A-star atmospheric model. We used the spectrum right before telluric correction to improve the wavelength solution and provide a zero point for the RVs.
We extracted the radial velocities of TOI 1227 using the IGRINS RV code 9 ) with a 3000 K PHOENIX model (Husser et al. 2013) and the TelFit code to create a synthetic telluric spectrum (Gullikson et al. 2014). Barycentric-corrected RVs from each epoch are listed in Table 2. IGRINS RV provided an estimate of the rotational broadening (v sin i * = 16.65±0.24 km s −1 ) and the star's systemic velocity (13.3 ± 0.3 km s −1 ), the calibration of which is detailed in Stahl et al. (2021).

High-Contrast Imaging
We observed TOI 1227 with the Gemini South speckle imager, Zorro (Scott et al. 2018), on UT 2020 March 13 (program ID GS-2020A-Q-125). Zorro provided simultaneous two-color, diffraction-limited imaging, reaching angular resolutions of ∼0.02 in ideal conditions with a field of view of about 1.2 . TOI 1227 was observed in 17 sets of 1000x60 msec exposures by the 'Alopeke-Zorro visiting instrument team with the standard speckle imaging mode in the narrow-band 5620Å and 8320Å filters ([562] and [832]). All data were reduced with the pipeline described in Howell et al. (2011). No additional sources were detected within the sensitivity limits ( Figure 3). The limiting contrasts for binary companions are set by the redder [832] filter, ruling out equal-brightness companions at ρ = 0.03 and excluding companions fainter than the target star by ∆[832] < 1.5 mag at ρ = 0.05 , ∆[832] < 4.5 mag at ρ = 0.10 , and ∆[832] < 5.0 mag at ρ = 0.2 . Converting M [832] to M G using the color-magnitude relations of Kraus et al. (in prep), the corresponding mass and physical scale limits implied using the 10 Myr isochrones of Baraffe et al. (2015) can rule out equal-mass companions at ρ = 3 AU, M > 50M Jup at ρ > 5 AU, M > 15M Jup at ρ > 10 AU, and M > 14M Jup AU at ρ > 20 AU.

Limits on Wide Companions from Gaia EDR3 imaging
To identify potential wide binary companions that might exert secular dynamical influences on the TOI 1227 planetary system, we searched the Gaia EDR3 catalog (Gaia Collaboration et al. 2021) for comoving and codistant neighbors. While several candidate neighbors are found at projected separations of ρ > 500 (Section 3.1), none are found at ρ < 400 (ρ 40, 000 AU), the typical outer bound for observed binary companions  for Sun-like stars (Raghavan et al. 2010) and well beyond the typical separation of mid-M binaries (Winters et al. 2019).
The nominal brightness limit of Gaia is G lim < 21 mag, which at an age of 10 Myr approximately corresponds to a limit of M lim 14M Jup . Ziegler et al. (2018) and Brandeker & Cataldi (2019) have shown that Gaia is sensitive to companions with ∆G = 6 mag at ρ > 2 and ∆G = 4 mag at ρ > 1 . Thus, Gaia would have identified almost all companions down to the separation regime where it meets the contrast limit of the Zorro speckle imaging (Section 2.6).

Archival Spectra from ESO
With the goal of better characterizing TOI 1227's age, we retrieved the publicly available reduced spectra from the European Southern Observatory (ESO) archive for stars in the same population as TOI 1227 (Section 3). We downloaded any spectra of candidate members that included the Li line. In total, 11 objects had spectra from HARPS at the 3.6m telescope (La Silla Observatory), FEROS at the 2.2 m telescope (La Silla Observatory), and/or X-shooter at VLT (Cerro Paranal Observatory).

Archival Photometry, Astrometry, and Velocities
For TOI 1227 and all candidate members of the parent population (Section 3.1) we retrieved parallaxes, positions, proper motions, and photometry from Gaia Early Data Release 3 (EDR3; Gaia Collaboration et al. 2021; Lindegren et al. 2021;Riello et al. 2021). From EDR3 we also retrieved the Renormalised Unit Weight Error (RUWE) for all stars. RUWE value is effectively an astrometric reduced χ 2 value, normalized to correct for color and brightness dependent effects 10 . RUWE should be around 1 for well-behaved sources, and higher values (RUWE 1.3) suggests with the presence of a stellar companion (Ziegler et al. 2020;Wood et al. 2021).
We downloaded velocities for candidate population members from a general Vizier/SIMBAD search, taking the most precise value for stars with multiple measurements. This yielded velocities for 21 stars. The full list of sources for these velocities is included in Table 5.

THE AGE AND MEMBERSHIP OF TOI1227
Our Goodman spectrum of TOI 1227 showed strong Hα emission and lithium absorption. Combined with TOI 1227's high position on the color-magnitude diagram (CMD), this set an upper age limit of < 30 Myr. However, we were able to derive a more precise estimate of the age of TOI 1227 from its parent population. To this end, we first identified likely members of the same association, then estimated the group age using rotation, lithium abundance, and fitting isochrones to the color-magnitude diagram.

The parent association of TOI 1227
To identify any co-moving and coeval population to TOI 1227, we first used the BANYAN-Σ tool (Gagné et al. 2018) 11 , providing as input the IGRINS systemic velocity and the Gaia EDR3 position, parallax, and proper motion. This yielded a membership probability of 99.3% for Cha, with 0.4% for LCC and 0.3% for the field. However, TOI 1227 is 25 pc away from the core of Cha, while nearly all known members of Cha are packed in a sphere 5 pc across (Murphy et al. 2013, also see Figure 4). The BANYAN algorithm preferred placing Cha over LCC likely because of better agreement in U V W despite poor XY Z agreement. However, the current BANYAN model did not account for more recent findings of significant velocity substructure in LCC (Goldman et al. 2018;Kerr et al. 2021), which changes the U V W model for LCC and hence the membership probabilities.
TOI 1227 was more recently identified as a member of the "B" sub-population of LCC by Kerr et al. (2021) and the A0 sub-group of the Crux Moving Group (CMG) by Goldman et al. (2018). These two groups are effectively the same, but use different selection methods and treatment of sub-groups. Kerr et al. (2021) considered this group (LCC-B) part of the LCC substructure, while Goldman et al. (2018) considered the CMG a moving group associated with LCC with multiple subgroups (A0, A, B, and C). The XY Z positions, proper motions, and available radial velocities of LCC-B/CMG-A0 members are similar to TOI 1227 (better so than Cha; Figure 4), so we adopted the LCC-B/CMG-A0 as the initial parent population.
We looked for additional members using the FriendFinder algorithm 12 (Tofflemire et al. 2021). FriendFinder uses the Gaia EDR3 astrometry and radial velocity of a source (we used our IGRINS systemic velocity) to compute Galactic U V W and XY Z, then projects the U V W motion into the tangential velocities that would be expected for nearby stars if they were to share the same space motion. We selected stars with separations <10 pc from TOI 1227 and tangential velocities <1 km s −1 from the expected values calculated by FriendFinder. These cuts were based on the fact that 1 km s −1 1 pc Myr −1 and an estimated age for the population of 10 Myr. We estimated these cuts would yield 2 field interlopers based on the background population of Gaia stars.
The three candidate membership lists-those from Goldman et al. (2018), Kerr et al. (2021), and this workhave significant overlap in both Galactic position and proper motion (Figures 4 and 5). Of the 108 stars in any of the three lists, 40 were in all three, and 27 were in two of the three. Further, most of the stars missing from Goldman et al. (2018) were missing precise parallaxes in Gaia DR2 (the former used DR2 while the other two selections used EDR3). Many of the stars in the FriendFinder list but not in Kerr et al. (2021) were removed by one of the quality cuts imposed by the latter work (e.g., BP/RP flux excess). The FriendFinder list was also missing a small number of objects slightly further away from TOI 1227 because of our separation cut.
12 https://github.com/adamkraus/Comove To be inclusive, we adopted the combination of all three lists as the membership list for TOI 1227's association; these 108 stars are listed in Table 5.
The resulting population has two (contradictory) names in the literature, both of which are difficult to remember. Given the presence of a young transiting planet, it deserves a more memorable name than 'A0' or 'B'. Most of the members land within the Musca constellation, so we refer to the merged group of stars as the Musca group (or just Musca).
TOI 1227 resides at the center of Musca, both based on Galactic position ( Figure 4) and proper motion (Figure 5). The mean velocity of the members with literature RVs (12.8km s −1 ) is also in good agreement with the value for TOI 1227. A randomly selected star that matches Musca in XY Z has a 1% probability of matching Musca in proper motion and radial velocity by chance. The CMD also indicates a common age for stars in the population ( Figure 6), with TOI 1227 matching the sequence. Combined with the presence of lithium in its spectrum, TOI 1227's membership in Musca is unambiguous.

Rotation
To better measure the age of Musca (and hence TOI 1227), we measured rotation periods for all candidate members using TESS Full-Frame Images (FFIs) downloaded from the Mikulski Archive for Space Telescopes (MAST). First, we created initial light curves from the FFI cutouts centered on each candidate member. After background subtraction, we used the unpopular package (Hattori et al. 2021) to generate a Causal Pixel Model (CPM) of the telescope systematics for each star, which we subtracted from the initial light curve. We used the CPM curves because it does a better job preserving long-period signals than PDCSAP (Stello et al. 2016;Wang et al. 2017).
We searched the resulting single-sector light curves for rotation periods between 0.1-30 days using the Lomb-Scargle algorithm (Horne & Baliunas 1986), repeating for each available sector. We selected the rotation period from the sector that returned the largest Lomb-Scargle power. As a check, we phase-folded the resulting light curves along that rotation period to inspect each measurement by eye. A few stars had multiple peaks in their Lomb-Scargle periodograms, which might indicate an unresolved binary. In such cases, we took the stronger signal (which was always the shorter period). Each rotation measurement was assigned a quality score during the visual inspection following Rampalli et al. (2021), with Q0 indicating an obvious rotation signal, Q1 a questionable signal, Q2 a spurious detection, and Q3 The planet host, TOI 1227, is marked with a red X, and is in all three lists. The significant overlap between these three selections re-affirms that these are all the same population, just selected using slightly different methods. For reference, we included a green circle that includes most of the known Cha members, highlighting that Cha is a separate and more compact population. a light curve dominated by noise. In total, we measured usable rotation periods (Q0 or Q1) for 90 stars out of a sample of 108. Of the remaining 18, 11 were too faint to retrieve a reliable light curve, one showed a clear dipper pattern (and was identified as such by Tajiri et al. 2020), three had high flux contamination from nearby stars, and the remaining three had no significant period detection. Rotation period measurements are included in Table 5. The high detection rate is consistent with a young population and suggests a low rate of field-star interlopers in our selection. The rotation period distribution ( Figure 7) is consistent with the spread seen in the 10 Myr Upper Scorpius association (Rebull et al. 2018), and marginally consistent with the tighter rotation sequence at 40-60 Myr from Gagné et al. (2020). TOI 1227's rotation period (1.65 days) was consistent with the Musca sequence. Spectral Type Figure 6. Color-magnitude diagram of stars in Musca (stars associated with TOI 1227), with symbols by the selection method or reference. We removed stars with poor astrometric fits from Gaia (RUWE> 1.5; Gaia Collaboration et al. 2021). TOI 1227 is marked as an red "X", and lands within the tight sequence. The blue line indicates an empirical main-sequence for reference (Pecaut & Mamajek 2013).
This effectively sets a limit of 60 Myr for the group and further confirms TOI 1227 as a member.

Lithium
Lithium (Li) is a powerful indicator of age in young stars (Chabrier et al. 1996). M dwarfs deplete their lithium over 10-200 Myr at a rate that depends on their spectral type; this forms a region where stars have no surface lithium (i.e., the lithium depletion boundary or LDB). The location and size of the LDB are strongly sensitive to the association's age, and LDB ages are largely independent of the isochronal age (Soderblom et al. 2014). At the youngest ages (<20 Myr), lithium is only partially depleted, leading to a 'dip' in the Li levels short of a full boundary ( Figure 8; also see Rizzuto et al. 2015). The location and depth of the lithium dip also depend strongly on age. For our Li determinations, we measured the equivalent width of the Li 6708Å line for 22 stars using our highresolution ESO archival (11 targets) and Goodman (13 targets) spectra (two stars overlap). To account for the variations in resolution, v sin i * , and velocity between targets and the instrument used, we first fit nearby atomic lines with a Gaussian profile (e.g., iron lines for warmer stars and potassium lines for the M dwarfs). We used the width from these fits to define the bounds of the Li line. To estimate the pseudo-continuum, we iteratively fit the 6990-6720Å region excluding the Li line, each time removing regions >4σ below the fit (there were no emission lines in this region). We did not attempt to correct for contamination from the Fe line at 6707.44Å or broad molecular contamination in the cooler stars, which likely set a limit on the precision of our equivalent widths at the 10% level.
A single star, TIC 359357695, had two clear sets of nearly equal-depth lines (an SB2). Interestingly, this star is a known TOI (1880), indicating a roughly equalmass eclipsing binary. For the Li equivalent width, we measured each line individually with a manually-applied offset. We then combined the two equivalent widths.
In Figure 8b we compare the Li sequence for Musca to that from β Pic ( 24 Myr; Shkolnik et al. 2017) and Cha (3-5 Myr; Murphy et al. 2013). The Cha cluster has high Li levels over the full sequence, while β Pic showed a full depletion around M3-M5. Musca resides between these two, with a dip in Li levels around M3 but not full depletion; this effectively bounds the age of Musca between the two groups. Based on Li predictions from the Dartmouth Stellar Evolution Program (DSEP; Dotter et al. 2008) with magnetic enhancement (Feiden & Chaboyer 2012), we estimated the Li age to be 8-14 Myr (Figure 8).

Comparison to theoretical isochrones
We compared the candidate members of Musca to the PARSECv1.2S stellar isochrones (Bressan et al. 2012). To handle contamination from binaries and nonmember interlopers, we used a mixture model described in detail in Appendix A. We also removed stars with RUWE> 1.3 (likely binaries) and stars outside the colors covered by the isochrones (i.e., stars with M G > 11.8 and G − R P > 1.45). Many tight pairs were resolved in Gaia, but not in 2MASS or similar ground-based surveys. The mixture model is not set up to handle these cases; a data point has a single outlier probability independent of the band. Thus, we performed this comparison using Gaia magnitudes only. An inspection of the sequence using 2MASS magnitudes suggested that this would not change the derived age in any significant way.
The fit, shown in Figure 9, yielded an age of 11.6±0.5 Myr, consistent with 11.8 Myr from Goldman et al. (2018), 13.0±1.4 Myr from Kerr et al. (2021), and 12 Myr for the lower end of LCC from Pecaut & Mamajek (2016). The age errors were likely underestimated, as our fit did not fully account for model systematics (often 1-3 Myr, e.g., Bell et al. 2015). Indeed, the fit slightly underpredicted the luminosity of K2-K5 dwarfs and overpredicted the value for M4 and later. As an additional test, we repeated the fit using magnetic DSEP, which yielded a consistent age of 11.1±1.4 Myr, suggesting model differences are comparable to our measurement errors.

The age of Musca
Given the constraints from the isochrones, lithium, rotation, and the literature, we adopt an age of 11 Myr with a conservative uncertainty of ±2 Myr. We used this age for TOI 1227 b as well, as any age spreads are likely to be similar to, or smaller than, the adopted uncertainties.

HOST STAR ANALYSIS
We summarize properties of TOI 1227 in Table 3, the details of which we provide below. We first fit the spectral-energy-distribution (SED) following the methodology from Mann et al. (2016a) and highlighted in Figure 10. To briefly summarize, we com-pared the observed photometry and Goodman spectrum to a grid of template spectra drawn from nearby (likely reddening-free) young moving groups. The observed data was simultaneously fit with Solar-metallicity BT-SETTL CIFIST models (Baraffe et al. 2015). The atmospheric models were used both to estimate T eff and fill in gaps in the template spectra (e.g., beyond 2.3 µm). We combined this full SED with the Gaia EDR3 parallax to estimate both the bolometric flux (F bol ) and the total luminosity (L * ). With T eff and L * , we calculated the stellar radius (R * ) using the Stefan-Boltzmann relation. The fit included six total free parameters: the choice of template, A V , T eff , and three parameters that account for flux and wavelength calibration offsets between the Goodman spectra, stellar templates, and model spectra. To account for variability in the star, we added (in quadrature) 0.03 mags to the errors of all optical photometry. The resulting fit yielded T eff = 3072 ± 84 K, F bol = (7.87 ± 0.53) × 10 −11 erg cm −2 s −1 , −0.09 mag, and R * = 0.56±0.03 R . Swapping to the PHOENIX model grid from Husser et al. (2013) yielded nearly identical F bol and a higher (but consistent) T eff of 3145±67 K.
During the fit, the surface-brightness predictions of the atmospheric models were scaled to match the observed photometry. The scale factor is equivalent to R 2 * /D 2 , which combined with the Gaia EDR3 paral- Residual (σ) Figure 10. Example fit from comparing BT-SETTL models and young un-reddened templates to the observed photometry and Goodman spectrum. The red points are the observed photometry converted to fluxes using the appropriate zeropoint and uncertainty, with the vertical errors denoting the uncertainties and the horizontal errors the approximate filter width. Green points are synthetic photometry derived from combining the spectrum with the relevant filter profile. The bottom panel shows the residual in standard deviations. TOI 1227 was not detected in WISE W 4 and the W 3 point is not shown for clarity (but is used in the fit). The template/model shown here (black) is an M4.5-M5 star in the β Pic association and a 3075 K model (blue).
lax, provided an additional estimate of R * . The technique is similar to the infrared flux method (IRFM; Blackwell & Shallis 1977). IRFM-based radii of lowmass stars have not consistently agreed with empirical measurements (Casagrande et al. 2008;Boyajian et al. 2012), but more recent efforts have been more successful (Morrell & Naylor 2019), and our resulting R * estimate (0.57 ± 0.03 R ) was consistent with our estimate using the Stefan-Boltzmann relation.
We also derive T eff using our IGRINS spectra. While we consider the SED-based estimate to be reliable, we were concerned about significant spot coverage impacting our derived T eff (and hence R * ). Spot coverage would manifest as a cooler temperature at redder wavelengths, where the (cooler) spots have a larger impact on the total flux. Hot spots would have the opposite effect, which would also be visible as a wavelength-dependent T eff . The IGRINS data offers another advantage here; T eff determined from the low-resolution Goodman spectra was driven by the molecular bands, while the highresolution H-and K-band data from IGRINS are less sensitive to missing or erroneous molecular opacities. Thus, while the two fits use the same BT-SETTL models, the fits using IGRINS data were sensitive to different systematics in those models (Rajpurohit et al. 2018).
We fit the highest SNR IGRINS spectra obtained. We fit each order separately, but only included orders with low telluric contamination (determined by eye) and SNR >80. For each fit, we explored six free parameters: T eff , log g, radial velocity, a broadening parameter (including instrumental and rotational effects), the coefficients of the first two Chebyshev polynomials (used to correct for flux calibration offsets), and a factor that describes missing uncertainties in the data as a fraction of the observed spectrum. We fixed [Fe/H] to Solar, consistent with measurements of young regions around Sco-Cen (James et al. 2006). We fit all free parameters using the Monte Carlo Markov Chain code emcee (Foreman-Mackey et al. 2013). Each order was run with 30 walkers for 50,000 steps after an initial burn-in of 10,000 steps and with uniform priors bounded only by the model grid or physical limits.
Parameters besides T eff were effectively nuisance parameters in this analysis, as most other parameters (e.g., v sin i * ) were better determined with more empirical methods. We summarize the T eff results in Figure 11. T eff measurements from a given order were often precise (typical errors of 15-50 K), but varied between orders by 50-100 K; we considered the latter a more accurate reflection of the true errors, as systematics in the models change with wavelength. The T eff for each order was generally consistent with our fit to the SED and optical spectrum, indicating that our derived T eff was reliable and spots (while still present) did not significantly impact our derived T eff .
We repeated our fit to the IGRINS data after coadding all 11 spectra. Because the spectra were taken across the rotational phase, these likely sample different spot coverage patterns. However, the resulting spectrum was significantly higher SNR (> 200). The results were broadly consistent, with smaller errors in each order (10-30 K) but a similar variation between orders (50-100 K). The variation between orders was not consistent with the effect of spots, further suggesting that the spread in inferred T eff was driven by systematics in the models, rather than observational noise or surface inhomogeneities on the host.
To determine M * and as an additional check on our stellar parameters above, we compared the observed photometry to Solar-metallicity magnetic DSEP evolution models. We compared Gaia and 2MASS absolute magnitudes to the model predictions with an MCMC framework using emcee. We fit for age, A V , and M * . An additional parameter (f , in magnitudes) captured underestimated uncertainties in the data or models. We used a hybrid interpolation method that found the nearest age in the model grid and then performed linear one-  498 1.51 1.522 1.535 1.548 1.561 1.574 1.587 1.601 1.658 1.674 1.689 1. 2.097 2.121 2.146 2.172 2.198 2.225 2.253 2.281 2.311 2.34 ( m) 2900 3000 3100 3200 3300 T eff (K) Figure 11. Marginal probability distributions on T eff ('violin' plot; Waskom et al. 2020) for 24 IGRINS orders fit independently over the H-(left) and K-band (right). Each violin has an equal area distributed according to the MCMC posterior, so narrow/tall violins represent larger uncertainties than wider/shorter ones. The blue bar represents the 1σ distribution from our optical spectrum and SED fit.
dimensional interpolation in mass to obtain stellar parameters and synthetic model photometry. To account for errors from our nearest-age approach and ensure uniform sampling in age, we pre-interpolated bilinearly in age and mass using the isochrones package (Morton 2015). The interpolated grid was far denser (0.1 Myr and 0.01M ) than expected errors. To redden the model photometry, we used synphot (Lim 2020), following the extinction law from Cardelli et al. (1989). We placed Gaussian priors for age (11±2 Myr; derived from the population), A V (0.21±0.1 mag; from the SED fit), and T eff (3072±74 K; from the SED fit). All other parameters evolved under uniform priors. From the best-fit posteriors, we were able to interpolate additional posteriors on other stellar parameters from the evolutionary model, including R * and T eff .
The fit yielded age = 11.5 ± 1.0 Myr, M * = 0.165 ± 0.010 M , A V = 0.148 ± 0.074 mag, R * = 0.554 ± 0.007 R , and T eff = 3050±20 K. Changing the T eff prior to uniform yielded consistent, but less precise parameters (e.g., M * = 0.170 ± 0.015M ). As an additional test, we repeated the analysis using the PARSEC models, which yielded parameters somewhat discrepant with our SED and population analysis (age = 14.9 ± 2.1 Myr, T eff = 2900±30 K, and M * = 0.20±0.01 M ). PARSEC did not reproduce the observed colors of the coolest stars in Musca, which manifested as systematically high luminosities for a fixed color past M4 (see Figure 9). Thus, we prefer the magnetic DSEP fit. However, we adopted the more conservative mass (M * = 0.170 ± 0.015 M ), which was 2σ consistent with the PARSEC value.

Rotation period, rotational broadening, and stellar inclination
Canto Martins et al. (2020) reported a rotation period (P rot ) of 1.663±0.028 days for TOI 1227 using TESS data. Our analysis of rotation periods in Section 3 yielded a consistent 1.65±0.04 days, which we adopted for our analysis. We computed v sin i * as part of extracting radial velocities from IGRINS/Gemini and HRS/SALT spectra.
The mean value from the IGRINS data was 16.65±0.24 km s −1 , while the SALT data yielded a marginally inconsistent value of 17.8±0.3 km s −1 . We adopted the former value, as the IGRINS data are significantly higher SNR and less impacted by spots or molecular line contamination.
We used the combination of v sin i * , P rot , and R * to estimate the stellar inclination (i * ), and hence test whether the stellar spin and planetary orbit are consistent with alignment. A basic version of this calculation can be done by estimating the V term in v sin i * using V = 2πR * /P rot , although in practice it requires additional statistical corrections, including the fact that we can only measure alignment projected onto the sky. To this end, we followed the formalism from Masuda & Winn (2020). The resulting stellar inclination was consistent with alignment with the planet, yielding a limit on inclination of i * > 73 • at 95% confidence, and i * > 77 • at 68% confidence.

X-Ray Luminosity
TOI 1227 was detected in X-rays in a ROSAT pointing of the globular cluster NGC 4372 in 1993 (Johnston & Verbunt 1996). The X-ray source was listed as #6 in the NGC 4372 (IDed in SIMBAD as [JVH96] NGC 4372 6), however, the quoted values were not usable here, as they quote an X-ray flux assuming that the source is at the distance to NGC 4372 and only gave upper limits on the hardness ratios. The X-ray detection was reanalyzed by Voges et al. (2000) in the 2nd ROSAT PSPC Catalog, which IDs the X-ray source as 2RXP J122703.8-722702, situated 4.9 away from TOI 1227 with a positional error of ±9 (i.e. consistent with TOI 1227 being the source of the X-rays). They list a soft X-ray count rate of f X = (1.967 ± 0.626) × 10 −3 ct s −1 with hardness ratios HR1 = 0.06 ± 0.32 and HR2 = 0.58 ± 0.32 with an effective exposure time of 7399 s observed over UT 4-6 September 1993. Using the energy conversion factor equation of Fleming et al. (1995), this count rate and HR1 hardness ratio translates to a soft X-ray energy flux of f X = 1.70 ± 10 −14 erg s −1 cm −2 , which at the distance to TOI 1227 translates to a soft X-ray luminosity of L X = 10 28.32 erg s −1 .
Given our estimate of the star's bolometric luminosity (Table 3), this translates to a fractional X-ray luminosity of log(L X /L bol ) = −2.66. This is within the range of activity levels stars in the saturated regime display. However, X-ray levels tell us little about TOI 1227's age; M5V stars remain in the saturated well into field ages (West et al. 2015;Newton et al. 2017). We used MISTTBORN to fit for the five regular transit parameters: time of inferior conjunction (T 0 ), orbital period of the planet (P ), planet-to-star radius ratio (R p /R ), impact parameter (b), and stellar density (ρ ). For each wavelength observed, we included two linear and quadratic limb-darkening coefficients (q 1 , q 2 ) following the triangular sampling prescription of Kipping (2013). We included data from four unique bands: TESS, g', i', and z s , requiring eight limb-darkening parameters in total (r' was not used for reasons detailed below). Gas drag and gravitational interactions are expected to dampen out eccentricities and inclinations of extremely young planets like TOI 1227 b (Tanaka & Ward 2004), so we locked the eccentricity to zero (although we check this assumption later).
To model stellar variations, MISTTBORN included a Gaussian Process (GP) regression module, utilizing the celerite code (Foreman-Mackey et al. 2017). We initially followed the procedure in Foreman-Mackey et al. (2017), using a mixture of two stochastically driven damped simple harmonic oscillators (SHOs) at periods P GP (primary) and 0.5P GP (secondary). However, we found the second signal was poorly constrained, suggest-13 https://github.com/captain-exoplanet/misttborn ing a single SHO was sufficient. Instead, we adopted a fit that included three GP terms: the dominant period (ln P GP ), amplitude (ln Amp), and the decay timescale for the variability (quality factor, ln Q).
Stellar variation from spots is wavelength-dependent. This raised a problem when fitting multiple transits over such a wide wavelength range with a single GP amplitude. The most robust solution would be to fit the data using multiple GPs with a common period (e.g., GPFlow; Matthews et al. 2017). However, most of the ground-based transits were partials and many lacked significant out-of-transit baseline, so the GP is poorly constrained from an individual transit. Further, our GP kernel was able to fit nonsimultaneous data of multiple wavelengths, because the rotation signal can evolve during the 27 days between transits. In cases where we had simultaneous data from multiple filters (3 transits), we used only the more precise dataset. This cut excluded all r' data but still left us with nine observed transits (6 from the ground). We used the excluded observations for our false-positive analysis (specifically the test of chromaticity; Section 6).
When handling crowded regions and faint sources like TOI 1227 (see Figure 1), the SPOC reduction of TESS data has been known to over subtract the sky background (Burt et al. 2020), leading to a deeper transit depth. This was a potential issue for data collected before Sector 27, after which SPOC applied a correction to their reduction. 14 An initial fit of our data yielded TESS depths somewhat deeper than the ground-based transits, and the Sector 11-12 data deeper than the Sector 38 data. To correct for this, we re-analyzed the target pixel files and estimate that the PCDSAP light curves required an additive correction of 37e −1 pixel −1 for Sector 11 and 75e −1 pixel −1 for Section 12. No significant bias was seen for the Sector 38 data. We applied this correction to the light curve before running MISTTBORN, and to account for uncertainties we fit for two dilution terms (C S11 and C S12 ) applied to the TESS Sector 11 and 12 data following Newton et al. (2019): where LC undiluted was the model light curve generated by MISTTBORN. A negative dilution would correspond to ovsersubtraction of the background (i.e., if the corrections above were underestimated). b Derived parameters calculated using stellar parameters in Table 3.
We applied Gaussian priors on the limb-darkening coefficients based on the values from the LDTK toolkit (Parviainen & Aigrain 2015), with errors accounting for errors in stellar parameters and the difference between models used (which differ by 0.04-0.08). All other parameters were sampled uniformly with physically motivated boundaries (e.g., |b| < 1 + R P /R * , 0 < R P /R * < 1, and ρ * > 0).
We ran the MCMC using 50 walkers for 500,000 steps including a burn-in of 50,000 steps. This run was more than 50 times the autocorrelation time for all parameters, indicating it was more than sufficient for convergence.
All output parameters from the MISTTBORN analysis are listed in Table 4, with a subset of the parameter correlations in Figure 12. Combining the derived R P /R with our estimated stellar parameters from Section 4 yielded a planet radius of R P = 0.854 +0.067 −0.052 R J . The resulting fit from MISTTBORN reproduced the individual transits (Figure 13) as well as the overall light curve variations seen in TESS (Figure 14). The dilution parameters were both consistent with zero, suggesting our corrections to the PDCSAP curves were reasonable. The transit depth and stellar radius gave a planet radius of 0.859 +0.065 −0.052 R J (9.6 ± 0.7 R ⊕ ). The posterior showed a tail at larger radii, corresponding to a higher impact parameter (Figure 12). At 99.7% confidence TOI 1227 b is <1.2 R J (13.4 R ⊕ ).
Including the GP was computationally expensive, which led us to use the (binned) 120 s TESS data from Sector 38 (although 20 s data was available). As discussed above, a single GP fitting chromatic variability over multiple wavelengths and a long time period may be misleading. As a test, we refit the transit without the GP, removing the stellar variability before running the MCMC fit. For the TESS data, we fit the shape of the transits and the low-frequency variability simultaneously following Pepper et al. (2020), then divided out the variability model. For the ground-based transits, we fit each transit individually including the GP model as above (such that each transit can have a unique GP amplitude) and divided out the best-fit GP model from the light curve. Because most of the ground-based photometry had a limited out-of-transit baseline, only four of the transits were used.
We ran the GP-free fit through MISTTBORN for 200,000 steps after a burn-in of 20,000 steps. Priors were the same as from our earlier fit.
The GP-free fit was largely consistent with our GP fit, but with smaller uncertainties. The resulting fit parameters are listed in Table 4. We also show the phase-folded transit data in Figure 15.
We assumed zero eccentricity for both fits, but any nonzero eccentricity might (depending on the argument of periastron) manifest difference between our stellar density derived in Section 4 and the value from the fit to the transit. The transit-fit density agreed with the spectroscopic value. However, the spectroscopic and isochronal parameters are far less precise and eccentricities inferred this way are degenerate with the argument of periastron. Thus, we can only say that a low eccentricity (< 0.2) is preferred.
The transit-fit density is more than a factor of two more precise than the value estimated from our SED fit and stellar models (Section 4). If the e = 0 assumption is valid, the transit-fit density could be used to improve the estimated R * (e.g., Seager & Mallén-Ornelas 2003;Sozzetti et al. 2007). However, this provided only a marginal improvement on the R * uncertainty (0.0250R vs 0.030R ) due to a large error on M * . Since the M * estimate is model-dependent, we opted to keep the stellar parameters from Section 4.

FALSE-POSITIVE ANALYSIS
We considered the three most common false positive (non-planetary) scenarios to explain the transit signal: an eclipsing binary, a background eclipsing binary, and a hierarchical (bound) eclipsing binary. Other scenarios that may apply specifically to young stars (e.g., stellar variability) were quickly dismissed due to the depth, duration, and shape of the transit, as well its consistency over more than a year.

Eclipsing Binary
We compared the radial velocity (RV) data from IGRINS (Table 2) to the predicted velocity curve of a planet or binary at the orbital period of the outer planet (27.4 d). We did not use HRS data for this analysis. There may be zero-point offsets between velocities from IGRINS and HRS. HRS velocities also showed significantly higher velocity scatter (>500 m s −1 ), likely due to higher stellar variability in the optical.  A Jovian-mass planet would induce an RV signal of 220 m s −1 , while the scatter in the IGRINS velocities was 130 m s −1 (Figure 16), ruling out even a highly eccentric brown dwarf. There was a marginal detection of a ∼0.5 M J planet, particularly if we allow the planet's orbit to be slightly eccentric. However, a modest stellar jitter expected from rotation and long-term variation in the spot pattern could explain all of this variation. Denser sampling and a full radial velocity model including stellar noise would be more likely to yield a detection, but we can still place an upper limit of M P < 1.7 M J at >95% confidence. This rules out any possible brown dwarf companion.

Background Eclipsing Binary
Given a maximum eclipse depth of 50% and an observed transit depth of 2.5%, the star producing the observed signal must be within ∆m < 3.26 mags of TOI 1227. Since the depth is consistent across all filters, this constraint applies to all griz magnitudes. The multi-band speckle data (Section 2.6) ruled out such a companion or background star down to 90 mas. Vanderburg et al. (2019) detail a similar metric based on the ratio of the ingress time (T 12 ) to the time between first and third contact (T 13 ). However, because TOI 1227 b has a high impact parameter, this metric gives marginally worse constraints than the transit depth alone.
The proper motion of TOI 1227 was large enough to see 'behind' the source using digitized images (DSS) from the Palomar Observatory Sky Survey (POSS-I, POSS-II) and the Southern ESO Schmidt (SERC) Survey (patient imaging; Muirhead et al. 2012;Rodriguez et al. 2018 Figure 15. Phase-folded transit data from TESS and four ground-based transits after removing stellar variability. The TESS and SOAR i data have been binned for clarity. Dilution on TESS data has been removed. The transit shape varies slightly due to differences in limb darkening, but the transit parameters generally agree. transit imaging ( 0.8 ) and speckle imaging (0.09 ), and somewhat larger than the resolution of the POSS plates (1.7 ). No source was present in the DSS images at the modern location of TOI 1227. After subtracting the source using a nonparametric model of the PSF from StarFinder (Diolaiti et al. 2000), we were able to rule out any source <3 magnitudes fainter than TOI 1227. Additionally, both USNO-A2.0 and USNO-B reported no other nearby sources data (Monet et al. 2003), but were sensitive to stars >3 magnitudes fainter than TOI 1227 (in R).

Hierarchical Eclipsing Binary
To check for a bound companion eclipsing binary, we first ran the MOLUSC code (Wood et al. 2021). MOLUSC simulated binary companions following an empiricallymotivated random distribution in binary parameters. We limited the mass of the synthetic companions to between 10 M J to 0.17 M , as the goal was only to identify any unseen stellar or brown dwarf companion. For each synthetic binary, MOLUSC computed the corresponding velocity curve, brightness, and (sky-projected) separation, as well as enhancement to the noise that would be measured by Gaia. The results were compared directly to our observed high-resolution images, radial velocities, and the Gaia astrometry. When running MOLUSC, we included all radial velocity data, high-resolution imaging, Gaia astrometry, and limits implied by TOI 1227's CMD position compared to the population (the latter effectively rules out companions <0.3 mag fainter than TOI 1227 at Gaia bandpass and resolution). In total, MOLUSC generated 5 million synthetic companions and determined which could not be ruled out by the data (survivors).
MOLUSC found that only 9% of the generated companions survived (Figure 17). Most of the survivors were faint or stars on wider orbits that happen to have orbital parameters that put them behind or in front of TOI 1227 during all observations. Furthermore, when we restricted the survivors to those that could reproduce the observed transit (unresolved in imaging and sufficiently bright for the observed transit), only 1% of the possible companions remained.
We also used the transit depths measured independently in each of the observed filters to measure the color of any potential undetected companions. At redder wavelengths, the contrast ratio of the primary and companion approaches unity, so a transit or eclipse from the redder companion would appear deeper in the longerwavelength data. Désert et al. (2015) and Tofflemire et al. (2019) showed how to convert the ratio of the transit depths for any two bands (δ 2 /δ 1 ) into the difference in colors between the target and any potential companion: where m 1,x − m 2,x is the color of the primary (p) or companion (c). We excluded TESS from this comparison because of the ambiguity introduced by the dilution term (Sec-  Figure 17. Distribution of surviving companions from MOLUSC in period (P ), eccentricity (e), cosine of the inclination (cos(i)) and mass ratio (q). The distribution represents only 9% of the total generated companions, the rest of which were ruled out by the observational data.
tion 5). We also excluded the LCO g data observed on 2021 Apr 24 (see Section 6.4 for more details). We split the data into four datasets corresponding to the four observed filters; we combined the ASTEP R c data with the SDSS r' data for simplicity, although our coverage in this is poor and the result did not have a significant impact when compared to the effect of other observations taken at other wavelengths. We then fit the dataset for each filter independent of the others. Our fit followed the same method outlined in Section 5, but we placed priors on the wavelength-independent parameters (ρ, b, P , T 0 , and ln(P GP )) drawn from our global fit. As we show in Figure 18, all transits yielded consistent depths. The deepest transits were the two bluest: g and r (although r' depth was very uncertain). A bound eclipsing binary would have yielded a shallower transit at the bluest wavelengths. Taking the 95% confidence of the depth ratio, any companion harboring the transit/eclipse signal must be <0.06 mag redder than TOI 1227 in g − z. To satisfy these color criteria, the companion would have to be approximately equal mass to TOI 1227, making the two appear brighter by >0.3 mag in G. However, TOI 1227 sat within 0.2 mag of the group CMD (Figure 9). More quantitatively, the resulting tight color constraints eliminated all surviving companions from the MOLUSC analysis, ruling out this scenario and validating the signal as planetary.

The unusual April 2021 transit
We excluded the g LCO transit from 2021 Apr 24 in our above analysis. The observations missed ingress but covered most of the transit, all of the egress, and about 1.5 h post-transit. Before color corrections, no clear transit was detected, and even with variability and color corrections the transit was weaker than expected ( Figure 19). While the scatter in the data was large (1.5%), the data suggested a maximum transit depth of 1%. That same transit was observed simultaneously in z s , which showed the expected shape; and in R c , which showed an egress as expected. Although the large errors and partial coverage in the in R c data could not rule out a weaker transit similar to that seen in g'.
Taken alone, the chromatic transit seen in the three data sets from 2021 Apr 24 suggests the transit-like signal is from a background or bound eclipsing binary where the eclipsing system is much redder than the source. However, this conclusion was strongly contradicted by our much more precise g', i', and z s transits from multiple other nights, which we used in our falsepositive analysis (Section 6).  Figure 19. The Apr 2021 g' transit data (grey points) was consistent with a non-detection or at least a significantly weaker transit. The model derived from the fit to all g'band data are shown as a teal shaded region encompassing the uncertainties, and the GP without the transit is shown as a dashed line. The large uncertainties meant that this particular transit had a small weight.
can explain both a non-detection in g' on Apr 24, and the clear detections on other nights, particularly since we have clear detections on both even and odd transits (see Table 1). A more likely explanation is that the 2021 Apr 24 transit was impacted by spots or flares.
Flares from the 24 Myr M dwarf AU Mic have been known to distort, weaken, and nearly erase transits of the two planets, even at TESS wavelengths (Martioli et al. 2021;Gilbert et al. 2021). These effects should be much larger in g , which could explain why z s was not impacted significantly.
If the planet crosses a heavily spotted region, it will block less light than when crossing the warmer region. For such cool stars, small changes in T eff have a large effect on the blue end of the spectrum and a relatively small effect where the SED peaks (in z s ); the result is that modest spots can make a g' transit extremely weak without impacting the z s transit. Such spot crossings can easily happen in one transit and not others, as the rotation period and orbital period are not harmonics of each other and spots on young M dwarfs are known to appear and disappear on month timescales (Rampalli et al. 2021).
To test this, we fit each of the three transits from 2021 Apr 24 individually, placing priors on the period, T 0 , b, ρ, and limb-darkening. As we explain below, the GP kernel can fit out much of the discrepancy (i.e., a flexible GP can get the transits to match), so we did not use a GP model in this test. Instead, we fit a linear trend to the out-of-transit data. This is a conservative assumption, as it will make the errors on the transit depth smaller and hence harder to explain. Fitting for a spot in the transit would be challenging with the data, so we instead fit for the possible range of spots that could explain these depths following the basic methodology of Rackham et al. (2018), but forcing the transit to cross the spotted region rather than a pristine transit cord and a spotted star. There were three free parameters: the spot temperature (T spot ), the fraction of the star that is spotted (f S ), and the transit depth had there been no spots (D). We restricted 5% < f S < 80%, because smaller fractions would not be sufficient to fill the transit cord and larger would be noticeable as a change in T eff with wavelength (see Section 4 and Figure 11). The surface temperature was locked to the value from Section 4. The measurements were the three transit depths and the global fit depth (which only constrains D). We wrapped this in an MCMC framework using emcee, running for 50,000 steps with 50 walkers.
We show the resulting posterior in Figure 20. The spots could not be much cooler than 2700 K unless the spot coverage was 60% because colder spots would change the z s transit. Spot coverage fractions >60% are rare even in young stars (Morris 2020), although not unheard of (e.g., Gully-Santiago et al. 2017). Such high fractions are also unlikely based on the stellar variability and spectrum. On the high end, the spot temperatures were limited by the stellar surface temperature. However, modest (5-30%) spot fractions consistent with ex- The constraints are the three transits (g, Rc, and zs) as well as the overall depth from the global fit. The planet is assumed to be crossing a spotted region, which results in a weaker transit at bluer wavelengths. Although the fit favors a high spot fraction (fS > 50%), this is unlikely given the stellar variability and wavelength-independent temperature. The region of the posterior with modest (< 30%) spot fractions and temperatures similar to the surface (3070 K) can fully explain the anomalous g-band transit and are consistent with all other observables.
pectations for young stars can explain all observations. Since spot temperatures close to the surface temperature are allowed, even moderate ( 30%) spot factions might not impact the spectroscopic analysis significantly (see Section 4 and Figure 11). We conclude that spots could easily explain the anomalous transit.
As an additional test, we re-ran all g' transit data through MISTTBORN. This was effectively a repeat of our analysis from Section 6 including the 2021 Apr 24 data; as in previous fits, a single GP was used to model stellar variability. The goal was to see if including this data could change our interpretation. The resulting transit depth was 2.27 +0.33 −0.26 %, consistent with our earlier fits as well as the other g-band data (Figure 18). The depth change was insignificant because the SOAR g'-band data are far more precise, and significant stellar variability can explain a lot of the apparent discrepancy (which the GP can fit; see Figure 19). The updated depth was still sufficient to rule out all possible binaries from the MOLUSC analysis.  Berta et al. 2013), but such programs found only upper limits in this regime. The planet experiences stellar irradiance compared to most known planets, and should not be inflated. Instead, by >100 Myr we expect TOI 1227 b to more closely resemble the common sub-Neptunes seen around older stars. We explore this further below.

Evolutionary Modeling
To determine the possible current internal structure and future evolution of TOI 1227 b we compared its current properties to evolutionary models. We followed the framework of Owen (2020) where we constrained TOI 1227 b's core mass, initial hydrogen dominated atmosphere mass fraction and initial (Kelvin-Helmholtz) cooling timescale to match its current age and radius. We considered planets with core masses <25 M ⊕ and restricted ourselves to initial cooling timescales <1 Gyr and applied uniform priors in core-mass, log initial envelope mass fraction and log initial cooling timescale. The evolutionary tracks were computed using mesa (Paxton et al. 2011(Paxton et al. , 2013 and included mass-loss due to photoevaporation and irradiation from an evolving host star (Owen & Wu 2013;Owen & Morton 2016). The methodology also accounted for uncertainty in the protoplanetary disk dispersal timescale. Figure 22 shows the resulting marginalized probability distribution for core-mass and initial envelope mass fraction. Normalised Probability Figure 22. The marginalized probability distribution for TOI 1227 b's core mass and initial hydrogen dominated envelope mass fraction that explains its current properties.
With initial envelope mass fractions <1, our results confirmed that TOI 1227 b is likely a young inflated planet that is actively cooling, eventually becoming one of the ubiquitous sub-Neptunes. We note the slight preference for higher core mass is likely artificial because the models with higher core masses are strongly biased towards longer initial cooling times. While we selected an upper limit of 1 Gyr in our modeling, this was likely an overestimate as the longest initial cooling timescales expected predicted by formation models that include an early "boil-off" phase are of order a few 100 Myr. However, our modeling indicated that core masses 5 M ⊕ are strongly ruled out.
We followed the future evolutionary pathways of a subset of the model planets by selecting three representative cooling tracks and appropriate initial cooling timescales. These results are presented in Figure 23, showing that for plausible initial planet masses and initial cooling timescales they evolve into a sub-Neptune with a size 3-5 R ⊕ at billion-year ages. Such planets are still rare around M dwarfs, but much more consistent with the locus of detections ( Figure 21). on LCC has identified substructures, i.e. smaller populations, each with slightly different Galactic motion, position, and age. TOI 1227 b was flagged as a member of the A0 group by Goldman et al. (2018) and LCC-B by Kerr et al. (2021), which are effectively the same group. Because of the contradictory and easily forgotten naming, we denote this group Musca after the constellation containing most of the members. As part of our effort to better measure the age of TOI 1227 b, we used rotation, lithium, and stellar evolution models to determine an age of 11±2 Myr for Musca, and hence TOI 1227 b. The young age placed TOI 1227 b as one of only two systems with transiting planets younger than 15 Myr (K2-33 b; Mann et al. 2016a;David et al. 2016), and one of only eleven transiting planets younger than 100 Myr (V1298 Tau, DS Tuc, AU Mic, TOI 837, HIP 67522, and TOI 942 David et al. 2019a;Newton et al. 2019;Benatti et al. 2019;Plavchan et al. 2020;Rizzuto et al. 2020;Bouma et al. 2020;Zhou et al. 2020;Carleo et al. 2021). Even including radial velocity detections only adds a few (Johns-Krull et al. 2016;Yu et al. 2017;Donati et al. 2017), some of which are still considered candidates (Donati et al. 2020;Damasso et al. 2020).
Like K2-33 b, TOI 1227 b is inflated, with a radius much larger than known transiting planets orbiting similar-mass stars ( Figure 21). TOI 1227 b's radius is closer to Jupiter than to the more common 1-3 R ⊕ planets seen around mid-to-late M dwarfs (Berta et al. 2013;Hardegree-Ullman et al. 2019). Such a large planet would be obvious in surveys of older M dwarfs, especially given that the equivalent mature M dwarfs are both less variable and a factor of a few smaller than their premain-sequence counterparts. Thus, TOI 1227 b is evidence of significant radius evolution, especially when combined with similar earlier discoveries.
Interestingly, evolutionary models favor a final radius for TOI 1227 b of >3 R ⊕ , which would still be abnormally large (Figure 21) for the host star's mass. However, that discrepancy is easier to explain as observational bias. Furthermore, the evolutionary tracks are poorly constrained by the existing young planet population, especially given that AU Mic b is the only <50 Myr planet with both mass and radius measurements (Klein et al. 2021;Cale et al. 2021). The lack of masses remains a significant challenge for these kinds of comparisons (Owen 2020). TOI 1227 b is too faint for radial velocity follow-up in the optical, but might be within the reach of high-precision NIR monitoring. Our IGRINS monitoring achieved 40-50 m s −1 precision, which is below the level from the stellar jitter. More dedicated monitoring and careful fitting of stellar and planetary signals should be able to detect a planet below 40 M ⊕ .
A less challenging follow-up would be spin-orbit alignment through Rossiter-Mclaughlin, where the signal is expected to be >200 m s −1 and the transit timescale ( 5 h) is much shorter than the rotational jitter (1.65 days). Recent measurements of stellar obliquity in young planetary systems have found that they are generally aligned with their host stars (e.g., Rodriguez et al. 2019;Zhou et al. 2020;Wirth et al. 2021), suggesting that these planets formed in situ or migrated through their protoplanetary disk. However, the number of misaligned systems even around older stars is 10% (varying with spectral type and planet mass; Fabrycky & Winn 2009; Morton & Winn 2014;Campante et al. 2016). We likely need a larger sample of obliquity measurements in young systems for a significant statistical comparison.
Follow-up and characterization of TOI 1227 b highlight the difficulties of characterizing and validating such young planets. We initially mischaracterized both the host star and planet, and ground-based transits yielded some contradictory results. The host star (TOI 1227) was initially thought to be part of the 5 Myr Cha cluster, based on the BANYAN code. At that age, the stellar host would be even larger, yielding an even larger radius for the planet and creating a more complex set of false-positive scenarios (e.g., the planet radius would be more consistent with a young brown dwarf). Without Gaia, it is unlikely that Musca would have been recognized, and assigning TOI 1227 to a given population without a parallax would have been extremely challenging.
Similarly, the THYME team originally flagged TOI 1227 b as a likely eclipsing binary based on its V-shaped transit and large depth. Only with the SALT/HRS spectra and a consistent transit depth from Goodman/SOAR did we pursue further follow-up. This motivated both spectroscopic monitoring with IGRINS and the suite of multi-band transit photometry observed over 1.5 years to better characterize the planet. While the majority of the transit photometry painted a clear picture that the signal was from a planet, a single transit showed an inconsistent transit depth, suggesting a false positive. We argued in Section 6.4 that regular flaring, the planet occulting spots, or simply low-precision photometry, could easily explain the unusual transit. Spots could also explain why our other g' transits were deeper than the i and z s data ( Figure 18); if the transit cord has fewer spots than the rest of the stellar surface, the bluest transits appear deeper (this was seen in K2-25; Thao et al. 2020).
While the case for TOI 1227 b as a planet is strong, the initially confusing results raise concerns about the future follow-up of young planets. Radial velocity detection remains challenging. As with K2-25 b, the discovery transit for TOI 1227 b was V-shaped, making it hard to classify on transit depth and shape alone. Transits from young planets are likely to show some chromaticity even if the planet is real given the presence of spots, and smaller planets are not as amenable to the large suite of ground-based follow-up used here. Earlier studies often made use of Spitzer, which provided the advantage of a wide wavelength range while operating in a regime where the effects of spots and flares are significantly smaller. With the end of Spitzer, we may need to focus on improving the NIR photometric precision from the ground. For now, we encourage caution when rejecting young planets based on metrics tuned for older systems.

ACKNOWLEDGMENTS
The authors thank the anonymous referee for providing insightful comments and suggestions on the paper. AWM would like to thank Bandit, who sat directly on AWM's keyboard whenever MISTTBORN was running, preventing him from working on this manuscript too much. The THYME collaboration also wants to thank Halee, Wally, Maizie, Dudley, Charlie, Marley, and Edmund for their thoughtful discussions and emotional support. AWM was supported through NASA's As-         Ages derived from direct comparison to theoretical evolutionary models are impacted by a wide range of systematics. Binaries, for example, make the association look younger; an effect that varies with stellar mass (Sullivan & Kraus 2021). The effect can be mitigated by explicitly including binaries in their model likelihood assuming an initial mass function and binary fraction (Kerr et al. 2021). This is less effective if the sample is biased towards or against binaries, which could easily arise from Gaia-based astrometry until later releases that account for binary motion. For young associations like those considered here, field interlopers (non-members) also bias the result towards older ages. Specific to populations like Musca is the presence of interlopers from other populations in LCC, which could bias the age in a direction that depends on the nearby populations.
Our solution is to fit the data with two populations simultaneously (a mixture model). The first population is the (ideally) single-age and single-star population of interest. The second is the 'outlier' population, which may contain nonmembers (field or nearby populations), binaries, stars with inaccurate photometry or parallaxes, or simply regions of the CMD that are poorly captured by the model grid. Following Hogg et al. (2010), the likelihood would be a sum of the two populations weighted by an amplitude: where M x,i is the absolute magnitude of the ith star in band x, f (age, E(B − V ), col i ) is the predicted absolute magnitude from the model isochrones for a given age, reddening (E(B − V )), and color (col i ), and σ Mx,i is the total error in the star's absolute magnitude and the propagated error from the star's color. The parameters P b , V b , and Y b are free parameters that describe the second population: the amplitude of the outliers (or the fraction of the population that are outliers), the mean absolute magnitude of the outliers, and the variance in the absolute magnitude of the outliers, respectively. The simple likelihood above works well when most outliers are field M dwarfs, which are well described by a single Gaussian distribution in absolute magnitude versus color. However, we have both field and young interlopers that may occupy a large range of the CMD. Fortunately, such stars are still (mostly) restricted to the physically allowed regions of the CMD. Thus, we change the second term in the likelihood from a simple distribution in absolute magnitude to an offset from the population of interest: A downside of using an offset from the CMD to model the outliers is that the three outlier parameters can be biased by targets in regions of parameter space that are not well described by a simple offset, such as sources with poor parallaxes or M-dwarf white-dwarf binaries. This can be fixed by applying initial quality cuts. A more complicated problem is that the offset from the expected CMD likely varies with luminosity/mass, while we assumed a single parameter.
As a test, we ran our mixture model on the LCC-C (Crux-S) group and LCC-B groups from Kerr et al. (2021), which are similar cases to Musca. In both cases, we used the PARSECv1.2S isochrones (Bressan et al. 2012), assumed Solar metallicity, and cut out stars with RUWE>1.2. We restricted the comparison to Gaia magnitudes, as those were available for all stars, but tests adding in JHK magnitudes yielded similar results. We show the resulting CMD and fit in Figure 24. For Crux-S, we derived an age of 12.6 +1.6 −0.8 Myr, consistent with the 14.6±0.8 Myr from Kerr et al. (2021). For LCC-B, we estimated an age of 11.2 +1.1 −0.8 Myr, consistent with 13.0±1.4 from Kerr et al. (2021). We also performed similar tests with older associations, such as MELANGE-1 from Tofflemire et al. (2021), which gave us an age of 210±70 Myr, an excellent match to the Li-and gyro-based ages (250 +50 −70 Myr) from Tofflemire et al. (2021).  Kerr et al. (2021). Points including in our fit are labeled as blue circles, with blue squares representing stars excluded (mostly because of a high RUWE). Each point is shaded by the probability of being an outlier (i.e., not well-described by the model). This may include non-members, nearby association interlopers, binaries, or simply stars with poor parallaxes/photometry.
One advantage of the mixture model approach is the handling of regions where the models perform poorly. For example, in both test cases the PARSEC models fail to reproduce the population of M0-M2 dwarfs ( Figure 24). Given that the magnetic models from DSEP better reproduce this region, it is likely that the offset is due to errors in the models as opposed to a high rate of binarity in that part of the selected population. In the case of LCC-C, there are few sources in that regime, so they have little effect on the overall age. For Crux-S, the mixture model weights these down in the fit.
However, the ability of this method to discard points is also a problem with this method; for Crux-S, the fit identifies many of the high-mass stars as medium-probability outliers, likely because it is challenging to fit the low-mass stars simultaneous with the high-mass ones (Feiden 2016). A solution to this would be to run some systems with independent age estimates (such as from lithium depletion or traceback ages, Soderblom et al. 2014) to identify such systematics and apply corrections or weights. However, this requires a much more extensive set of associations with non-isochronal ages. A less elegant approach is to run the mixture code with multiple model grids.
Another drawback is when the population is poorly described by two models. However, the framework can be increased to explicitly model additional populations with sufficient terms. For example, one could use one term to handle the population of interest, a term to handle binaries with a prior or limits on the Y b -like term, and a term to model the (assumed) single-age population of young interlopers from other LCC groups. We found the simpler model to be sufficient for our purposes here, but hope to explore this in future iterations.