The SPHINX M-dwarf Spectral Grid. I. Benchmarking New Model Atmospheres to Derive Fundamental M-dwarf Properties

Our 1D radiative–convective equilibrium model builds upon that from the ScCHIMERA tool primarily developed for extrasolar planet atmospheres described in Arcangeli et al. (2018), Bonnefoy et al. (2018), Kreidberg et al. (2018), Piskorz et al. (2018), Gharib-Nezhad & Line (2019), Gharib-Nezhad et al. (2021a), and Mansfield et al. (2021). The model iteratively solves for vertical profiles of volume mixing ratio, cloud/condensate properties, thermal structure, and the top-of-atmosphere disk-integrated stellar spectrum. We then generate a grid of model spectra and atmospheres, ⁷ spanning effective temperatures (T_eff) from 2000 to 4000 K in steps of 100 K, surface gravity (log g) from 4.0 to 5.25 in steps of 0.25, metallicity ([M/H]) ranging from −1 to 1 in steps of 0.25, and carbon-to-oxygen ratio (C/O) from 0.3 to 0.9 in steps of 0.2, all with uniform spacing. The components of our model pipeline are outlined below.

• 1.  
  Radiative Transfer. We follow the two-stream source function approximation described in Toon et al. (1989) and Petty (2006). We compute radiative transport through the atmosphere constrained by hydrostatic equilibrium in a plane-parallel geometry, solving for the temperature–opacity combination that results in net zero flux divergence across each atmospheric layer (Marley & Robinson 2014; Hubeny 2017) using the Newton–Raphson iterative scheme described in McKay et al. (1989).
• 2.  
  Equation of State. We use the NASA–CEA package (Gordon & Mcbride 1994) to compute equilibrium chemical abundances of molecules as well as neutral and singly ionized species. This package can compute equilibrium chemical abundances for thousands of species; however, we only use it to compute abundances for a subset of key molecules and elemental species characteristic of M-dwarf atmospheres (notably metal oxides, hydrides, and neutral and some doubly ionized species).
• 3.  
  Opacity Sources. We take into account opacities relevant to cool stellar atmospheres, including CH₄, CO, CO₂, NH₃, PH₃, HCN, C₂H₂, H₂S, H₂–H₂/He collision-induced absorption (Freedman et al. 2014; Lupu et al. 2022), H₂O (Polyansky et al. 2018), metal hydrides SiH (Yurchenko et al. 2017), TiH, FeH, CaH, CrH, AlH, MgH, and metal oxides TiO, VO, SiO (EXOPLINES, Gharib-Nezhad et al. 2021a, 2021b and references therein). We also include neutral atomic lines Na, K (Allard et al. 2016, 2019), Fe i, Mg i, Ca i, C i, Si i, Ti i, O i, ionized atoms Fe ii, Mg ii, Ti ii, C ii, Ca ii (Kurucz 1993, 2011; Piskunov et al. 1995), and bound–bound and bound–free H-continuum (John 1988). Our molecular opacities are continually updated based on the latest line lists from ExoMol (Tennyson et al. 2016, 2020). Opacity functions can be dense due to numerous transitions as a function of temperature leading to millions of lines within any spectral window, depending on the resultant resolution. Therefore, to overcome the computational burden caused by integrating over the opacity function in a tedious line-by-line manner, we use the correlated-k resort–rebin approach as described in Amundsen et al. (2017). Here, the precomputed cross sections are effectively sampled with line-by-line bins mapped onto a smooth cumulative distribution function of line strengths per bin while keeping the ranking of line strengths consistent through the atmospheric column (Marley et al. 1996; Burrows et al. 1997; Fortney et al. 2010; Grimm & Heng 2015). The line transmittance integral is computed by performing a Gaussian quadrature sum using the cumulative distribution transformation, choosing tens of quadrature points per spectral bin (we use about 20 points). This opacity transformation allows for faster computation of the transmission function for opacities without significant information loss (Marley et al. 1996; Burrows et al. 1997; Showman et al. 2009; Fortney et al. 2010; Line et al. 2013, 2014, 2015). In Figure 2, we show molecular, atomic, and continuum opacities included in our model suite at 3200 K temperature and 1 bar pressure. In fact, the choice of line list used for computing key M-dwarf molecular opacities can show a difference in spectral shape by over a factor of 2 with residual differences in spectral shape of up to 125% and a difference in the thermal profile of up to 60 K (Gharib-Nezhad et al. 2021a). Therefore, it is absolutely crucial to use the most up-to-date line lists with appropriate pressure-broadening treatment (see Gharib-Nezhad et al. 2021a), particularly for key molecules that dominate M-dwarf atmospheres such as metal oxides (TiO, VO, H₂O) and metal hydrides (MgH, CrH, FeH, CaH, SiH, TiH), which is the highlight of our synthesis model.
