TOI-4201: An Early M Dwarf Hosting a Massive Transiting Jupiter Stretching Theories of Core Accretion

We confirm TOI-4201 b as a transiting Jovian-mass planet orbiting an early M dwarf discovered by the Transiting Exoplanet Survey Satellite. Using ground-based photometry and precise radial velocities from NEID and the Planet Finder Spectrograph, we measure a planet mass of 2.59−0.06+0.07 M J, making this one of the most massive planets transiting an M dwarf. The planet is ∼0.4% of the mass of its 0.63 M ⊙ host and may have a heavy-element mass comparable to the total dust mass contained in a typical class II disk. TOI-4201 b stretches our understanding of core accretion during the protoplanetary phase and the disk mass budget, necessitating giant planet formation to take place either much earlier in the disk lifetime or perhaps through alternative mechanisms like gravitational instability.


INTRODUCTION
The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) observes millions of stars across the entire sky searching for transiting candidates.It has been instrumental in finding planets around M-dwarf stars, the most common type of star in the Galaxy.Among the confirmed TESS candidate planets orbiting M-dwarfs are ∼15 giant planets in short period orbits (e.g., Gan et al. 2022;Kanodia et al. 2023;Hobson et al. 2023;Kagetani et al. 2023), which seem to challenge current theories of planet formation through core accretion (Laughlin et al. 2004;Ida & Lin 2005).
Under the core-accretion model, the formation of giant exoplanets around M-dwarf stars (GEMS) is difficult for two main reasons.First, we expect disk mass, particularly the solid surface density, to scale with host star mass (Andrews et al. 2013;Pascucci et al. 2016).Under typical opacity assumptions runaway gas accretion can only be triggered after a core of ∼ 10 M ⊕ is formed.Since 10 M ⊕ is a large percentage of the solid mass typically available in a disk around an M-dwarf (Ansdell et al. 2016;Tazzari et al. 2021;Manara et al. 2022), forming this core would be a challenge.Second, due to their lower host star masses, the Keplerian orbital time-scales are much longer for M-dwarfs compared to solar-type stars, which coupled with their lower disk surface densities should make it harder to form a massive enough core in time to initiate runaway gaseous accretion before the disk dissipates (Laughlin et al. 2004).
The discovery and characterization of these GEMS is already challenging our understanding of giant planet formation, and despite their enhanced detection probabilities with transit and radial velocity observations (large planet to stellar mass and radius ratios), they remain rare.While currently difficult to perform sample level comparisons on the existing population (< 15), efforts are underway to significantly increase this sample size and to enable such analysis in the future.
In this letter, we present the confirmation of a massive super-Jupiter around the early M-dwarf TOI-4201.While previous gas giants have been difficult to explain from a mass budget perspective (Kanodia et al. 2023;Hobson et al. 2023), one proposed explanation is that planet formation begins early on in disk lifetime (< 1 Myr), when the disks are much more massive.The discovery of a super-Jupiter around an early M-dwarf with † NASA Postdoctoral Fellow ‡ Flatiron Research Fellow a planet to host mass ratio of ∼0.4% stretches this even further.
To characterize TOI-4201 and confirm the TESS signal as a planet, we use precision radial velocities (RVs) from the Planet Finder Spectrograph (PFS) and NEID, high contrast speckle imaging from WIYN/NESSI, and ground based photometry from TMMT and LCRO, as well as archival photometry from LCOGT.We detail these observations in Section 2 and describe the stellar characterization in Section 3. Section 4 covers the data analysis, including details of the joint modeling of RVs and photometry.In Section 5, we compare TOI-4201 b with other giant planets around M-dwarfs and consider potential formation scenarios.We summarize our findings in Section 6.

