TESS Hunt for Young and Maturing Exoplanets (THYME). VI. An 11 Myr Giant Planet Transiting a Very-low-mass Star in Lower Centaurus Crux

TOI 1227 was observed by the TESS mission (Ricker et al. 2014) from UT 2019 April 22 through 2019 June 19 (Sectors 11 and 12) and then again from 2021 April 28 through 2021 May 26 (Sector 38). In all three sectors, the target fell on Camera 3. The first two sectors had 2-minute 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 TESS object of interest (TOI) at that phase (G03141; PI Newton). We initially used the 2-minute cadence for all analysis for computational efficiency and TESS cosmic-ray mitigation ³³ (not included in 20 s data), but we 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; Smith et al. 2012; Stumpe et al. 2012, 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). ³⁴

The planet signal, TOI 1227.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). ³⁵ 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 67522b (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 three 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.

Zoom In Zoom Out Reset image size

The photometry is summarized in Table 1. We provide details by instrument below.

Table 1. Log of Transit Observations

Start	Telescope	Filter	Transit	T ${}_{\exp }$	Obs Duration
(UT)			#	(s)	(hr)
2019 Apr 22 a , b	TESS S11	TESS	1	120	N/A
2019 Apr 22 a , b	TESS S12	TESS	2	120	N/A
2020 Jan 15 b	SOAR	i'	9	120	8.1
2020 Jan 15	LCO	i'	9	200	2.8
2020 May 03 b	LCO	i'	13	200	1.9
2020 Jul 24	LCO	r'	16	200	2.5
2020 Jul 24 b	LCO	zs	16	200	2.5
2021 Mar 28 b	SOAR	g'	25	120	7.6
2021 Apr 24	ASTEP	Rc'	26	200	5.7
2021 Apr 24	LCO	g'	26	300	6.2
2021 Apr 24 b	LCO	zs	26	210	6.2
2021 Apr 28 a , b	TESS S38	TESS	27	20	N/A
2021 Jun 17 b	LCO	g'	28	300	4.1

Notes.

^a TESS observed TOI 1227 for one transit in each of Sectors 11, 12, and 38. ^b Included in the global fit.

Download table as:  ASCIITypeset image

2.3.1. Goodman/SOAR

2.3.2. LCO Photometry

2.3.3. ASTEP Photometry

Since atmospheric extinction is strongly color dependent, changes in air mass 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 air mass. However, TOI 1227 was unusual because of both its long orbital period (27.3 days) and its red color; of stars <5' from the target and differing by <1 mag 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 second-order 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). ³⁷ We then combined the relevant BT-SETTL models (assuming solar metallicity) with an air-mass-dependent 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 data sets.

2.5.1. Goodman/SOAR

We observed TOI 1227 with the Goodman spectrograph (Clemens et al. 2004) on the Southern Astrophysical Research (SOAR) 4.1 m telescope located at Cerro Pachón, Chile, on two nights. We used the red camera, the 1200 line mm^–1 grating in the M5 setup, and the 046 slit rotated to the parallactic angle. This setup gave a resolution of R ≃ 5880 spanning 6150–7500 Å.

We obtained the first spectrum on 2019 December 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 February 23 (UT) under clear conditions. Our goal was a spectrum with a signal-to-noise ratio (S/N) ≫ 100 pixel^–1 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.

Zoom In Zoom Out Reset image size

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 2021 February 23 and 24 April 2021 using an identical setup to the observations for TOI 1227. Exposure times were set to ensure an S/N greater than 50 (per pixel) around the Li line and varied from 90 to 420 s per exposure, with at least five 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 fifth-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%.

2.5.2. HRS/SALT

2.5.3. IGRINS/Gemini-S

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 ∼002 in ideal conditions with a field of view of about 12. TOI 1227 was observed in 17 sets of 1000 × 60 ms exposures by the 'Alopeke-Zorro visiting instrument team with the standard speckle imaging mode in the narrowband 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).

Zoom In Zoom Out Reset image size

The limiting contrasts for binary companions are set by the redder [832] filter, ruling out equal-brightness companions at ρ = 003 and excluding companions fainter than the target star by Δ[832] < 1.5 mag at ρ = 005, Δ[832] < 4.5 mag at ρ = 010, and Δ[832] < 5.0 mag at ρ = 02. Converting M_[832] to M_G using the color–magnitude relations of Kraus et al. (in preparation), 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 > 50 M_J at ρ > 5 au, M > 15 M_J at ρ > 10 au, and M > 14 M_J at ρ > 20 au.

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}_{\mathrm{lim}}\lt 21$ mag, which at an age of ≃10 Myr approximately corresponds to a limit of ${M}_{\mathrm{lim}}\simeq 14\,{M}_{J}$. 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).

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.6 m telescope (La Silla Observatory), FEROS at the 2.2 m telescope (La Silla Observatory), and/or X-shooter at VLT (Cerro Paranal Observatory).

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; Lindegren et al. 2021; Gaia Collaboration et al. 2021; Riello et al. 2021). From EDR3 we also retrieved the renormalized unit weight error (RUWE) for all stars. The RUWE value is effectively an astrometric reduced χ² value, normalized to correct for color- and brightness-dependent effects. ⁴² RUWE should be around 1 for well-behaved sources, and higher values (RUWE ≳ 1.3) suggest the presence of a stellar companion (Ziegler et al. 1919; 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.

To aid with stellar characterization of TOI 1227, we also pulled photometry from the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006), the Wide-field Infrared Survey Explorer (WISE; Cutri et al. 2014), and the SkyMapper Second Data Release (SkyMapper DR2; Onken et al. 2019).