• 4.  
  Convection. Our convection model follows mixing length theory (Böhm-Vitense 1958) to compute heat diffusion through a turbulent medium. We solve for turbulent diffusivity and upward transport velocity using a nonzero convective heat flux in the adiabatic region of the atmosphere. For simplicity, we assume a mixing length parameter value of 1. We use the Schwarzschild criterion to assess whether conditions for convection are met and therefore solve for the convective flux through each stratum and self-consistently update the thermal profile for regions with steep gradients. Mixing length theory also allows for a continuous transition from a fully radiative to a fully convective solution, providing better control over P–T profile convergence and permitting detached convection zones (Hubeny 2017).
• 5.  
  Dust/Clouds. Dust/cloud particles form when solid grains or liquid droplets condense due to the gas vapor pressure surpassing the local saturation vapor pressure (Ackerman & Marley 2001; Helling et al. 2008). Contributions from condensate opacity are calculated from the total absorption and scattering efficiencies given the indices of refraction and particle size distributions. Additionally, the particle sizes as a function of altitude are assumed to follow a log-normal distribution (similar to Ackerman & Marley 2001; Sharp & Burrows 2007). For the current grid-model spectra, we assume the M-dwarf atmospheres to be cloud-free; however, this model capability is readily available to expand to late and/or low-surface-gravity M dwarfs, with atmospheres including (but not limited to) Mg₂SiO₄, MgSiO₃, CaTiO₃, Na₂S, KCl, and Fe.

Zoom In Zoom Out Reset image size

In Figure 3, we show slices of our model grid versus PHOENIX-ACES ⁸ (Husser et al. 2013) (smoothed to R ∼ 250) varying over T_eff, log g, and [M/H], which we define as the sum of all elemental abundances relative to hydrogen and helium, with the solar value being 0.0 in log space. Overall, we demonstrate that our models are consistent when compared to the widely used stellar model grids that are available to the community today, with noticeable differences, particularly for the M dwarfs with greater than solar metallicity and early M dwarfs with T_eff > 3500 K, a likely consequence of using the most recent line lists.

Zoom In Zoom Out Reset image size

We analyze M-dwarf data in two groups. The first is widely separated M+G binaries (WBS) with empirically calibrated metallicities ([M/H]) from the host stars (Mann et al. 2013). The spectra obtained were observed with the SpeX spectrograph (Rayner et al. 2003) and SNIFS (Lantz et al. 2004) and cover a combined visible + near-IR (NIR) range spanning 0.32–2.4 μm (Mann et al. 2013). Both sets of observations were carried out by Mann et al. (2013) with the SpeX using the NASA Infrared Telescope Facility (IRTF) on Maunakea with the SXD mode (03 × 15'' slit at R ∼ 2000) and the SNIFS observations with the University of Hawaii 2.2 m telescope (with R ∼ 1300). Mann et al. (2013) then constructed the combined spectrum we use for our first group of M-dwarf benchmark targets by stitching both sets of data within the overlapping region of 0.81–0.96 μm and performing appropriate corrections for radial velocity (RV) such that the offsets between the two sets are comparable to random errors in their RV measurements (Mann et al. 2013). Our second group focuses on nearby M dwarfs from an interferometric sample (IS) with measured angular diameters from Boyajian et al. (2012). For this group, we obtain data observed by Mann et al. (2015) using SNIFS (Aldering et al. 2002; Lantz et al. 2004) and STIS (Bohlin et al. 2001; Bohlin & Gilliland 2004) on the HST that cover a total optical bandpass of 0.32–0.97 μm with R ∼ 1000 and NIR Spex (Rayner et al. 2003) covering 0.8–2.4 μm at R ∼ 2000 (Mann et al. 2015). The SNIFS, SpeX, and STIS data were again stitched, H-band flux-calibrated, and renormalized to minimize the relative offsets between the overlapping regions (see Table 1, and Mann et al. 2015 for further details on the observations and data reduction procedures).