TMMT & LCRO
We observed a transit of TOI-4201 b on the night of 2021 December 28 with the Three-hundred MilliMeter Telescope (TMMT; Monson et al. 2017) at Las Campanas Observatory, with 120 s exposures in Bessell I. Observations were taken as the target set from an airmass of 1.19 to 2.48.
We then observed a second transit on the night of 2022 January 15, with TMMT in Bessell V and the 0.6 m Las Campanas Remote Observatory (LCRO) in SDSS i , both using 300 s exposures.The target set from an airmass of 1.04 to 1.66 over the course of the night.
We reduced the data from all three observations (bias correction, flat-fielding, cosmic/bad-pixel removal, etc.) following the procedure outlined in Monson et al. (2017).Then, we performed differential aperture photometry on all images using a python script based on Monson et al. (2017).Final light curves are plotted in Figure 1.

1.0 m LCOGT Archival Data
We also retrieve a full transit for TOI-4201 b from the Las Cumbres Observatory global telescope network (LCOGT; Brown et al. 2013) public data archive.This transit was observed in both SDSS g and i on the night of 2021 October 3 (proposal ID: KEY2020B-005, PI: Shporer, A.) from the Sinistro imaging camera on the 1m LCOGT Dome B telescope at Cerro Tololo Inter-American Observatory.These observations were taken mildly defocused (FWHM ∼2.5 ) at exposure times of 300 s in g and 180 s in i , with the target rising from an airmass of ∼2.3 to 1.06 over the course of the night.The raw data were automatically processed by the BANZAI pipeline (McCully et al. 2018).We then perform aperture photometry on the processed images using AstroImageJ (AIJ; Collins et al. 2017).We see a strong slope in both the g and i lightcurves, so we detrend both datasets, in time, prior to our analysis.These transits are shown in Panel e and f of Figure 1.

NESSI
We observed TOI-4201 on the night of 2022 April 18 with the NN-Explore Exoplanet Stellar Speckle Imager (NESSI; Scott et al. 2018) on the WIYN1 3.5 m telescope at Kitt Peak National Observatory to rule out stellar companions.We took diffraction-limited exposures using the red camera and the Sloan z filter at a 40 ms cadence for 9 minutes and reconstructed the speckle image following the methods described by Howell et al. (2011).We compute 5σ contrast limits as a function of separation, ∆θ, from the primary source and find no evidence for nearby or background sources with ∆z < 3.5 at ∆θ = 0.5" and ∆z < 3.9 at ∆θ = 1.2".

NEID
We observed TOI-4201 using NEID, a fiber-fed, environmentally stabilized spectrograph on the WIYN 3.5 m telescope (Halverson et al. 2016;Schwab et al. 2016;Robertson et al. 2019) for three epochs between 2022 November 15 and 2022 November 20.NEID covers the wavelengths 380 to 930 nm with a spectral resolution of R ∼ 110, 000.Each visit consisted of a single 30 minute exposure.We use the standard NEID data reduction pipeline followed by the SERVAL algorithm developed by Zechmeister et al. (2018), and adapted for NEID by Stefànsson et al. (2022).At 850 nm, the SNR for each of the three RV points is 6.0, 7.6, and 5.5; we attribute the high jitter seen in the fit to the low SNR.The final NEID RVs are included in Table 1 and as a machine readable table with the manuscript.

Planet Finder Spectrograph
We observed TOI-4201 with the Planet Finder Spectrograph (PFS; Crane et al. 2006Crane et al. , 2008Crane et al. , 2010) ) on the 6.5 m Magellan II (Clay) telescope at Las Campanas Observatory.Between 2022 November 8 and 2023 April 12 we obtained 6 visits, each consisting of three exposures of 1200 seconds 2 in 3 x 3 CCD binning mode with a 0.3 slit.This data was taken with the iodine gas absorption cell in the light path, which imprints a dense forest of molecular lines (Hatzes 2019) between 5000 to 6200 Å.We also obtained one template spectrum without the iodine cell consisting of six exposures of 1200 seconds.The RVs were derived following the methodology of Butler et al. (1996).As noted by Hartman et al. (2015) and Bakos et al. (2020), due to the faintness of the target, and the optical region for the iodine region, the formal errors on the PFS RVs are likely underestimated.The final PFS RVs are included in Table 1 and as a machine readable table with the manuscript.