To identify any comoving and coeval population to TOI 1227, we first used the BANYAN-Σ tool (Gagné et al. 2018), ⁴³ 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; see also Figure 4). The BANYAN algorithm preferred placing Cha over LCC likely because of better agreement in UVW despite poor XYZ 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 UVW model for LCC and hence the membership probabilities.

Zoom In Zoom Out Reset image size

TOI 1227 was more recently identified as a member of the "B" subpopulation of LCC by Kerr et al. (2021) and the A0 subgroup 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 subgroups. 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 XYZ positions, proper motions, and available RVs 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 ⁴⁴ (Tofflemire et al. 2021). FriendFinder uses the Gaia EDR3 astrometry and RV of a source (we used our IGRINS systemic velocity) to compute Galactic UVW and XYZ and then projects the UVW 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⁻¹ from the expected values calculated by FriendFinder. These cuts were based on the fact that 1 km s⁻¹ ≃ 1 pc Myr⁻¹ and an estimated age for the population of ≃10 Myr. We estimated that 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 work—have 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 farther away from TOI 1227 because of our separation cut.

Zoom In Zoom Out Reset image size

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, based on both Galactic position (Figure 4) and proper motion (Figure 5). The mean velocity of the members with literature RVs (12.8 km s⁻¹) is also in good agreement with the value for TOI 1227. A randomly selected star that matches Musca in XYZ has a ≪1% probability of matching Musca in proper motion and RV 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.

Zoom In Zoom Out Reset image size

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 and 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 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. (1919). TOI 1227's rotation period (1.65 days) was consistent with the Musca sequence. This effectively sets a limit of ≲60 Myr for the group and further confirms TOI 1227 as a member.

Zoom In Zoom Out Reset image size

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; see also Rizzuto et al. 2015). The location and depth of the lithium dip also depend strongly on age.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For our Li determinations, we measured the equivalent width of the Li λ6708 line for 22 stars using our high-resolution 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 equal mass 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 8(b) 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).

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 appendix. 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 that model differences are comparable to our measurement errors.

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 1227b as well, as any age spreads are likely to be similar to, or smaller than, the adopted uncertainties.

We first fit the SED following the methodology from Mann et al. (2016b) and highlighted in Figure 10. To briefly summarize, we compared the observed photometry and Goodman spectrum to a grid of template spectra drawn from nearby (likely reddening-free) young moving groups. The observed data were simultaneously fit with solar-metallicity BT-SETTL CIFIST models (Baraffe et al. 2015). The atmospheric models were used both to estimate T_eff and to 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 mag to the errors of all optical photometry. The resulting fit yielded T_eff = 3072 ± 84 K, F_bol = (7.87 ± 0.53) × 10⁻¹¹ erg cm⁻² s⁻¹, L_* = (2.51 ± 0.17) × 10⁻³ L_⊙, ${A}_{V}={0.21}_{-0.09}^{+0.11}$ 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.

Zoom In Zoom Out Reset image size

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 parallax, provided an additional estimate of R_*. The technique is similar to the infrared flux method (IRFM; Blackwell & Shallis 1977). IRFM-based radii of low-mass 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 offer another advantage here; T_eff determined from the low-resolution Goodman spectra was driven by the molecular bands, while the high-resolution 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-S/N IGRINS spectra obtained. We fit each order separately but only included orders with low telluric contamination (determined by eye) and S/N > 80. For each fit, we explored six free parameters: T_eff, $\mathrm{log}\,g$, RV, 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 Markov Chain Monte Carlo 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.

Zoom In Zoom Out Reset image size

We repeated our fit to the IGRINS data after co-adding all 11 spectra. Because the spectra were taken across the rotational phase, these likely sample different spot coverage patterns. However, the resulting spectrum had significantly higher S/N (>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-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.01 M_⊙) than expected errors. To redden the model photometry, we used synphot (Lim 1919), 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.015 M_⊙). 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.