For both groups of M dwarfs, we focus purely on the optical/NIR wavelength range as it is ideal for inferring molecular abundances and deriving bulk chemical compositions seamlessly (Line et al. 2015). Furthermore, we restrict ourselves only to the low-resolution regime, which provides a way to leverage broadband molecular features over a wide wavelength range without incorporating non-LTE effects or microturbulence—a concern when working with higher-resolution spectra. Additionally, Jofré et al. (2014) showed that varying instrument resolution resulted in metallicity differences of less than 0.05 dex for benchmark FGK stars they analyzed and less than 0.1 dex for stars with T_eff below 5000 K. Such a paradigm shift to lower resolutions has also been demonstrated in Galactic abundance synthesis models and brown dwarfs (Line et al. 2015, 2017; Ting et al. 2018; Zalesky et al. 2019, 2022; Zhang et al. 2021a, 2021b). Our model grid spectra are computed at R ∼ 250 (see Figure 3); however, we smooth both the data and model to R ∼ 120 for this analysis to match a generic SpeX Prism Library spectrum (Burgasser 2014) as previously demonstrated for benchmark T dwarfs by Zhang et al. (2021a) and mainly to demonstrate the strength in our methodology in the regime of low spectral resolution.

We use the Bayesian inference tool Starfish ⁹ (Czekala et al. 2015b, 2018) as our primary grid-model interpolation method, as successfully used for M/G-type stars (Czekala et al. 2015b, 2018; Gully-Santiago et al. 2017), benchmark brown dwarfs (Zhang et al. 2021a, 2021b), and protostars (Greene et al. 2018). We construct a spectral emulator with Starfish, specifying the range of our grid of models: 2000 K ≤ T_eff ≤ 4000 K, −1 ≤ [M/H] ≤ 1, 0.3 ≤ C/O ≤ 0.7, and 4.0 ≤ log g ≤ 5.5. Similar to Zhang et al. (2021a), we reduce the spectral resolution of our model spectra for all grid points using a Gaussian kernel matching the resolving power delivered by SpeX in Prism mode with a 03 slit (Rayner et al. 2003) to be R ∼ 120. The primary function of the emulator is to interpolate by first deconstructing the convolved model grid into eigenspectra using principal component analysis (PCA). In this step, the models are evaluated with a training set and a testing set to generate eigenspectra of PCA components, the mean spectrum over the entire grid, the standard deviation spectrum, as well as associated weights capturing model grid systematics (Czekala et al. 2015b, 2018). Then, for any given set of parameters spanning our grid, a reconstructed (or interpolated) spectrum is computed as a linear combination of these. Emulator reconstruction systematics do not tend to dominate the forward modeling error budget, as demonstrated by Zhang et al. (2021a) in the case of Sonora–Bobcat models (Marley et al. 2021). A Gaussian process kernel is used to solve for an emulator covariance matrix comprising the standardized weights for each PCA eigenspectrum, providing a probability distribution of the weights, a distribution of the interpolated spectra, as well as the propagated interpolation uncertainties (see Czekala et al. 2015b, 2018 for further details).

The key feature of Starfish is the construction of a composite covariance matrix that also accounts for the total error budget with both parameter estimation and model evaluation. Components of the full covariance matrix are inferred via retrievals (referred to as hyperparameters) to allow for realistic error analysis driven by the data. In addition to the emulator covariance matrix described above, Starfish (Czekala et al. 2015a) also computes a noise covariance matrix, which accounts for differences in the form of data–model uncertainties and observed flux uncertainties within and between wavelength pixels, inferred as parameter a_N . Next, a Gaussian process global covariance matrix is computed, which uses a squared exponential kernel accounting for correlated residuals caused by oversampling the instrument line-spread function or deficiencies in the atmosphere model as well as systematics such as insufficient opacity data are all inferred as parameters log(a_G ) and log(l) (Czekala et al. 2015b, 2018; Zhang et al. 2021a). With this method, the traditional evaluation of the likelihood function within an inference problem is modified as a linear combination of these covariance matrices, ultimately leaving us to solve for the following eight parameters in our grid-model retrievals: T_eff, log g, [M/H], C/O, R_*, as physical parameters, and log(a_G ), log(l), and a_N as hyperparameters (see Figure 4). Note that we essentially fit for the (R_*/d)² scaling parameter as a proxy for the radius R_* constrained by distances measured from the Gaia Collaboration et al. (2016) and Gaia Collaboration (2018).