STELLAR PARAMETERS
The stellar parameters presented in Table 2 are derived using the available broadband photometry and Gaia astrometry.The stellar metallicity is estimated as [Fe/H] = 0.30 ± 0.15 using WISE and Gaia colors3 (Equation 4from Duque-Arribas et al. 2023).This photometric relationship was determined by Duque-Arribas et al. ( 2023) to be the most accurate photometric calibration when compared to a well-characterized (i) spectroscopic sample of M-dwarfs (Birky et al. 2020) and (ii) M-dwarfs in binary systems (Montes et al. 2018).The stellar radius is calculated as R = 0.62 ± 0.02 R using Equation 5from Mann et al. (2015).This value was then used to determine a mass of M = 0.63 ± 0.02 M with Equation 6from Schweitzer et al. (2019), which also agrees with the stellar mass estimate from Mann et al. (2019).The stellar effective temperature is derived as 3920 ± 50 K with the empirical calibration from Equation 7 in Rabus et al. (2019), that was derived using interferometric observations of well-characterized M-dwarfs.The adopted stellar parameters are consistent at the 1σ level with the (i) Bayesian stellar parameters ([Fe/H] = 0.300 +0.002 −0.388 , M = 0.65 +0.002 −0.050 M , T eff = 3955 +20 −49 K) derived using the StarHorse code and a combination of multi-wavelength photometry and Gaia parallaxes (Anders et al. 2022) and (ii) stellar parameters derived from a fit to the spectral energy distribution using EXOFASTv2 (using the broadband photometry in Table 2; Eastman et al. 2019

Galactic kinematics
We calculate the UVW velocities in the barycentric frame and relative to the local standard of rest from Schönrich et al. (2010) using GALPY (Bovy 2015)4 .These velocities are reported in Table 2 and are used to clas-sify TOI-4201 as a thin disk star using the criteria from Bensby et al. (2014).TOI-4201 is also determined to be a field star (> 99% change of membership) using the BANYAN Σ tool (Gagné et al. 2018).

JOINT FITTING OF PHOTOMETRY AND RVS
We perform a joint fit of the radial velocity time series and photometry using exoplanet (Foreman-Mackey et al. 2021a), which utilizes PyMC3 to perform Hamiltonian Monte Carlo (HMC) posterior sampling algorithm (Salvatier et al. 2016).We follow a prescription similar to Kanodia et al. (2023), where we allow eccentricity to float, and include a Gaussian Process (GP) kernel for the TESS photometry to detrend out instrument systematics.Similar to previous works, we also include a dilution term to correct the TESS photometry based on the ground-based transits.
We performed joint fitting of the RV timeseries and the photometric curves using exoplanet (Foreman-Mackey et al. 2021a), a software package that utilizes PyMC3, a Hamiltonian Monte Carlo (HMC) posterior sampling algorithm (Salvatier et al. 2016).HMC is a Markov chain Monte Carlo method that uses first-order gradients to avoid random walk behavior, that is further improved by the implementation of the No U Turn Sampler (NUTS) which reduces sensitivity to user specified parameters (Hoffman & Gelman 2014).We use a Keplerian orbit to model the RVs, leaving the eccentricity as a free parameter.A linear trend and RV offsets for each instrument were fit to the data to account for longterm changes from both astrophysical and instrumental causes.
We modeled the transits using starry (Luger et al. 2019;Agol et al. 2020), which uses the analytic formulae derived in Mandel & Agol (2002) to compute the light curves based on model stellar atmospheres and relies on separate quadratic limb darkening parameters for each filter.During the joint fit with all RV and photometric instruments, we fit each phased transit with separate limb-darkening coefficients.A jitter term was fit for each data set as a white noise model that was then added in quadrature to the uncertainty of the data sets.We used celerite2 (Foreman-Mackey et al. 2017; Foreman-Mackey 2018) to include a Gaussian Process (GP) kernel in the likelihood function of the TESS data to model the quasi-periodic signal.
A dilution term, D TESS is included for each of the two TESS transits as the larger pixel scale can lead to contamination from background stars that alter the transit signal.We assume the ground based transits do not ex- perience flux contamination due to their higher spatial resolution and use them to correct the TESS photometry.We fit a separate term for each sector due to variations in the placement of the target and background stars on the detector pixels.A linear trend was seen in the residuals of the g band transit from the LCOGT, which was fit out prior to inclusion in the joint fit.
The final derived planet parameters from the joint fit are included in Table 3 and the phased transits and RVs are shown in Fig 1 .While we estimate a non-zero eccentricity, we note that this is dominated by the PFS observations as evinced by the low RV jitter compared to NEID.Previous work on HATS-6 b also derived a non-zero eccentricity using PFS which was inconsistent with stellar parameters and it was suggested there may be an underestimation of the formal RV errors when it comes to faint targets (Hartman et al. 2015).While we formally adopt the eccentric fit here in the interest of completeness, we caution against over-interpretation of this tentative eccentricity detection.−0.06 M J , and density 2.0±0.2g cm −3 .Orbiting a metal-rich M-dwarf host star (M = 0.63±0.02M ), it joins a small but growing number of GEMS with precise masses and radii.Additionally, at a mass ratio of ∼ 0.39%, it is just below the new limit of 0.4% for exoplanets that has been proposed by the IAU (Lecavelier des Etangs & Lissauer 2022).In Figure 2, we plot TOI-4201 b together with all known transiting giant planets (R > 4 R ⊕ ) around M-dwarfs, up to T eff <4100 K to account for host stars on the late K/early M border in the first panel.Subsequent panels use a cutoff of T eff <4000 K.
As shown in Figure 2, the closest grouping in massradius space to TOI-4201 b consists of HATS-76 b and HATS-77 b (Jordán et al. 2022), both of which orbit stars that are on the edge between late K and early M-dwarfs.Of particular note is the similarity in radius between TOI-4201 b and HATS-77 b; despite both planets orbiting old inactive stars5 , and neither planet experiencing sufficient irradiation to inflate the radius (Demory & Seager 2011), both have somewhat larger radii than models would suggest (Mordasini et al. 2012).In order to characterize the interior of TOI-4201 b and quantify the degree of inflation, we use the giant planet evolution models from Müller & Helled (2021) and calculate the cooling for different possible heavy-element masses.Using the derived planetary parameters, these models suggest that the planet is inflated by ∼ 5-10% beyond what would be expected for a planet with a pure H/He composition as shown in Figure 3.
Previous work has found that the assumptions underlying interior models of giant planets can cause variations in the derived radius on the order of a few percent (for a review, see Müller & Helled 2023).The main culprits appear to be uncertainties in the H-He equations of state (Müller et al. 2020a;Howard & Guillot 2023) and the opacity (Müller et al. 2020a).Additionally, current evolution models of giant exoplanets assume that their interiors are fully convective.However, there is evidence that Jupiter and Saturn have regions that are not fully convective today due to composition gradients (Debras et al. 2021;Mankovich & Fuller 2021).This could suggest that giant planets do not cool predominantly by large-scale convection, and therefore their interiors may be hotter, leading to inflated radii of up to about 10% past 1 Gyr (Kurokawa & Inutsuka 2015).However, it is currently unclear whether these composition gradients could be sustained over a few Gyr (Müller et al. 2020b).To explain the apparent inflation of TOI-4201 b, either its cooling must be slowed by the aforementioned mechanisms, or there must be another process heating the interior.However, a detailed investigation of this would require next-generation evolution models, and is beyond the scope of this work.

Planet formation and migration
Core accretion is a model of planet formation by which small solid particles coagulate and gradually grow to planetary embryos either through pebble or planetesimal formation.These embryos can be massive enough to trigger runaway gas accretion and allow for planets to retain large H/He dominated atmospheres (Pollack et al. 1996).While this is the favored model for close in planets, the decreased solid mass available (Andrews et al. 2013;Pascucci et al. 2016) and the increased Keplerian orbital timescales (Laughlin et al. 2004) around M-dwarfs make reaching runaway accretion and forming giant planets a challenge.
We use a basic mass budget framework to determine the possibility of TOI-4201 b forming through coreaccretion.TOI-4201 b seems to be inflated beyond what current theoretical models for giant planets support, so we are unable to use available planetary inte- a The reported values refer to the 16-50-84% percentile of the posteriors.
b In addition to the "Absolute RV" from Table 2.
c Jitter (per observation) added in quadrature to photometric instrument error.f We assume the planet to be a black body with zero albedo and perfect energy redistribution to estimate the equilibrium temperature.rior and evolution models such as planetsynth meaningfully to find a heavy-element content or bulk metallicity (Müller & Helled 2021).Instead, we make a lower-bound estimate of the heavy-element content by assuming the planet will have approximately the same metallicity as its host star.The composition of the Sun is approximately 1.39±0.06%metals by mass (Asplund et al. 2021).Combining this with a metallicity of 0.30±0.15dex for TOI-4201 yields an estimated heavy-element mass of 23±8 M ⊕ for TOI-4201 b.Conservatively, if we assume 10% formation efficiency (Liu et al. 2019) this would require a minimum dust mass of ∼ 230 ± 80 M ⊕ to have been present for planet formation in the Class II disk.In the Lupus association, the dust available in Class II disks around M dwarfs ranges from 1-50 M ⊕ (Manara et al. 2022), suggesting it is unlikely there would be sufficient material for TOI-4201 b to form through core accretion in these disks.
There are caveats to this simple argument; the efficiency of the core accretion mechanism is poorly constrained and our understanding of the dust masses contained in disks when planet formation begins is incom- Figure 3.The radius evolution of a planet with the same mass and equilibrium temperature as TOI-4201 b assuming different core (heavy-element) masses, following the model described in Müller & Helled (2021).The grey shaded region represents the area of this parameter space in which our analysis places TOI-4201 b, a region substantially above the radius expected for a planet composed purely of H/He.As there are no indications this is a particularly young system, we include a dashed line at 1 Gyr to suggest a lower limit on the age of the planet.
plete.In our argument we have assumed a 10% efficiency, but simulations have shown this value is strongly dependent on both the turbulence of the disk and the fragmentation velocity of individual particles.Under reasonable disk conditions this efficiency may range from ∼1% up to 40% (Guillot et al. 2014;Chachan & Lee 2023).Ring structures in the protoplanetary disk may also increase formation efficiency and leave behind wide planetesimal belts with distinct profiles exterior to the planets in a system (Jiang & Ormel 2023).While a debris belt is unlikely to be detected around an older star like TOI-4201, the detection of younger GEMS may present an opportunity to look for this profile as a clue towards the formation environment of similar objects.
The dust mass of a disk is typically derived by assuming blackbody emission and using single wavelength continuum flux measurements, often in the millimeter/submillimeter regime, combined with an assumption that the emission is optically thin (Hildebrand 1983).Recent results have shown this is likely to underestimate the solid mass in the disk; full radiative transfer modeling across multiple wavelengths has shown continuum emission is likely optically thick (Michel et al. 2022;Xin et al. 2023) and dust masses derived by SED fitting are greater than the analytical estimate by a factor of 1.5 -5 (Rilinger et al. 2023).One possibility is that a significant amount of the dust mass is contained in larger bodies to which millimeter observations are not sensitive (Najita & Kenyon 2014).The exact amount contained in such bodies is unknown, but the upcoming Next Generation Very Large Array (ngVLA) will be sensitive to cm-sized grains and its higher resolution will allow us to probe disk substructure caused by low mass planets (Selina et al. 2018(Selina et al. , 2022)).
Due to the gaps in our understanding, we cannot completely rule out core accretion as a formation mechanism for TOI-4201 b.However, if we assume a substantial amount of the dust in a Class II disk is contained in larger bodies, this could be indicative of planet formation beginning earlier in the lifetime of the disk (Tsukamoto et al. 2017).Formation beginning in the Class I/0 phase can also decrease the accretion timescale, forming a core large enough to trigger runaway gas accretion within 0.5 Myr (Tanaka & Tsukamoto 2019).Ongoing work with ALMA is looking for substructures in these early disks to better constrain when planet formation begins (Ohashi et al. 2023).
If planet formation efficiency is <10% then gravitational instability, whereby instabilities in massive circumstellar disks allow for the material to directly collapse into giant planet cores on short timescales (∼10 3 -10 4 years after the disk becomes unstable Boss 1997), could be a viable alternative formation mechanism for TOI-4201 b.While the hope of differentiating between the two mechanisms and therefore gaining insight into formation efficiency through atmospheric spectroscopy is desirable (Hobbs et al. 2022), the complex and uncertain nature of the models makes any individual results more qualitative than quantitative (Mollière et al. 2022).Atmospheric spectroscopy of GEMS is now beginning with JWST (and ARIEL in the future), and may help shed light on lines of attack, or lack thereof, to this problem.

SUMMARY
We present the discovery of TOI-4201 b, a Jovian exoplanet with an inflated radius orbiting an early M dwarf.The planet was first identified from TESS photometry and follow up observations consisting of ground based photometry, RVs, and speckle imaging constrained the orbital parameters and allowed for characterization of the planet.
The TOI-4201 b mass ratio of ∼0.39% is one of the highest known for transiting planets around M dwarfs.Assuming a 10% formation efficiency (Liu et al. 2019) and a stellar/sub-stellar atmospheric metallicity and the corresponding planetary heavy-element content of ∼20 M ⊕ would require a disk with a dust mass of ∼200 M ⊕ .Better estimates of disk dust mass demon-strate the most massive disks may reach this threshold, but most remain below the threshold.The existence of TOI-4201 b is suggestive of planet formation beginning before the Class II disk phase.(Oliphant 2007;Virtanen et al. 2020), SERVAL (Zechmeister et al. 2018), starry (Luger et al. 2019;Agol et al. 2020), Theano (The Theano Development Team et al. 2016), planetsynth (Müller & Helled 2021).
identified as hosting a transiting object of interest in the TESS Sector 6 longcadence (1800 s) light curve spanning 2018 December 11 to 2019 January 7, by the Quick Look Pipeline (QLP; Huang et al. 2020) during the Faint-Star Search (Kunimoto et al. 2022).It was re-observed by TESS during Sector 33 from 2020 December 17 to 2021 January 13 with 600 s candence.We extract a light curve for each sector by performing aperture photometry with eleanor (Feinstein et al. 2019) using the CORR_FLUX values, in which eleanor uses linear regression with pixel position, measured background, and time to remove signals correlated with these parameters.We show the TESS photometry in Figure 1.
Placing TOI-4201 b in context TOI-4201 b is a massive Jovian planet with a radius of 1.17±0.04R J , mass 2.59+0.07

d
Dilution due to presence of background stars in TESS aperture, not accounted for in the eleanor flux.e We use a Solar flux constant = 1360.8W/m 2 , to convert insolation to incident flux.

Figure 1 .
Figure 1.Figure including all photometric and RV observations used in the analysis of TOI-4201 b.(a) Timeseries of RV observations of TOI-4201 b with NEID (red) and PFS (green).The best fitting model from the joint fit is plotted in blue and residuals after subtracting the fit are included below.(b) NEID and PFS observations phase-folded on the best fit orbital period from the joint fit, with the model in blue and 16-84% confidence levels in light blue.(c)-(i) Photometric observations of TOI-4201 b; in all plots, the detrended data are in grey and the model and 16-18% confidence levels are in blue.In each figure, the point at x=-0.075 represents the median uncertainty in the photometric data.

Figure 2 .
Figure 2. Upper left: We include TOI-4201 b (circled in black) on a mass-radius plane, with planets colored by host star temperature.We include planets around FGK stars in the background and grey contours indicate bulk densities of 0.3, 1.0, and 3.0 g/cm 3 (left to right).Upper right: TOI-4201 b in an insolation-radius plane for the same sample of planets.Lower: Planet to star mass ratio vs. host star mass, colored by equilibrium temperature.Planets with a true mass measurement from transit observations are represented by circular points and planets with a true mass measurement from astrometry are squares, while triangles are minimum masses (RV only).Around M-dwarfs, TOI-4201 b has the highest mass ratio for transiting planets and the planet with the highest mass ratio overall is GJ 463 b (Endl et al. 2022; Sozzetti 2023).

Table 2 .
Summary of stellar parameters for TOI-4201.
. ACKNOWLEDGEMENTS Data presented herein were obtained at the WIYN Observatory from telescope time allocated to NN-EXPLORE through the scientific partnership of the National Aeronautics and Space Administration, the National Science Foundation, and NOIRLab.This work was supported by a NASA WIYN PI Data Award, administered by the NASA Exoplanet Science Institute.These results are based on observations obtained with NEID on the WIYN 3.5 m telescope at KPNO, NSF's NOIRLab under proposal 2022B-785506 (PI: S. Kanodia), managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the NSF.This work was performed for the Jet Propulsion Laboratory, California Institute of Technology, sponsored by the United States Government under the Prime Contract 80NM0018D0004 between Caltech and NASA.WIYN is a joint facility of the University of Wisconsin-Madison, Indiana University, NSF's NOIRLab, the Pennsylvania State University, Purdue University, University of California-Irvine, and the University of Missouri.The authors are honored to be permitted to conduct astronomical research on Iolkam Du'ag (Kitt Peak), a mountain with particular significance to the Tohono O'odham.Data presented herein were obtained at the WIYN Observatory from telescope time allocated to NN-EXPLORE through the scientific partnership of NASA, the NSF, and NOIRLab.Some of the observations in this paper made use of the NN-EXPLORE Exoplanet and Stellar Speckle Imager (NESSI).NESSI was funded by the NASA Exoplanet Exploration Program and the NASA Ames Research Center.NESSI was built at the Ames Research Center by Steve B. Howell, Nic Scott, Elliott P. Horch, and Emmett Quigley.This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia),processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium).Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.Some of the observations in this paper were obtained with the Samuel Oschin Telescope 48-inch and the 60inch Telescope at the Palomar Observatory as part of the ZTF project.ZTF is supported by the NSF under Grant No. AST-2034437 and a collaboration including Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the Uni-versity of Maryland, Deutsches Elektronen-Synchrotron and Humboldt University, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, Trinity College Dublin, Lawrence Livermore National Laboratories, and IN2P3, France.Operations are conducted by COO, IPAC, and UW.This work makes use of observations (Proposal ID: KEY2020B-005) from the Sinistro imaging camera on the 1m Dome B telescope at Cerro Tololo Inter-American Observatory, operated by the Las Cumbres Observatory global telescope network (LCOGT).Computations for this research were performed on the Pennsylvania State University's Institute for Computational and Data Sciences Advanced CyberInfrastructure (ICDS-ACI).This content is solely the responsibility of the authors and does not necessarily represent the views of the Institute for Computational and Data Sciences.The Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University, the Eberly College of Science, and the Pennsylvania Space Grant Consortium.Some of the data presented in this paper were obtained from MAST at STScI.Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts.This work includes data collected by the TESS mission, which are publicly available from MAST.Funding for the TESS mission is provided by the NASA Science Mission directorate.This research made use of the (i) NASA Exoplanet Archive, which is operated by Caltech, under contract with NASA under the Exoplanet Exploration Program, (ii) SIMBAD database, operated at CDS, Strasbourg, France, (iii) NASA's Astrophysics Data System Bibliographic Services, and (iv) data from 2MASS, a joint project of the University of Massachusetts and IPAC at Caltech, funded by NASA and the NSF.This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France, and NASA's Astrophysics Data System Bibliographic Services.This research has made use of the Exoplanet Followup Observation Program website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program CIC acknowledges support by NASA Headquarters through an appointment to the NASA Postdoctoral Program at the Goddard Space Flight Center, administered by USRA through a contract with NASA.GS acknowledges support provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51519.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.