Canto Martins et al. (1919) 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 RVs from IGRINS/Gemini and HRS/SALT spectra. The mean value from the IGRINS data was 16.65 ± 0.24 km s⁻¹, while the SALT data yielded a marginally inconsistent value of 17.8 ± 0.3 km s⁻¹. We adopted the former value, as the IGRINS data have significantly higher S/N and are 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.

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 NGC 4372 (identified 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 identifies the X-ray source as 2RXP J122703.8−722702, situated 49 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⁻³ counts s⁻¹ 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 1993 September 4–6. Using the energy conversion factor equation of Fleming et al. (1995), this count rate and HR1 hardness ratio translate to a soft X-ray energy flux of f_X = 1.70 ± 10⁻¹⁴ erg s⁻¹ cm⁻², which at the distance to TOI 1227 translates to a soft X-ray luminosity of L_X = 10^28.32 erg s⁻¹.

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 compared the 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 days). 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⁻¹), likely due to higher stellar variability in the optical.

A Jovian-mass planet would induce an RV signal of ≃220 m s⁻¹, while the scatter in the IGRINS velocities was ≃130 m s⁻¹ (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 RV 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.

Zoom In Zoom Out Reset image size

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 mag of TOI 1227. Since the depth is consistent across all filters, this constraint applies to all griz magnitudes. The multiband 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₁₂) to the time between first and third contact (T₁₃). However, because TOI 1227b 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). Between the blue POSS images (1976.3) and the most recent SOAR transit data (2021.2), TOI 1227 moved more than 18. This motion was much larger than the resolution of our SOAR transit imaging ( ≃ 08) and speckle imaging (009) and somewhat larger than the resolution of the POSS plates (17). 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 point-spread function from StarFinder (Diolaiti et al. 2000), we were able to rule out any source <3 mag fainter than TOI 1227. Additionally, both USNO-A2.0 and USNO-B reported no other nearby source data (Monet et al. 2003) but were sensitive to stars >3 mag fainter than TOI 1227 (in R).

To check for a bound companion eclipsing binary, we first ran the MOLUSC code (Wood et al. 2021). MOLUSC simulated binary companions following an empirically motivated random distribution in binary parameters. We limited the mass of the synthetic companions to between 10 M_J and 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, RVs, and the Gaia astrometry. When running MOLUSC, we included all RV 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.

Zoom In Zoom Out Reset image size

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 longer-wavelength 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 (δ₂/δ₁) into the difference in colors between the target and any potential companion:

$$\begin{eqnarray}&&{m}_{1,c}-{m}_{2,c}\lt {m}_{1,p}-{m}_{2,p}+2.5\mathrm{log}\left(\displaystyle \frac{{\delta }_{2}}{{\delta }_{2}}\right),\end{eqnarray} \tag{ 2 }$$

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 (Section 5). We also excluded the LCO g' data observed on 2021 April 24 (see Section 6.4 for more details). We split the data into four data sets 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 data set 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₀, and $\mathrm{ln}({P}_{\mathrm{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.

Zoom In Zoom Out Reset image size

We excluded the g' LCO transit from 2021 April 24 in our above analysis. The observations missed ingress but covered most of the transit, all of the egress, and about 1.5 hr 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. However, the large errors and partial coverage in the R_c data could not rule out a weaker transit similar to that seen in g'.

Zoom In Zoom Out Reset image size

Taken alone, the chromatic transit seen in the three data sets from 2021 April 24 suggests that 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 false-positive analysis (Section 6). No false-positive scenario can explain both a nondetection in g' on April 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 April 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 (Gilbert et al. 2021; Martioli 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 April 24 individually, placing priors on the period, T₀, 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 ones 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 expectations 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.

Zoom In Zoom Out Reset image size

As an additional test, we reran all g' transit data through MISTTBORN. This was effectively a repeat of our analysis from Section 6 including the 2021 April 24 data; as in previous fits, a single GP was used to model stellar variability. The goal was to see whether including these data could change our interpretation. The resulting transit depth was ${2.27}_{-0.26}^{+0.33} \%$, 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.

To determine the possible current internal structure and future evolution of TOI 1227b, we compared its current properties to evolutionary models. We followed the framework of Owen (2020), where we constrained TOI 1227b'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, 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.

Zoom In Zoom Out Reset image size

With initial envelope mass fractions <1, our results confirmed that TOI 1227b is likely a young inflated planet that is actively cooling, eventually becoming one of the ubiquitous sub-Neptunes. We note that the slight preference for higher core mass is likely artificial because the models with higher core masses are strongly biased toward 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 hundred 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).

Zoom In Zoom Out Reset image size