Our grid-model fitting pipeline allows for priors on both physical parameters and the hyperparameters. For the effective temperature, surface gravity, metallicity, carbon-to-oxygen ratio, and radius, we use uniform priors spanning the range of the grid (2000 K ≤ T_eff ≤ 4000 K, 4.0 ≤ log g ≤ 5.0, −1 ≤ [M/H] ≤ 1, 0.3 ≤ C/O ≤ 0.7 and 0.01 R_⊙ ≤ R_* ≤ 0.8 R_⊙). For efficiency, we truncate the grid ranges for each parameter axis informed by our linear interpolation routine (as done in Zalesky et al. 2019), which precedes the Starfish analysis, for each target in our sample. We found that this speeds up the matrix reconstruction by the Starfish emulator for a PCA training set over the grid points by about an order of magnitude. We define physical priors as uniform within the parameter ranges of our subgrid, specific to each target.

For hyperparameters, we follow recommendations from Zhang et al. (2021a) and Czekala et al. (2015b) and use a Gaussian prior on a_N with a mean of 1.0 and width of 0.25 over all positive values to sufficiently describe the "noise covariance matrix" that accounts for covariance in data–model residuals within each pixel. We use uniform priors on both "global covariance matrix" parameters l and log(a_G ) to be between [200, 2000 × 5] km s⁻¹ and (−∞, ∞) respectively for our stitched SNIFS+SpeX (Rayner et al. 2003; Lantz et al. 2004) spectrum (Mann et al. 2013, 2015), which describes the behavior of the covariance values between two wavelengths, the SXD-mode line-spread function, and the amplitude of the Gaussian process kernel to quantify correlated residuals from model systematics. Although we did not rigorously compute the optimum autocorrelation wavelength to infer the appropriate range for the wavelength scale parameter (l) given our choice of data as done in Zhang et al. (2021a), we did assume the average offset in the stitching of the two spectra (SNIFS+SpeX, Rayner et al. 2003; Lantz et al. 2004) to range between 30 and 50 Å as indicated by Mann et al. (2013, 2015) (which corresponds to around 429–829 km s⁻¹ using Equation (A3) of Zhang et al. 2021a). We simply choose uniform priors on l to comfortably accommodate this range as well as those prescribed by Zhang et al. (2021a) when working with SpeX 08 and 05 slits. In a future analysis, fine-tuning the inferences on the wavelength scale parameter (l) as well as the amplitude of the Gaussian process kernel log(a_G ) will help in quantifying the impact of other possible sources of systematic noise that could influence the residuals.

With a grid-model fitting process, it is straightforward to consider a likelihood function evaluated between the model spectrum and data assuming the residuals are pixel-independent, meaning that the covariance matrix describing the residuals only has diagonal terms, reflecting the measurement uncertainty in the pixel. The remainder of the uncertainty is captured with an additional free parameter simply inflating the error bars. The true uncertainty budget, however, can be properly characterized by using a nontrivial covariance matrix to capture correlated residuals within a likelihood function to be maximized, thus accounting for pixel-to-pixel variation and quantifying the probability of the data given the model, where large residuals are penalized. In essence, following Equation (7) of Czekala et al. (2015a), the log-likelihood function evaluated within Starfish is as follows:

$$\begin{eqnarray}&&\mathrm{ln}\ p(D| M)=-\displaystyle \frac{1}{2}({R}^{T}{C}^{-1}R+\mathrm{ln}\ \det \ C+{N}_{\mathrm{pix}}\mathrm{ln}\ 2\ \pi ).\end{eqnarray} \tag{ 1 }$$

Here, R is the residual spectrum, C is the N_pix × N_pix covariance matrix, which is a sum of two covariance matrices, and N_pix is the noise in each element of the residual spectrum. This first nontrivial covariance matrix comprises the measurement uncertainty in each pixel/element of the residual spectrum and the covariance between each pair of these pixels. In practice, Starfish constructs the covariance matrix by parameterizing with a Gaussian process (GP) kernel that describes the individual covariances between pairs of pixels or wavelength bins (via the noise and global covariance matrices with their associated priors described in Section 2.4). Additionally, model deficiencies from assumed physics are taken into account in Equation (1) via mapping differences spanning the model grid (through PCA decomposition) when constructing the emulator object for interpolation (described in Section 2.3). In other words, every model call that is a reconstructed/interpolated spectrum is defined in the following way, as per Equation (23) of Czekala et al. (2015a):

$$\begin{eqnarray}&&{f}_{\lambda }({{\theta }_{* }}_{i}^{\mathrm{grid}})\approx {\xi }_{\mu }+{\xi }_{\sigma }\displaystyle \sum _{k=1}^{m}\ {w}_{k}({{\theta }_{* }}_{i}^{\mathrm{grid}}){\xi }_{k}\end{eqnarray} \tag{ 2 }$$

where f_λ is an interpolated model, ξ_μ is the mean spectrum of the entire synthetic model grid, ξ_σ is the standard deviation spectrum, and w_k is the weight associated with the kth eigenspectrum. The interpolation/reconstruction uncertainty is thus enveloped in a second nontrivial covariance matrix, which is parameterized yet again with GP kernels that describe the probability distribution of weights associated with the PCA eigenspectra. These weights are modeled with joint multivariate Gaussian distribution kernels spanning each dimension of the synthetic grid (T_eff, log g, [M/H], and C/O in our case), thereby capturing reconstruction/interpolation residuals for each wavelength bin as well as between the bins, in addition to downweighting grid deficiencies. Finally, the likelihood function is evaluated over the defined priors for both the physical parameters and hyperparameters describing the properties of the GP kernel.

First, we assess the influence of photospheric heterogeneity in the retrieved abundances. To answer this question, we include a parameterization (Rackham et al. 2017, 2018; Apai et al. 2018; Pinhas et al. 2018; Zhang et al. 2018; Iyer & Line 2020) for the stellar spectrum to be a linear combination of two spectra—one representative of the calm photosphere and one representative of a cooler spot region on the photosphere, which are both weighted by the fractional coverage in area of the spots:

$$\begin{eqnarray}&&{F}_{* }={F}_{\mathrm{phot}}(1-{f}_{\mathrm{spot}})+{F}_{\mathrm{spot}}{f}_{\mathrm{spot}}\end{eqnarray} \tag{ 3 }$$

where F_* is the final stellar model spectrum, F_phot is the spectrum of the stellar photosphere region alone, F_spot is the spectrum of the stellar spot region alone, and f_spot is the fractional coverage in area of the spots. Here, both F_phot and F_spot are model spectra as functions of their respective T_eff, log g, [M/H], C/O, and scaled (R_*/d)² values evaluated (generated with the Starfish interpolation emulator, see Section 2.3) at each instantiation within our retrieval. All parameters are sampled as a single instance from their respective posteriors for each likelihood evaluation within our routine, except for the effective temperature T_eff alone, which is drawn twice: as T_phot for the temperature of the photosphere region and as T_spot for the temperature of the spot region, to ultimately produce the final model spectrum F_*. Following this, we end up with three additional parameters for our grid-model retrievals, namely, T_phot, T_spot, and f_spot. Both the temperature values are bound by uniform priors spanning the range of our model grid (but with the condition that T_spot is always lower than T_phot), and f_spot is taken as a uniform prior between 0 and 1 because we do not have abundant a priori knowledge of physically reasonable spot covering fractions on M-dwarf photospheres (Iyer & Line 2020).

We perform retrievals including this "mixed model" parameterization for two targets with the largest σ deviation from observations in Figure 5, NLTT 10349 (M1.0) and NLTT 11500 (M1.8). In Figure 11, NLTT 10349 shows an improvement from 3σ to 2σ in our ability to constrain [M/H] relative to the empirically derived value. At the same time, NLTT 11500 does not show any further improvements in [M/H] upon the inclusion of stellar spots. However, the retrieval does appear to be artificially pushing the spot covering fraction to be significantly high. We also see that the log g posterior is driven up to higher values, necessitating further investigations while potentially exploiting the mass/metallicity correlation for a mid-M dwarf (see Figure 11).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also test the influence of including stellar spot parameterization on the retrieved radii. For the case of Gl 725B (M3.5)—a target in our sample with the most significant deviation from the interferometric measurements (IS)—we find that the inclusion of stellar heterogeneity ever so slightly improves accuracy in the retrieved radius. The retrieved effective temperature, however, is significantly more overestimated than before, with a significant deviation from interferometry. This could be due to a simplified stellar spot parameterization causing degeneracies with other fundamental stellar properties, which do not always prove to be accurate for active M dwarfs (Iyer & Line 2020) (extra figures provided in our zenodo link). It is also well established by now that there is very little knowledge of physically reasonable coverage fractions of stellar photospheric heterogeneities or of their distribution throughout the photosphere of an M dwarf, which adds to numerous problems that exist in characterizing the active environment of such stars (Rackham et al. 2022). Spots occur in stars primarily due to pockets of magnetic flux tubules stretching outwards beyond the photosphere, suppressing the flux in that region as a result of convection processes that produce an overall cooler effect in this region. Therefore, efforts to model this more accurately, including complementary efforts from the observational side, could provide a great starting point for fine-tuning priors on such correlated parameters and improving parameterizations within a retrieval framework.

Here we explore the influence of the convective mixing length parameter (Böhm-Vitense 1958) on values retrieved from our stellar atmosphere model. We define the mixing length parameter as

$$\begin{eqnarray}&&\alpha \,=\,\displaystyle \frac{{l}_{\mathrm{mix}}}{{H}_{p}}\end{eqnarray} \tag{ 4 }$$

where l_mix is the convective mixing length and H_p is the pressure scale height; the ratio α is a free parameter that is often based on solar-calibrated values (Chabrier & Baraffe 1997). All our models in the grid in Section 2 are calculated assuming α = 1 for simplicity. In Figure 8, we find that most of our retrieved radii are consistent with the interferometry values, except that for Gl 725B (M3.5) is underestimated. In the same figure (right panel), we also find that their effective temperatures are overestimated relative to the interferometry values. To address this, we perform a test by changing the mixing length value to α = 0.5 and generate a mini-grid of models to assess whether it could improve our predictions and explain the systematic deviations seen in both stellar radius and effective temperature for this target. We choose this value as recommended by Mann et al. (2015), who found that reducing the mixing length parameter yielded a reduction in their T_eff values of about 2% and an increase in R_* of 2.5% for low-mass stars <0.3 M_⊙. This is also consistent with findings of Barklem et al. (2002) using the MARCS model (Gustafsson et al. 2008), where simultaneous fits to Hα and Hβ in spectra of metal-poor dwarf stars were sensitive to the mixing length value assumed in the model (α = 0.5 instead of 1.5), which varied the effective temperature by up to 40 K.

We also repeat this exercise with a value of α = 2 following Baraffe et al. (2015), who primarily favor mixing length values above 1.6 for cooler and denser stars with larger convective envelopes. [M/H] could potentially be sensitive to the choice of mixing length based on relevant model calibrations; however, Song et al. (2020) find that the typical scatter in [Fe/H] values due to α is well within the uncertainty levels of metallicity itself.

We also pick Gl 699 for this exercise (Barnard's star, which is known to have subsolar metallicity, [Fe/H] = −0.36, Rojas-Ayala et al. 2012), where our model appears to derive a supersolar [M/H] (0.17) and overestimated T_eff (3450 K, see Figure 8). In Figure 12, the model grid assuming α = 0.5 fares better than the α = 2 models in terms of the retrieved effective temperature, getting it closer to the interferometry value within a shift of 1σ and 2σ from the median respectively for Gl 725B and Gl 699. We do not see any changes in the retrieved radii for either case; however, the retrieved [M/H] value changes to better match the empirically derived value from Mann et al. (2015) under 1σ for Gl 725B and 2σ for Gl 699 compared to observations. In particular, the median retrieved [M/H] with the α = 0.5 model results in a subsolar metallicity for Barnard's star as expected. We also do not see a significant difference in the retrieved surface gravity parameter, log g, for either case, consistent with findings from Mann et al. (2015) where the mixing length parameter did not appear to perturb derived masses significantly. Hence, for very low-mass targets (<0.3 M_⊙) such as Gl 725B and Gl 699, we find that effective temperature and metallicity are sensitive to the mixing length parameter assumed. Reducing α from the solar-calibrated value (α = 1.88, Mann et al. 2015) to a lower value (α = 0.5) is in essence equivalent to reducing the convective envelope and increasing the radiative interior, at the same time lowering the surface temperature and influencing the molecular opacities in accordance. Consistently, we see a slight reduction in model-derived effective temperatures and metallicities for both Gl 725B and Gl 699 that fall within the fully convective regime (<0.3 M_⊙, Figure 12). For the >0.3 M_⊙ target Gl 436, reducing the mixing length parameter showed a dramatic reduction in effective temperature and metallicity, with no noteworthy change in the radius (see supplementary figures). In line with findings of Mann et al. (2015, see their Table 4), we also find that the relative difference in retrieved effective temperature due to reducing the convective parameter is likely larger for >0.3 M_⊙ M dwarfs. However, analyzing a larger sample and the resulting population-level trends will allow for a more thorough assessment (forthcoming in the second paper of this series).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also investigated the effects of stellar spots for the case of Gl 725B alone while keeping the mixing length parameter constant at 1. We found that the model-derived radius is slightly increased by an additional 1σ, thereby shifting the histogram closer to the observational value (supplementary figures). No significant changes occurred in the other derived parameters. Therefore, we recommend that including mixed models that consider the effects of both convection and stellar photospheric heterogeneity could be one way of understanding the nuances that influence retrieved fundamental properties of M dwarfs, particularly when moving toward later spectral types.

The carbon-to-oxygen ratio is an extremely valuable composition that can inform about the balance of carbon and oxygen species chemistry or refractory elements in atmospheres of cool stars and their possible formation scenarios such as fragmentation (Offner et al. 2010) or core collapse (Gagné et al. 2018). C/O ratios also shed light on the density, formation, and migration history of rocky planets or planetary-mass companions harbored by M dwarfs (Gaidos 2000). Several retrieval techniques using low-resolution NIR spectra have acquired robust C/O constraints for cold brown dwarfs (Line et al. 2015, 2017; Zalesky et al. 2019, 2022; Gonzales et al. 2020); in addition, pEW (Nakajima & Sorahana 2016) and synthetic spectral model techniques (Gizis et al. 2016; Veyette et al. 2016) have been employed for M dwarfs. This work, however, is a first in the simultaneous and systematic determination of [M/H] and C/O constraints of benchmark M dwarfs from low-resolution spectra. In Figure 13, we overlay our grid-model-derived [M/H] and C/O values with FGK stars in the local population (Hinkel et al. 2014) (gray), showing overall consistency. Most of our low-mass star sample falls below the solar C/O value of 0.55. We also find that our average C/O constraints of ∼10% have not been previously achieved even for main-sequence G stars, specifically in low-resolution (R ∼ 120) spectra.

Zoom In Zoom Out Reset image size

Deriving the C/O ratios of FGK stars has been challenging yet feasible due to limited blended lines in the spectra and not accounting for non-LTE corrections (Fortney 2012). Delgado Mena et al. (2021) have computed C/O values for over 1000 FGK stars and found that C/O ratios derived from high-resolution (R > 100,000) spectra are extremely sensitive to the choice of oxygen line or stellar atmosphere models used to compute them, leading to systematic shifts that need to be scaled with a solar reference value, particularly for cool stars (K type). In even cooler M dwarfs, the bulk of C goes into the CO opacity and the O into metal oxides (TiO/VO) with dominant absorption features. Given that this falls typically in the optical/NIR region of the spectra that is more susceptible to model assumptions as seen from the above analysis, we plan in the future to perform forward model Bayesian retrievals (Gonzales et al. 2020; Zalesky et al. 2022) to directly solve for the radiative–convective atmospheres of M dwarfs without having to rely on self-consistent grid-model assumptions.