The M-dwarf Ultraviolet Spectroscopic Sample. I. Determining Stellar Parameters for Field Stars

Understanding the physics of low-mass stars, M dwarfs, is of critical importance to current astronomy research. These objects represent the most numerous stars in our Galaxy (Bochanski et al. 2010), and the majority of our nearest stellar neighbors are M dwarfs (e.g., Henry et al. 2018). They are frequent targets of exoplanet searches, revealing an abundant population of rocky planets (e.g., Dressing & Charbonneau 2015), with ongoing searches often dedicated to M-dwarf hosts (e.g., Ricker et al. 2015; Reiners et al. 2018). The fundamental stellar properties, mass, radius, and luminosity, are essential parameters in assessing stellar structure, evolution, and their interplay with other observable phenomena, like magnetic activity. Additionally, these are the stellar parameters necessary to accurately characterize the properties and atmospheres of exoplanetary systems.

Unlike warmer stars, M dwarf properties have been difficult to model, revealing large discrepancies between different models and observed properties; for example, field M dwarfs are known to display larger radii than predicted by evolutionary models (e.g., López-Morales 2007; Kraus et al. 2011), including across the fully convective boundary (Jackson et al. 2018; Kesseli et al. 2018). In the cool M-dwarf atmospheres displaying a large number of molecular lines, model spectra can sometimes fail to reproduce the observed features. This is especially the case at blue optical wavelengths, a discrepancy often attributed to missing line opacities (e.g., Allard et al. 2000, 2012). Over the years, this evidence has led to a focus on utilizing empirical methods to estimate M-dwarf properties, including mass–luminosity relations (e.g., Delfosse et al. 2000; Benedict et al. 2016; Mann et al. 2019), mass–radius relations (e.g., Boyajian et al. 2012; Rabus et al. 2019; Schweitzer et al. 2019), and ever more precise photometric bolometric corrections (e.g., Tinney et al. 1993; Mann et al. 2015). These often rely on eclipsing binaries, or interferometrically observed stars with precise bolometric fluxes, and benefit from improved distance and stellar metallicity estimates. The statistical scatter in these relationships can be calculated and is the basis for accurate values with reliable uncertainties in the estimates of M-dwarf fundamental properties.

Because of the varied development of such relations, the literature estimates of M-dwarf properties contain an assortment of underlying assumptions, in aggregate producing populations that incorporate the uncertainties inherent to the different methods. This approach makes it difficult to assess systematic effects from the property determinations when analyzing populations of stars and exoplanets. To reveal potential biases introduced by the chosen set of stellar parameters, it is important to understand the underlying assumptions using consistently determined stellar fundamental properties, and to accurately characterize their often correlated statistical uncertainties. Analysis of samples with consistently determined properties can thus point to the underlying causes of observed trends or which of the known sources of error need to be improved to reveal new findings. The statistically correlated uncertainties between properties is a feature often ignored in the literature. Regression analyses relying on such estimates are likely to yield biased uncertainties that are enlarged relative to the underlying data, potentially obscuring important astrophysical relationships.

In this paper, we demonstrate how to self-consistently apply multiple empirical relations to simultaneously determine the fundamental stellar properties of field M dwarfs, with reliable uncertainties. Our Monte Carlo methods readily enable the use of our results in additional analyses to account for the statistical correlations in the parameter uncertainties. These methods are flexible to observational inputs, allowing multiple and distinct measurements to inform the stellar properties. Our methods also avoid the approximations inherent to standard error propagation through nonlinear empirical relationships, which can fail to reproduce the distribution tails—the basis for statistically significant results. As part of this effort, we define a new field M-dwarf mass–radius relation, using the literature results of eclipsing binaries.

To illustrate our methods, we introduce the M-dwarf Ultraviolet Spectroscopic Sample (MUSS), a compilation of low-mass stars generally defined as being previously targeted for ultraviolet spectra by the Hubble Space Telescope, but including a few other M dwarfs of interest and M dwarfs with interferometric radii. We treat the field stars of this sample in this work, and address how to assess the properties of young stars in our companion paper (J. S. Pineda et al. 2021, in preparation). Applying our methodology to this sample, we create a set of M dwarfs with consistently determined properties that are already of interest to the astronomy community as planet hosts or as interesting/useful objects in probing stellar magnetic activity. We discuss this sample in Section 2. We define the new field mass–radius relation in Section 3. We then demonstrate and use our unified Bayesian framework for field M dwarfs in Section 4, with applications for the MUSS ensemble stellar properties. We discuss our conclusions in Section 5, and summarize our findings in Section 6.

The Hubble Space Telescope is the only current observatory capable of spectroscopically observing across the entire ultraviolet regime. This unique capability has been used to study low-mass stars in a variety of ways, with two major themes in recent years: (1) the high-energy spectral energy distribution and its impact on planetary atmospheres (e.g., France et al. 2013, 2016, 2020; Bourrier et al. 2017b; Froning et al. 2019), and (2) the nature of these emissions and their temporal variability on long and short timescales as probes of magnetic activity and stellar upper atmospheric structure (e.g., Hawley et al. 2003, 2007; Osten et al. 2006; Loyd et al. 2018a, 2018b, 2021; Kowalski et al. 2019; Pineda et al. 2021).

We chose these HST targets across a variety of programs as the focus of our study of M-dwarf stellar properties. These programs include MUSCLES (HST-12464, 13650, PI-France), Mega-MUSCLES (HST-15071, PI—Froning), FUMES (HST-14640, PI—Pineda), and HAZMAT (HST-14784, PI-Shkolnik), among miscellaneous others cataloged in Table 1. These objects are of interest to the community, and future analyses will benefit from consistently determined stellar properties.

We dub the compilation of these objects the MUSS. The targets are enumerated in Table 1. The sample fundamental properties and their determination are described in Section 4 for field stars, and we defer the discussion of young stars to Pineda et al. (2021). Many of these targets host exoplanetary systems for which we provide consistently determined stellar properties. With uniformly determined stellar properties, this work enables new consistent ensemble statistical examinations of planet property trends—for example, the planet radius valley around low-mass stars (e.g., Cloutier & Menou 2020). Similarly, the methods presented here have already informed the activity analyses of Pineda et al. (2021) and Youngblood et al. (2021). We also show the metallicity distribution of the MUSS field stars used in this paper in Figure 1. Those metallicities are taken from the literature when available on the [Fe/H] scale defined by the work of Mann et al. (2013), or converted to that scale as necessary (Gaidos & Mann 2014; Neves et al. 2014; Mann et al. 2015; Terrien et al. 2015).

Zoom In Zoom Out Reset image size

Although the focus for the MUSS objects is targets with spectroscopic UV data, we included additional stars that may be of interest to the community (exoplanet hosts with close-in or habitable zone planets), as well as a sample of M dwarfs with interferometric angular diameter measurements to provide an additional consistency check for the accuracy of our methods, with certain objects selected from existing catalogs with measured bolometric fluxes to help fill out the entire M-dwarf mass range. Notes on individual stars can be found in Appendix A, with the observational properties of the sample compiled in Table 5 therein. The MUSS stars include 67 targets treated as field objects, and 25 targets with well-determined ages. The latter are treated in our companion paper (J. S. Pineda et al. 2021, in preparation).

Empirical mass–radius relationships have been used considerably in the literature to both assess stellar properties and test stellar models in the low-mass regime (e.g., Chen et al. 2014). We build on previous work to derive a new mass–radius relation for low-mass stars, utilizing an expanded data set and a thorough statistical treatment. This sample is distinct from the MUSS stars, but we will use the results discussed in this section to analyze the MUSS stars in Section 4.

Recently, Schweitzer et al. (2019) published an empirical relationship based on a sample of 55 eclipsing binaries spanning 0.092–0.73 M_⊙. We show their best-fit relation for radius as a function of mass and their reported scatter (0.02 R_⊙, i.e., >3% for ${ \mathcal R }\lt 0.6$ R_⊙) in Figure 2 against interferometric radius measurements and eclipsing spectroscopic binaries at the end of the main sequence from von Boetticher et al. (2019). Their relation is visually consistent with available interferometric measurements of single stars and the eclipsing binaries at the lowest stellar masses, where there are few targets in the Schweitzer et al. (2019) sample (only five stars of ${ \mathcal M }\lt 0.2$ M_⊙). This supports previous findings indicating that there is little difference in the observed masses and radii of single stars versus those in binary systems. Kesseli et al. (2018) stated this conclusion, as both single stars and binary stars show the same extent of radius inflation (∼10%) relative to model predictions. These result demonstrate confidence in the use of such an empirical relation for single stars in the field as well as binaries.

Zoom In Zoom Out Reset image size

We derive a new ${ \mathcal M }$–${ \mathcal R }$ relationship by supplementing the sample used in Schweitzer et al. (2019), the measurements as compiled in their Table B.2, with the 10 additional targets from von Boetticher et al. (2019) to fill in the parameter space for ${ \mathcal M }\lt 0.2$ M_⊙. Although these additional measurements depend on the model-based masses of the FG primaries (von Boetticher et al. 2019), those models do not suffer from large discrepancies relative to observations for Sun-like stars, and because they are consistent with the extant ${ \mathcal M }$–${ \mathcal R }$ relationships, we conclude that the addition of these points improves the reliability of our relation at the lowest masses. We do not include the interferometrically observed single stars, because we will use them as a check on our results and omitting them allows us to confine underlying systematics strictly to assumptions inherent to the eclipsing binary measurements.

The need for a new ${ \mathcal M }$–${ \mathcal R }$ relation is motivated by our desire to assess potential correlations, which were not reported in Schweitzer et al. (2019), in the best-fit parameters, and model the scatter about that fit relationship with the expanded data set. We followed Kelly (2007) in our regression model for stellar radius as a function of mass, accounting for uncertainty in both the mass and radius measurements. That is, each observed data point in mass and radius is jointly sampled from a distribution defined by their observed measurements and uncertainties. The underlying true masses and radii, ${{ \mathcal M }}_{i}$ and ${{ \mathcal R }}_{i}$, are then constrained to follow our model relationship

$$\begin{eqnarray}&&{{ \mathcal R }}_{m}({{ \mathcal M }}_{i})=\sum _{j=0}^{n}{b}_{j}{{ \mathcal M }}_{i}^{j}={b}_{0}+{b}_{1}{{ \mathcal M }}_{i}+...+{b}_{n}{{ \mathcal M }}_{i}^{n},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{{ \mathcal R }}_{i}\,| \,{{ \mathcal M }}_{i},\{b\},{\sigma }_{{ \mathcal R }}\,\sim \,{ \mathcal N }({{ \mathcal R }}_{m}({{ \mathcal M }}_{i}),{\sigma }_{{ \mathcal R }}),\end{eqnarray} \tag{ 2 }$$

where i corresponds to each measured data pair, with predicted radius, ${{ \mathcal R }}_{m}$, at a given mass. The notation of Equation (2) indicates that the random variable ${{ \mathcal R }}_{i}$ follows the distribution ${ \mathcal N }(\mu ,\sigma )$, a normal distribution of mean μ and standard deviation σ, {b} denotes the polynomial coefficients of the empirical relation, and ${\sigma }_{{ \mathcal R }}$ is the intrinsic scatter for the regression model. The polynomial coefficients are treated with uniform priors and ${\sigma }_{{ \mathcal R }}$ uses a log-uniform prior. The likelihood function is thus the product of Equation (2) over the data set, indexed by i, conditioned on the fit parameters {b} and ${\sigma }_{{ \mathcal R }}$. We tested linear and quadratic relations, finding like Schweitzer et al. (2019) that the linear relationship characterized the data set well and that the extra freedom in the higher-order polynomial was not warranted by the data using BIC-like comparisons (see below). We thus focus our discussion on comparing two linear models, one with a constant scatter value, and one with a fractional scatter relative to the predicted radius, i.e., a scatter at a fixed fraction of ${{ \mathcal R }}_{m}$.

For the literature eclipsing M-dwarf binary measurements, we treated the masses and radii as uncorrelated bivariate normal distributions consistent with their measured errors. ³ Incorporating the measurements of von Boetticher et al. (2019), for which many of their reported masses and radii show asymmetric confidence intervals about the modal value, we treated those results within our model as skew normal distributions following Appendix B, to better capture the shapes of the asymmetric posterior distributions. Without the additional 10 data points from von Boetticher et al. (2019), our Bayesian regression model in the constant scatter case gave results consistent within errors to the values determined by Schweitzer et al. (2019).

We sampled the posterior distribution for the fit parameters of our regression model within PyMC3 (Salvatier et al. 2016) using the Hamiltonian NUTS sampler. We tested convergence by examining the trace for each variable, making sure the autocorrelation functions vanish and computing the Gelman–Rubin statistic, which should be very close to unity. We report our best-fit parameters (median/modes) in Table 2 for both the fractional and constant scatter models, with their 68% confidence intervals. We plot our preferred fractional scatter model against the data in Figure 3.

Zoom In Zoom Out Reset image size

Table 2.  ${ \mathcal M }$–${ \mathcal R }$ Relation Model Parameters and Comparison ^a

Model b	b0	b1	${\sigma }_{{ \mathcal R }}$	WAIC	pWAIC	dWAIC	SE
Fractional Scatter	0.0375 ± 0.0038	0.918 ± 0.012	${{ \mathcal R }}_{m}\times$ 0.031 ${\pm }_{0.003}^{0.006}$	−374.43	21.19	0	10.7
Constant Scatter	0.0315 ± 0.0052	0.932 ± 0.013	$0.0143{\pm }_{0.0015}^{0.0023}$	−341.19	19.04	33.23	9.16

Notes.

^a Best model parameters for our ${ \mathcal M }$–${ \mathcal R }$ relation are listed as medians when symmetric and modes when asymmetric, with marginalized 68% confidence intervals. The two coefficients are highly correlated (see Figure 4), and when using our linear relationship they should be sampled from their joint posterior to account for the uncertainty in the relation. ^b We also tested a quadratic model for both fractional and constant scatter models, which yielded WAIC values only marginally better than the linear model, consistent well within their standard errors (SE). We choose the simpler linear models as more representative of the astrophysical ${ \mathcal M }$–${ \mathcal R }$ relation.

Download table as:  ASCIITypeset image

We used the Watanabe–Akaike Information Criteria (WAIC; Watanabe 2010; Gelman et al. 2014) for model comparison. WAIC is useful for assessing how well the models can predict data not included within the regression analysis, accounting for the effective number of free parameters in each model. ⁴ We show this comparison in Table 2, including the WAIC score, pWAIC, which reflects the number of parameters (higher means a more flexible model), dWAIC as the difference relative to the best model, and SE as the standard error on the WAIC assessment. These scores are also consistent with leave-one-out cross validation comparison of the models (see Gelman et al. 2014). Our results show that the fractional scatter model is preferred. In Figure 4, we show the pairwise joint posteriors of this fractional scatter model fit. We use this relation, which has a scatter of 3.1% in radius at fixed stellar mass, for our parameter estimations in Section 4. We compare our new relation to currently available relations in the literature from Boyajian et al. (2012), Schweitzer et al. (2019), and Rabus et al. (2019) in Figure 5. The relations largely agree, and are consistent with the interferometric radius measurements.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In both cases, the best-fit parameters b₀ and b₁ show correlations (see Figure 4). Thus, in practice, when accounting for the uncertainty in our ${ \mathcal M }$–${ \mathcal R }$ relations, these parameters should be drawn from their joint distributions (see Appendix C). We have made these posterior distributions publicly available for both the linear and fractional scatter models, with accompanying code to facilitate their use. ⁵ The fitted sample spans masses of 0.084–0.73 M_⊙, but with few points below 0.09 M_⊙ and above 0.7 M_⊙. We thus recommend 0.09–0.7 M_⊙ as the secure range of applicability for input masses in the field (also see residuals in Figure 3).

While our result agrees with the recent literature, our mass–radius relation has more rigorously characterized scatter, and with the publicly available code, it allows for readers to estimate stellar radii from masses properly accounting for the fit correlations. As a final comparison, we consider how our result compares to the semi-empirical radius relation of Mann et al. (2015). While that work relies on effective temperatures based on atmospheric models, their calibration between absolute K_S -band magnitude and radius is observationally accessible, and uses methods tuned to match empirical radii. It has moreover been used widely for radius estimates in the literature (e.g., Feinstein et al. 2019; Gilbert et al. 2020; Waalkes et al. 2021). To compare, we first converted the absolute magnitudes to masses, and then used our mass–radius relation to plot ${ \mathcal R }({M}_{{K}_{S}})$. That comparison is illustrated in Figure 6 with central 68% confidence intervals accounting for our uncertainties on the ${ \mathcal M }$–${ \mathcal R }$ relation and error propagation through the mass–absolute K_S relation of Mann et al. (2019).

Zoom In Zoom Out Reset image size

The curves match very closely, usually well within the uncertainty. The total errors in determining radii in this fashion are ∼4% (${M}_{{K}_{S}}\to { \mathcal M }\to { \mathcal R }$), slightly larger than the 2.89% quoted by Mann et al. (2015) for just using their semi-empirical method (${M}_{{K}_{S}}\to { \mathcal R }$), although that omits the impact of any uncertainty in their fit parameters, as it is not reported in their article.

These results can be used on their own; however, in concert with other expressions, we reveal a self-consistent and powerful means to fully characterize field M dwarfs in the next section.

4.1. Framework

Estimating stellar properties in the low-mass star regime has heavily relied on empirical relations connecting quantities of interest, like mass and radius, to observables like absolute magnitude (if the distance is known) or potentially spectroscopic features. For example, M-dwarf masses are often estimated from mass–luminosity relationships (e.g., Delfosse et al. 2000) calibrated on binary systems with measured masses. These calibrations are designed to take the observed quantities (e.g., absolute K-band magnitude, ${M}_{{K}_{S}}$) and convert them to the physical attribute of interest (e.g., mass, ${ \mathcal M }$).

However, in Bayesian modeling, the reverse is often desired—we want to be able to predict an observable quantity from the underlying physical properties, for comparison through a likelihood function on the data. Generally, the empirical relations are not invertible, but show a scatter at fixed observable. Even if the relation is invertible, the resulting probability density function deviates from normal. We demonstrate here a method (and use it in Section 4.2) to incorporate these empirical relationships into Bayesian modeling of stellar properties, which can be readily folded into any likelihood function.

To illustrate our methods, we consider a generic empirical mass–radius relationship that gives radius (${ \mathcal R }$) as a polynomial function of mass (${ \mathcal M }$), like Equation (1); however, we rewrite it as

$$\begin{eqnarray}&&f({ \mathcal M },{ \mathcal R })={ \mathcal R }-\sum _{j=0}^{n}{b}_{j}{{ \mathcal M }}^{j}=0.\end{eqnarray} \tag{ 3 }$$

The empirical relationship is then taken to have a measured scatter, ${\sigma }_{{ \mathcal R }}$. This scatter can itself be a function of ${ \mathcal R }$ or ${ \mathcal M }$, as it would be for a fractional scatter. Equation (3) states that empirically this combination of ${ \mathcal M }$ and ${ \mathcal R }$, $f({ \mathcal M },{ \mathcal R })$, is observed to vanish. Thus, f can be taken to follow the probability density function,

$$\begin{eqnarray}&&f({ \mathcal M },{ \mathcal R })\sim { \mathcal N }(\mu =0,\sigma ={\sigma }_{{ \mathcal R }}),\end{eqnarray} \tag{ 4 }$$

where ${ \mathcal N }$ denotes a normal distribution of mean μ and standard deviation σ. We are thus assuming that $f({ \mathcal M },{ \mathcal R })$ vanishes probabilistically within a known scatter, ${\sigma }_{{ \mathcal R }}$. This distribution (Equation (4)) can then be treated as a product within a likelihood function (see below). If uncertainties in the fit coefficients b_j are also known, then these parameters can be drawn from their own respective distributions when sampling ${ \mathcal M }$ and ${ \mathcal R }$. In practice, the b_j are likely highly correlated and should be drawn from their joint posterior. This requires having the output posterior, as we provide for the mass–radius relation we derive in Section 3.

In Figure 7, we show an example of sampling from the single star mass–radius relationship of Boyajian et al. (2012) using Equation (4), with a 5% fractional scatter. The constant fractional scatter of the assumed relation is seen in the narrowing of the sample distribution at lower masses. Combined with an additional observational constraint—for example, the stellar density from a planetary transit light curve—the mass and radius can be constrained individually, accounting for all the observational uncertainty and the scatter of the empirical relation. Although we considered in particular a mass–radius relation, any useful empirical relationship can be recast in the form of Equation (3) and a corresponding probability density function (Equation (4)) written down if an uncertainty about the best empirical relationship is well-defined. This treatment also considers the two quantities, ${ \mathcal M }$ and ${ \mathcal R }$, as separate random variables constrained within the measured empirical relation, not treated deterministically, so that each may vary within the intrinsic scatter of these relationships, e.g., as a result of unaccounted for effects like magnetic activity or metallicity. Recast as in Equation (4), multiple empirical relations can be multiplied within a single likelihood function, assuming they are independent, to simultaneously apply the available empirical constraints; we take this approach when determining the properties of MUSS field stars in the next section.

Zoom In Zoom Out Reset image size

4.2. Determining MUSS Field Properties

We are interested in estimating the stellar properties, mass (${ \mathcal M }$), radius (${ \mathcal R }$), bolometric luminosity (L_bol), and effective temperature (T_eff) of field M dwarfs in a consistent and uniform way. In general, the most widely available data is photometric, and we thus focus our approach on those relations most accessible in their application to any field M dwarf. Additional spectroscopic relations exist in the literature (e.g., Newton et al. 2015; Terrien et al. 2015) that, when independent, can readily supplement the methods used here. As these data are not always available, we proceed with largely photometric empirical relations.

We determined the physical properties within a Bayesian framework through Markov chain Monte Carlo (MCMC) sampling by using a variety of empirical relations for low-mass M dwarfs, outlined below. This allows us to jointly and accurately constrain the stellar ${ \mathcal M }$, ${ \mathcal R }$, and L_bol. We then determine the T_eff from the combination of ${ \mathcal R }$ and L_bol. ⁶ Although T_eff can be estimated spectroscopically, such estimates are often tied to atmospheric models, which potentially introduce new systematics. Our approach for T_eff is grounded in empirical measurements, and directly applies the definition of effective temperature. We used the methods outlined in Section 4.1 to combine multiple empirical relations in order to jointly constrain stellar properties.

Mass–Luminosity: To constrain the mass of our stellar sample, we relied upon the empirical mass–luminosity relationship of Mann et al. (2019), which defines ${ \mathcal M }$ as a polynomial function of absolute K_S -band magnitude, ${M}_{{K}_{S}}$, based on resolved photometry and the measured orbits of a large sample of low-mass binary systems determined from accumulated community efforts. As has been discussed much in the literature, the choice of K-band is due to its relative insensitivity to metal content, more specifically, C and O abundances (Veyette et al. 2016), among stars in the solar neighborhood, providing a tight correlation for ${ \mathcal M }$–${M}_{{K}_{S}}$ (e.g., Delfosse et al. 2000; Benedict et al. 2016; Mann et al. 2019).

Although Mann et al. (2019) also developed a metallicity-based calibration, they do not find the extra variable to significantly improve the fit to their data, and following their recommendations, their fifth-order polynomial is favored for near-solar-metallicity main-sequence stars in the range of 0.08–0.7 M_⊙, yielding ∼2% errors on mass at fixed ${M}_{{K}_{S}}$. Following Section 4.1, we rearranged the results of Mann et al. (2019) to define a normally distributed variable ${f}_{{ \mathcal M }}$ as a function of ${ \mathcal M }$ and ${M}_{{K}_{S}}$, which is empirically observed to vanish with a measured dispersion. We can thus use the probability distribution function for ${f}_{{ \mathcal M }}$ within our Bayesian framework:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}m({M}_{{K}_{S}})=\sum _{i=0}^{5}{a}_{i}\,{\left({M}_{{K}_{S}}-7.5\right)}^{i},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{f}_{{ \mathcal M }}({ \mathcal M },{M}_{{K}_{S}})={ \mathcal M }-m,\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{f}_{{ \mathcal M }}({ \mathcal M },{M}_{{K}_{S}})\sim {{ \mathcal N }}_{{ \mathcal M }}(0,m\,{\sigma }_{m}),\end{eqnarray} \tag{ 7 }$$

where $m({M}_{{K}_{S}})$ is the stellar mass in solar units implied by the empirical relationship, and σ_m = 0.02. We incorporated the probability distribution for ${f}_{{ \mathcal M }}$ as a product within our total likelihood function (see below). We use this empirical relation from Mann et al. (2019) because it is precise and allows us to account for fit coefficient uncertainties and their correlations. Relative to other ${ \mathcal M }$–${M}_{{K}_{S}}$ relations from Delfosse et al. (2000) or Benedict et al. (2016), for early M dwarfs, the Mann et al. (2019) relation yields smaller masses at fixed ${M}_{{K}_{S}}$. The effect is small (∼5%) but systematic; see the discussion within Mann et al. (2019).

Mass–Radius. For typical isolated stars in the field, we cannot determine their radii unless they are sufficiently close for interferometric measurements (e.g., Boyajian et al. 2012; von Braun et al. 2014). As such, to estimate radii, we employed the empirical mass–radius relationship developed in Section 3 using the fraction scatter model. Like the mass–luminosity relation for ${f}_{{ \mathcal M }}$ above, we summarize the mass–radius results (see Section 3) defining the normally distributed random variable ${f}_{{ \mathcal R }}$ as a function of ${ \mathcal M }$ and ${ \mathcal R }$:

$$\begin{eqnarray}&&{r}_{m}({ \mathcal M })={b}_{0}+{b}_{1}{ \mathcal M },\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{f}_{{ \mathcal R }}({ \mathcal M },{ \mathcal R })={ \mathcal R }-{r}_{m}({ \mathcal M }),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{f}_{{ \mathcal R }}({ \mathcal M },{ \mathcal R })\sim {{ \mathcal N }}_{{ \mathcal R }}(0,{\sigma }_{{ \mathcal R }}({ \mathcal M })),\end{eqnarray} \tag{ 10 }$$

where ${r}_{m}({ \mathcal M })$ is the stellar radius implied by the empirical relations at fixed mass, and ${\sigma }_{{ \mathcal R }}=0.031\times {r}_{m}({ \mathcal M })$, using the best value of the scatter parameter. We also sample the coefficients b₀ and b₁ from their joint posterior to account for the correlated fit parameter uncertainties. ⁷ Although this relationship was determined from eclipsing binary measurements, as shown in Figure 5 our result is consistent with the measured interferometric radii of single M dwarfs.

Color–Bolometric correction. As an additional empirical result, we used bolometric correction as a function of color. This provides a way to constrain the stellar bolometric luminosity from widely available photometry, such as the Johnson V-band or SDSS r-band. We used the calibration published by Mann et al. (2015) for J-band bolometric correction as a function of V − J or r − J, the two expressions from their work that provide the best precision. ⁸ Their results are expressed, following Section 4.1, as

$$\begin{eqnarray}&&{ \mathcal C }=(V-J)\,\mathrm{or}\,(r-J),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{ \mathcal B }({ \mathcal C })=\sum _{i=0}^{2\,\mathrm{or}\,3}{c}_{i}\,{{ \mathcal C }}^{i},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{f}_{L}({L}_{\mathrm{bol}},{M}_{J},{ \mathcal C })=-2.5\mathrm{log}[{L}_{\mathrm{bol}}/{L}_{0}]-{M}_{J}-{ \mathcal B }({ \mathcal C }),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{f}_{L}({L}_{\mathrm{bol}},{M}_{J},{ \mathcal C })\sim {{ \mathcal N }}_{L}(0,{\sigma }_{L}),\end{eqnarray} \tag{ 14 }$$

where σ_L = 0.016, and L₀ = 3.0128 × 10³⁵ erg s⁻¹ is the luminosity at zero bolometric magnitude. Mann et al. (2015) also demonstrated systemic residuals as a function of metallicity, [Fe/H]. Thus, for objects that may deviate significantly from solar metallicity, we included a term in Equation (12) that is dependent on [Fe/H], which reduces the scatter to σ_L = 0.012 (see Table 3 of Mann et al. 2015). ⁹ However, this improved precision is lost if the [Fe/H] uncertainty is large.

Likelihood. In general, we expressed Bayes' theorem for the probability distribution describing our parameters of interest as

$$\begin{eqnarray}&&{ \mathcal P }(\alpha \,| \,{ \mathcal D })\propto { \mathcal L }({ \mathcal D }\,| \,\alpha )\,p(\alpha ),\end{eqnarray} \tag{ 15 }$$

where ${ \mathcal P }$ is the posterior distribution for the set of parameters $\alpha =\{{ \mathcal M },{ \mathcal R },{L}_{\mathrm{bol}},d,{M}_{{K}_{S}},{M}_{J},[\mathrm{Fe}/{\rm{H}}],{ \mathcal C }\}$ given the observational constraints from the data, ${ \mathcal D }$, for the likelihood function, ${ \mathcal L }$. By sampling α, we can generate the observable constraints; for example, inverting the distance, d, to obtain the observed parallax, or combining d and M_J to obtain the observed J magnitude. Here, p(α) is the prior on these parameters, taken to be uniformly uninformative within the range of properties for which the empirical relations are applicable. Although [Fe/H] and ${ \mathcal C }$ are randomly sampled, they are used to connect the observables to the physical stellar properties through the empirical relations as needed. They are in effect inputs, and are not further constrained by our likelihood function. For example, when used, the [Fe/H] or color, ${ \mathcal C }$, are simply drawn from their observed distributions and used in the calibration for bolometric correction. The metallicity thus presents an instance in which spectroscopic data can provide an informative constraint.

The likelihood for the data is the product of the observational constraints, taken here to be normal distributions based on the observed values and uncertainties:

$$\begin{eqnarray}&&{\pi }_{\mathrm{obs}}\,| \,{\sigma }_{\pi },d\sim {{ \mathcal N }}_{\pi }(\pi (d),{\sigma }_{\pi }),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{J}_{\mathrm{obs}}\,| \,{\sigma }_{J},d,{M}_{J}\sim {{ \mathcal N }}_{J}(J(d,{M}_{J}),{\sigma }_{J}),\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{K}_{\mathrm{obs}}\,| \,{\sigma }_{{K}_{S}},d,{M}_{{K}_{S}}\sim {{ \mathcal N }}_{K}(K(d,{M}_{{K}_{S}}),{\sigma }_{{K}_{S}}),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{ \mathcal L }({ \mathcal D }| \alpha )={{ \mathcal N }}_{\pi }{{ \mathcal N }}_{J}{{ \mathcal N }}_{K}\prod _{j\in \{{ \mathcal M },{ \mathcal R },L\}}{N}_{j},\end{eqnarray} \tag{ 19 }$$

where the means of the distributions are determined from the set of random variables, α. We then multiply these by the distributions N_j representing the empirical relations for ${f}_{{ \mathcal M }}$, ${f}_{{ \mathcal R }}$, and f_L given in Equations (7), (10), and (14). Additional empirical relations can be employed in the same manner, assuming they are independent, as new multiplicative factors in Equation (19). We also sampled the distributions for the fit coefficients in Equations (5) and (8) (mass–luminosity and mass–radius relations, respectively) to account for that uncertainty. The former is made possible by the publicly available posterior from Mann et al. (2019), ¹⁰ and for the latter we sample from the posterior reported in Section 3. The bolometric correction relation does not have published coefficient uncertainties, and we can only include the scatter on the relation (Mann et al. 2015).

In the most general case, we employed all three empirical relations. To summarize this process: in effect, the parallax and K_S -band magnitude constrain the mass with the ${ \mathcal M }$–${M}_{{K}_{S}}$ relation, which in turn through the ${ \mathcal M }$–${ \mathcal R }$ relation constrains the radius. The color, [Fe/H], and J-band magnitude plus parallax constrain the bolometric luminosity, which combined with the radius estimate defines the effective temperature. By combining all the relations simultaneously, we can account for correlations between the parameters and properly include all of the uncertainty in observables (e.g., the parallax measurement is used at multiple stages) and empirical relations in our final parameter estimates.

However, for some of the targets in the literature, there are additional observational constraints that we can use. For example, certain targets have measurements of their bolometric flux, in which case we omitted Equation (14) for the bolometric correction calibration and simply used the observed normally distributed F_bol = L_bol/4π d² to constrain the luminosity. Similarly, measured angular diameters, θ = 2R/d, from interferometry (e.g., Boyajian et al. 2012), directly constrain the radius, and through the ${ \mathcal M }$–${ \mathcal R }$ relation refine the mass estimate. In Table 5, we display the set of observational constraints used for each field star. The use of our Bayesian framework allows us to flexibly incorporate these additional observational constraints. ¹¹

We sampled our posterior distribution (Equation (15)) within PyMC3 (Salvatier et al. 2016) using the Hamiltonian NUTS sampler, similarly to the methods discussed in Section 3. We then applied this process to the sample of field MUSS stars shown in Table 1, yielding marginalized parameter distributions and uncertainties listed in Table 3. We show an example for one of our targets, GJ 49, of the joint posterior distributions for our main stellar properties of interest in Figure 8. ¹² The results for each star are qualitatively identical, with more precise parameters available for those objects with observed bolometric fluxes and/or angular diameters.

Zoom In Zoom Out Reset image size

Table 3. Physical Properties of MUSS Field Stars ^a

Name b	Lbol (1031 erg s−1)	Mass (M⊙)	Radius (R⊙)	Teff (K)
MUSCLES
Prox. Cen.*	$0.5997{\pm }_{0.0075}^{0.0076}$	0.123 ± 0.003	0.147 ± 0.005	$2992{\pm }_{47}^{49}$
GJ 1061	0.628 ± 0.014	0.125 ± 0.003	0.152 ± 0.007	$2977{\pm }_{69}^{72}$
GJ 1214	$1.343{\pm }_{0.034}^{0.036}$	0.181 ± 0.005	$0.204{\pm }_{0.0084}^{0.0085}$	$3111{\pm }_{66}^{69}$
GJ 628*	$4.209{\pm }_{0.038}^{0.037}$	0.304 ± 0.007	0.319 ± 0.007	$3307{\pm }_{36}^{38}$
GJ 581*	4.562 ± 0.037	0.307 ± 0.007	0.310 ± 0.008	$3424{\pm }_{42}^{43}$
GJ 667C	5.51 ± 0.13	0.327 ± 0.008	0.337 ± 0.014	$3443{\pm }_{71}^{75}$
GJ 725A*	$5.783{\pm }_{0.068}^{0.069}$	0.336 ± 0.007	0.354 ± 0.003	$3401{\pm }_{17}^{18}$
GJ 876*	5.012 ± 0.040	0.346 ± 0.007	0.372 ± 0.004	$3201{\pm }_{19}^{20}$
GJ 436*	9.43 ± 0.11	0.425 ± 0.009	0.432 ± 0.011	$3477{\pm }_{44}^{46}$
GJ 832	10.56 ± 0.34	0.441 ± 0.011	0.442 ± 0.018	$3539{\pm }_{74}^{79}$
GJ 887*	14.08 ± 0.22	$0.479{\pm }_{0.010}^{0.011}$	0.474 ± 0.008	$3672{\pm }_{34}^{36}$
GJ 176*	13.46 ± 0.12	0.485 ± 0.012	0.474 ± 0.015	$3632{\pm }_{56}^{58}$
Mega-MUSCLES
TRAPPIST-1	$0.234{\pm }_{0.008}^{0.009}$	$0.090{\pm }_{0.002}^{0.003}$	0.120 ± 0.006	$2619{\pm }_{66}^{71}$
LP 756-18	$0.993{\pm }_{0.026}^{0.027}$	0.145 ± 0.004	0.170 ± 0.007	$3155{\pm }_{68}^{71}$
LHS 2686	1.078 ± 0.024	0.157 ± 0.004	0.182 ± 0.008	$3119{\pm }_{68}^{70}$
GJ 699*	$1.302{\pm }_{0.023}^{0.024}$	$0.1610{\pm }_{0.0035}^{0.0036}$	0.187 ± 0.001	3223 ± 17
GJ 729	1.537 ± 0.018	0.177 ± 0.004	0.200 ± 0.008	$3248{\pm }_{66}^{68}$
GJ 1132	$1.668{\pm }_{0.047}^{0.049}$	0.194 ± 0.005	0.215 ± 0.009	$3196{\pm }_{70}^{71}$
L 980-5	$2.488{\pm }_{0.078}^{0.079}$	0.232 ± 0.006	0.250 ± 0.010	$3278{\pm }_{70}^{74}$
GJ 674*	6.03 ± 0.14	0.353 ± 0.008	$0.361{\pm }_{0.011}^{0.012}$	$3404{\pm }_{57}^{59}$
GJ 15A*	8.608 ± 0.069	$0.393{\pm }_{0.008}^{0.009}$	0.385 ± 0.002	$3601{\pm }_{11}^{12}$
GJ 163	8.28 ± 0.24	0.405 ± 0.010	$0.409{\pm }_{0.016}^{0.017}$	$3460{\pm }_{74}^{76}$
GJ 849	$11.051{\pm }_{0.094}^{0.095}$	0.465 ± 0.011	0.464 ± 0.018	$3492{\pm }_{68}^{70}$
GJ 649*	16.74 ± 0.170	0.524 ± 0.012	0.531 ± 0.012	$3621{\pm }_{40}^{41}$
GJ 676A	$34.04{\pm }_{0.80}^{0.84}$	0.631 ± 0.017	$0.617{\pm }_{0.027}^{0.028}$	$4014{\pm }_{90}^{94}$
FUMES
G 249-11	2.587 ± 0.069	0.237 ± 0.006	0.255 ± 0.010	$3277{\pm }_{68}^{70}$
GJ 4334	$3.680{\pm }_{0.088}^{0.090}$	0.294 ± 0.007	0.307 ± 0.012	$3260{\pm }_{67}^{69}$
LP 55-41	8.10 ± 0.22	$0.413{\pm }_{0.010}^{0.011}$	0.416 ± 0.017	$3412{\pm }_{73}^{75}$
LP 247-13	12.68 ± 0.35	0.495 ± 0.013	$0.492{\pm }_{0.020}^{0.021}$	$3511{\pm }_{76}^{79}$
GJ 49 c	17.84 ± 0.44	$0.541{\pm }_{0.015}^{0.014}$	0.534 ± 0.023	$3670{\pm }_{80}^{86}$
GJ 410	$21.34{\pm }_{0.51}^{0.52}$	0.557 ± 0.015	0.549 ± 0.024	$3786{\pm }_{82}^{89}$
HAZMAT
GJ 173	12.57 ± 0.12	$0.470{\pm }_{0.011}^{0.012}$	0.469 ± 0.019	$3589{\pm }_{71}^{75}$
G 75-55	$22.28{\pm }_{0.56}^{0.58}$	0.560 ± 0.015	0.552 ± 0.024	$3818{\pm }_{85}^{87}$
Living with a Red Dwarf
GJ 213	$2.407{\pm }_{0.035}^{0.036}$	0.218 ± 0.005	0.238 ± 0.009	$3334{\pm }_{63}^{67}$
GJ 821	$6.977{\pm }_{0.087}^{0.089}$	0.355 ± 0.009	$0.363{\pm }_{0.014}^{0.015}$	$3518{\pm }_{69}^{73}$
Miscellaneous
LHS 2065	0.1271 ± 0.0007	0.082 ± 0.002	0.113 ± 0.006	$2317{\pm }_{56}^{61}$
VB 10	0.1911 ± 0.0013	$0.0881{\pm }_{0.0024}^{0.0026}$	$0.1183{\pm }_{0.0057}^{0.0059}$	$2508{\pm }_{60}^{63}$
VB 8	0.2468 ± 0.0016	$0.0914{\pm }_{0.0025}^{0.0026}$	$0.1214{\pm }_{0.0057}^{0.0060}$	$2640{\pm }_{64}^{65}$
LHS 3003	0.230 ± 0.0019	$0.0929{\pm }_{0.0025}^{0.0027}$	0.123 ± 0.006	$2580{\pm }_{61}^{64}$
GJ 406*	$0.406{\pm }_{0.005}^{0.006}$	0.110 ± 0.003	0.144 ± 0.004	$2749{\pm }_{41}^{44}$
LHS 1140	$1.615{\pm }_{0.039}^{0.041}$	0.183 ± 0.005	0.205 ± 0.008	$3250{\pm }_{69}^{71}$
GJ 447*	1.401 ± 0.019	0.176 ± 0.004	0.198 ± 0.007	$3189{\pm }_{53}^{55}$
GJ 1207	2.586 ± 0.042	0.235 ± 0.006	0.253 ± 0.010	$3288{\pm }_{65}^{67}$
GJ 1243	2.729 ± 0.090	0.240 ± 0.006	0.258 ± 0.010	$3303{\pm }_{68}^{72}$
GJ 545	$3.198{\pm }_{0.055}^{0.056}$	0.253 ± 0.007	0.269 ± 0.011	$3362{\pm }_{68}^{70}$
GJ 3378	3.292 ± 0.062	0.263 ± 0.006	$0.278{\pm }_{0.010}^{0.011}$	$3332{\pm }_{64}^{66}$
GJ 403	$3.292{\pm }_{0.052}^{0.051}$	0.265 ± 0.007	0.280 ± 0.011	$3320{\pm }_{64}^{69}$
GJ 402	3.057 ± 0.032	$0.268{\pm }_{0.006}^{0.007}$	0.284 ± 0.011	$3240{\pm }_{60}^{65}$
GJ 3991	$3.270{\pm }_{0.043}^{0.044}$	$0.2714{\pm }_{0.0065}^{0.0066}$	0.286 ± 0.011	$3279{\pm }_{63}^{65}$
GJ 3325	$3.595{\pm }_{0.043}^{0.042}$	0.274 ± 0.007	0.289 ± 0.011	$3341{\pm }_{65}^{69}$
GJ 273*	$4.143{\pm }_{0.057}^{0.058}$	$0.2982{\pm }_{0.0064}^{0.0067}$	0.319 ± 0.004	$3297{\pm }_{21}^{22}$
GJ 4367	$4.114{\pm }_{0.089}^{0.087}$	0.304 ± 0.008	0.316 ± 0.012	$3305{\pm }_{65}^{69}$
YZ CMi	$4.35{\pm }_{0.14}^{0.15}$	0.316 ± 0.008	0.328 ± 0.013	$3293{\pm }_{71}^{74}$
EV Lac	4.88 ± 0.15	0.320 ± 0.008	0.331 ± 0.013	$3370{\pm }_{71}^{75}$
GJ 411*	8.399 ± 0.081	0.389 ± 0.008	0.392 ± 0.004	3547 ± 18
GJ 687*	8.330 ± 0.070	0.411 ± 0.009	$0.4187{\pm }_{0.0063}^{0.0066}$	$3426{\pm }_{28}^{27}$
AD Leo	8.67 ± 0.10	0.426 ± 0.010	$0.428{\pm }_{0.016}^{0.017}$	3425 ± 68
GJ 588	$10.95{\pm }_{0.31}^{0.32}$	0.461 ± 0.012	0.460 ± 0.019	$3499{\pm }_{73}^{79}$
GJ 526*	14.18 ± 0.19	0.488 ± 0.011	0.487 ± 0.008	$3628{\pm }_{31}^{32}$
GJ 3470	$15.20{\pm }_{0.41}^{0.43}$	0.502 ± 0.013	0.499 ± 0.021	$3649{\pm }_{78}^{82}$
Kepler-138	$21.78{\pm }_{0.53}^{0.56}$	0.545 ± 0.014	$0.537{\pm }_{0.022}^{0.023}$	$3846{\pm }_{84}^{86}$
GJ 809*	$19.85{\pm }_{0.21}^{0.22}$	0.547 ± 0.012	0.546 ± 0.006	$3729{\pm }_{23}^{22}$
K2-3	$21.45{\pm }_{0.62}^{0.64}$	0.551 ± 0.015	$0.543{\pm }_{0.023}^{0.024}$	$3810{\pm }_{84}^{89}$
GJ 880*	$20.00{\pm }_{0.15}^{0.16}$	0.553 ± 0.012	0.549 ± 0.003	3724 ± 12
GJ 205*	24.65 ± 0.22	$0.590{\pm }_{0.013}^{0.014}$	$0.5780{\pm }_{0.0025}^{0.0024}$	$3824{\pm }_{11}^{12}$
GJ 338A	$29.22{\pm }_{0.28}^{0.27}$	0.599 ± 0.016	$0.587{\pm }_{0.024}^{0.025}$	$3959{\pm }_{80}^{86}$
TOI-1235	$33.51{\pm }_{0.78}^{0.77}$	0.617 ± 0.016	0.604 ± 0.025	$4041{\pm }_{85}^{89}$

Notes.

^a Stellar properties are as determined in Section 4, quoting medians and the central 68% confidence interval. Joint posteriors for parameters are available upon request. ^b Star names listed with an asterisk have measured angular diameters from interferometric measurements; see Table 5. ^c For GJ 49, we updated the input optical photometry used in the estimation relative to those used in Pineda et al. (2021), leading to a slight change in luminosity and temperature.

Download table as:  ASCIITypeset images: 1 2

4.3. Results and Applications

4.3.1. Comparison of Stellar Properties with the Literature

By applying our methodology to the MUSS field stars, we have consistently determined their properties with well-defined uncertainties and their correlations. To test how our methods compare against typical literature determinations, we have compiled the stellar parameters for the MUSS stars found in eclectic assessments across numerous works. We typically utilized estimates from the same work, using the radii and effective temperatures to set the bolometric luminosity. These effective temperatures were often spectroscopically determined (atmospheric model based or otherwise), with radii coming from semi-empirical relations like that of Mann et al. (2015). The literature masses were often determined from earlier iterations of the ${ \mathcal M }$–${M}_{{K}_{S}}$ relation. The literature sources are listed in Appendix D.

We provide this comparison in Figure 9, showing both luminosity versus effective temperature and effective temperature versus mass. Often the estimates, like those of the interferometric objects, are nearly identical. Nevertheless, in both planes, our best-fit stellar parameters provide ensemble M-dwarf sequences with reduced scatter. As not every work treated all these properties, the compilation as is typical in the literature can be inconsistent, leading to larger outliers.

Zoom In Zoom Out Reset image size

4.3.2. M-dwarf Main Sequences

The MUSS field targets span the entire M-dwarf regime, and with consistently determined properties, we now examine the main sequence of these low-mass stars. Using our stellar parameters with well-defined uncertainties, we applied Gaussian Process (GP) regression models to the ensemble to capture the shape of the stellar main sequence across different parameters: L_bol versus T_eff, L_bol versus ${ \mathcal M }$, and T_eff versus ${ \mathcal M }$. GP models enable a flexible means of fitting, and have been used extensively to fit data of unknown functional form (e.g., Aigrain et al. 2016). We used a Matern52 covariance kernel, k, and fit the hyperparameters of scale length, l, covariance amplitude, η, and added an intrinsic noise, σ, at fixed abscissa to the data set. Specifically, we used the covariance function

$$\begin{eqnarray}\begin{array}{rcl}k(x,x^{\prime} ) & = & {\eta }^{2}\left(1+\displaystyle \frac{\sqrt{5{\left(x-x^{\prime} \right)}^{2}}}{l}+\displaystyle \frac{5{\left(x-x^{\prime} \right)}^{2}}{3l}\right)\\ & & \times \exp \left(-\displaystyle \frac{\sqrt{5{\left(x-x^{\prime} \right)}^{2}}}{l}\right).\end{array}\end{eqnarray} \tag{ 20 }$$

This choice of kernel is general-purpose, flexible, and yields GP models that are smooth without being too smooth (i.e., infinitely differentiable).

In fitting a particular sequence, for example, L_bol–T_eff, the data set consisted of bivariate normal distributions for each star, determined by the covariance and medians of the joint posteriors for the parameters from the Monte Carlo sampling explained in Section 4.2. This allows our methods to both account for uncertainties in both dimensions in the fitting and include the correlations of those uncertainties. The fitting was done through Monte Carlo sampling within PyMC3. We present the results of these methods in Figures 10–12.

Zoom In Zoom Out Reset image size

We illustrate the traditional theoretical HR diagram with luminosity as a function of effective temperature in Figure 10. In the left panel, the data ensemble shows the shape of this main sequence, with the parameter uncertainties illustrated as 1σ error ellipses, conveying the size and correlation of the output parameters, with shading corresponding to the stellar radius. The bulk of the data pertain to early and mid-M dwarfs, with a handful of late-M dwarfs extending the sequence to effective temperatures as low as 2300 K. For comparison, we also include the 5 Gyr isochrones from the MESA Isochrones and Stellar Tracks (MIST; Dotter 2016). We show the results of the GP regression of L_bol–T_eff in the right panel of Figure 10, with data points displaying the central medians, and the best-fit GP curves sampled from the posterior distribution. The ensemble of lines presents the best GP solutions, with the residual plot showing the best-fit scatter about that solution. The result is smooth and fits the data well with an intrinsic scatter of ∼0.1 dex.

In contrast to L_bol–T_eff, the L_bol–${ \mathcal M }$ sequence shown in Figure 11 is very well-defined, with a tighter scatter about the mean GP regression model. Figure 11 is similar to Figure 10, but the data ellipses in the left panel now have shading to visualize the target effective temperatures. The sample posterior for the fit luminosity scatter peaks at ∼0.012 dex, but the distribution on the scatter extends down toward zero, because the data point uncertainties can largely account for the bulk of the scatter in the model regression. Accordingly, we describe the sampled scatter in luminosity at fixed mass with a 3σ upper limit of 0.03 dex in Table 4. The differences between the L_bol–T_eff and the L_bol–${ \mathcal M }$ sequences are largely a consequence of the mid-M dwarfs. Objects around 3300 K show a broad range of possible luminosities, whereas in L_bol–${ \mathcal M }$ space, the relationship exhibits far less scatter.

Zoom In Zoom Out Reset image size

Table 4. GP M-dwarf Main-sequence Models ^a

Variables b	l	η	σ
L–Teff	${ \mathcal H }{ \mathcal C }(100)$ K	${ \mathcal H }{ \mathcal C }(3)$ dex	${ \mathcal H }{ \mathcal C }(0.4)$ dex
	$2260{\pm }_{550}^{800}$	$2.4{\pm }_{0.8}^{1.7}$	$0.096{\pm }_{0.012}^{0.014}$
L–${ \mathcal M }$	${ \mathcal H }{ \mathcal C }(0.3)$ M⊙	${ \mathcal H }{ \mathcal C }(3)$ dex	${ \mathcal H }{ \mathcal C }(0.1)$ dex c
	$0.55{\pm }_{0.15}^{0.23}$	$4.6{\pm }_{1.8}^{4.0}$	<0.03
Teff–${ \mathcal M }$	${ \mathcal H }{ \mathcal C }(0.1)$ M⊙	${ \mathcal H }{ \mathcal C }(50)$ K	${ \mathcal H }{ \mathcal C }(50)$ K
	$0.57{\pm }_{0.14}^{0.21}$ M⊙	$3860{\pm }_{1430}^{3300}$ K	$57{\pm }_{9}^{10}$ K

Notes.

^a The rows for each variable pair show first the model priors, and underneath the median and central 68% confidence intervals of the fit parameters. ${ \mathcal H }{ \mathcal C }(\beta )$ indicates a half-Cauchy distribution with shape parameter β. ^b In practice, for the luminosity sequences, we fit in the space of log₁₀ of the bolometric luminosity relative to solar. ^c For the mass–luminosity sequence, we additionally bounded the scatter to be greater than 0.001 dex, to prevent the small scatter from causing the likelihood to diverge.

Download table as:  ASCIITypeset image

We show this further in the T_eff–${ \mathcal M }$ relationship of Figure 12, with the shading of the error ellipses in the left panel indicative of the stellar radius. Among the very low-mass stars, the higher-mass M dwarfs exhibit much higher effective temperatures, with differences of ∼600 K between 0.1 and 0.2 M_⊙ stars. The effective temperature sequence then plateaus in the mid-M dwarfs before increasing again in the early-to-mid-M-dwarf regime. For the mid-M dwarfs, on average across the sample between about 0.2 and 0.3 solar masses, the effective temperature is ostensibly constant. Furthermore, the GP regression model in T_eff–${ \mathcal M }$ reveals an intrinsic scatter of only ∼60 K, at fixed mass. The results of all our GP regression models are summarized in Table 4, and we further publicly provide the GP solutions with 68% confidence curves (see footnote 5).

Zoom In Zoom Out Reset image size

Collectively, these results show how, in M dwarfs, especially the mid-M dwarfs, effective temperature, when used to predict other stellar properties, is subject to a great deal of scatter. Although T_eff is a useful quantity and directly related to the inputs to model spectra, its use should, when possible, be tied to empirically determined temperatures, like those from interferometric measurements. This is especially true if the metallicity is unknown. Mass and luminosity are both more fundamental and show a very tight relationship, and are thus better suited as the basis for inferring additional stellar properties. Our sample averages over the possible metallicity differences, and being composed largely of near solar metallicity objects, is representative of the nearby stellar neighborhood. It is likely that the intrinsic scatter in the T_eff–${ \mathcal M }$ relation is defined by the distribution of metallicities in the MUSS field stars (see Figure 1).

In Figures 10–12, we also showed how the MIST isochrones at 5 Gyr compare to our empirical sequences. The L_bol–${ \mathcal M }$ plane shows that the solar metallicity track matches the bulk of the data quite well (the model does not extend to the lowest-mass stars in our sample). However, the model effective temperatures are systematically higher than the MUSS sample stars. With well-matched bolometric luminosities, this likely reflects the long-standing issue that the evolutionary models often underpredict M-dwarf stellar radii (e.g., López-Morales 2007; Kesseli et al. 2018) and thus need higher T_eff to match the luminosity.

4.3.3. An Empirical Boundary for Fully Convective Interiors?

There is an apparent plateau at low masses in the T_eff–${ \mathcal M }$ relationship (see Figure 12). Rabus et al. (2019) showed, using an ensemble of M dwarfs, that the ${ \mathcal R }$–T_eff relation appeared to have a "discontinuity" at effective temperatures between 3200 and 3340 K. Within this temperature regime, their relation was effectively vertical; we illustrate this feature in Figure 13 using the MUSS sample, and including the linear fits from Rabus et al. (2019) (their Equation (7)), not accounting for possible metallicity effects. The linear fits from Rabus et al. (2019) are consistent with our data.

Zoom In Zoom Out Reset image size

Given the largely one-to-one relationship between stellar mass and radius (see Section 3), recasting this relation as T_eff–${ \mathcal M }$ yields the results we have shown in Figure 12. The effective temperature plateau at low masses (∼0.25 M_⊙) corresponds exactly to this "discontinuity" in ${ \mathcal R }$–T_eff. Rabus et al. (2019) considered this feature within the M-dwarf population as indicative of the transition between partially and fully convective interiors, a claim we discuss further below. First, we consider the question: at what mass does this feature occur?

To identify the plateau feature, we determined where the derivative of the GP model approaches zero or is otherwise minimal. One benefit of the GP model is that not only can it predict the values of the functional relationship, but the derivative of the GP is itself a GP, and can therefore be used similarly to define the derivative of that function (Rasmussen & Williams 2006). Using the results of the regression model, we further calculate the derivative, $\partial {T}_{\mathrm{eff}}/\partial { \mathcal M }$, and show those results in Figure 14, sampling at ${\rm{\Delta }}{ \mathcal M }=0.001$ M_⊙. ¹³ The top panel replots the T_eff–${ \mathcal M }$ relation, with the bottom one showing its derivative. At the location of the plateau, the derivative indeed appears to dip to zero, but is otherwise positive.

Zoom In Zoom Out Reset image size

To define at which mass this feature exists, we took the ensemble of GP curves that represent the derivative, and defined a representative sequence from the median value of the derivative at each mass. The median derivative minimum occurs clearly at 0.256 M_⊙. To estimate a confidence interval, we restricted consideration of each individual GP derivative curve from the ensemble to masses below 0.35 M_⊙. Then, for each curve, we recorded within that range (${ \mathcal M }\lt 0.35$ M_⊙) where the curve was at its lowest point. Although individual curves can be variable, the ensemble produced a distribution of masses reflective of the plateau feature. The result with 68% confidence interval was ${{ \mathcal M }}_{p}=0.256{\pm }_{0.021}^{0.027}$ M_⊙. Along the median T_eff–${ \mathcal M }$ relation, this mass corresponds to 3305 K.

Rabus et al. (2019) suggested that this plateau was indicative of the transition between fully and partly convective interiors, with it taking place at masses just below that corresponding to the Gaia M-dwarf gap (∼0.35 M_⊙, Jao et al. 2018; Feiden et al. 2021). However, modeling of the impact of the ³He instability (van Saders & Pinsonneault 2012), the cause for the Jao Gap, indicates that the onset of this instability marks the transition between partly and fully convective interiors (Baraffe & Chabrier 2018; MacDonald & Gizis 2018; Feiden et al. 2021). The plateau mass (∼0.26 M_⊙) is too far from the masses corresponding to the Gaia M-dwarf gap to likely be associated with the convective transition.

It is unclear what effect may be leading to small or no changes in effective temperature between 0.2 and 0.3 M_⊙, although it may be related to the influence of electron degeneracy pressure in the core (Cassisi & Salaris 2019). For example, field age model solar metallicity isochrones do not exhibit such a shallow effective temperature derivative in this regime (e.g., Baraffe et al. 2015; Feiden 2016). We illustrate this by plotting our median GP curves against 5 Gyr MIST isochrones (Morton 2015; Dotter 2016) in Figure 15, zooming in on this mass regime around 0.3 M_⊙. Despite matching the L_bol–${ \mathcal M }$ sequence well, the model effective temperature sequence does not reproduce the extent of the plateau, and it generally overestimates T_eff. It is possible that the observed plateau feature may be an emergent effect of the stellar sample, which averages together a range of stellar metallicities (see Figure 1). However, the observed slope remains difficult to reproduce. Nevertheless, on average it appears observationally that the radial response to increasing luminosity over that mass range balances such that the effective temperature remains nearly constant.

Zoom In Zoom Out Reset image size

The desire to measure the plateau feature in the T_eff–${ \mathcal M }$ sequence prompted our derivative analysis. Interestingly, this derivative also reveals a local maximum at ∼0.35 M_⊙, which is easy to miss when simply examining the standard sequence (e.g., top panel of Figure 14). To define at which mass this local maximum takes place, we examined the median GP derivative, identifying the peak to be at 0.361 M_⊙. Considering the individual GP curves within the range 0.25–0.42 M_⊙ (roughly between the local derivative minima; see Figure 14), we recorded the location of the maximal value of the derivative curves, finding ${{ \mathcal M }}_{b}=0.361{\pm }_{0.024}^{0.036}$ M_⊙. This coincides with theoretical predictions for the solar metallicity mass boundary between partly and fully convective interiors in low-mass stars (e.g., Chabrier & Baraffe 2000). The corresponding temperature along the median T_eff–${ \mathcal M }$ sequence is 3414 K. These mass values, M_p and M_b , correspond to the temperatures between M3 and M4 spectral types (e.g., Boyajian et al. 2012; Pecaut & Mamajek 2013), where the convective transition has been thought to occur (Reid & Hawley 2005). We consider only the higher mass feature (M_b ) as indicative of the convective transition, as it is consistent with theoretical predictions.

To corroborate this finding, we also examined the derivative of the GP L_bol–${ \mathcal M }$ sequence. We show those results in Figure 16. The median derivative sequence shows a well-defined local maximum at 0.337 M_⊙, with the slope remaining positive throughout the masses corresponding to the temperature plateau. We defined the confidence interval for this feature from where the ensemble derivative curves were maximal on the interval 0.27–0.40 M_⊙, yielding a result of ${{ \mathcal M }}_{b}=0.337{\pm }_{0.026}^{0.013}$ M_⊙ (see Figure 17). Along the median L_bol–${ \mathcal M }$ sequence, this mass corresponds to a log luminosity of −1.843 relative to solar. This corresponds exactly to the stellar luminosity that near-solar-metallicity models predict will exhibit the ³He instability (Feiden et al. 2021)—the marker for the convective transition, and explanation for the Gaia M-dwarf gap.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although within uncertainties, this mass value is slightly smaller than that found in the T_eff–${ \mathcal M }$ sequence, possibly due to the effective temperature plane combining the effects on radius and luminosity across the convective transition (Baraffe & Chabrier 2018). In both luminosity and temperature planes, we see derivative undulations with changing stellar mass, indicating perturbations along the main stellar sequence. With these features occurring across the 0.3–0.4 M_⊙ range, we take the first feature (${{ \mathcal M }}_{b}=0.337{\pm }_{0.026}^{0.013}$ M_⊙) in the derivative sequence with increasing mass as indicative of the transition onset.

We consider the feature in the luminosity–mass sequence to be the more accurate empirical measure of the transition from fully to partially convective interiors, as this stellar sequence shows less scatter than the effective temperature–mass sequence. This boundary measurement not only agrees with theoretical predictions, its determination from the derivative of the stellar sequence, agrees well with the physical expectations of the ³He instability. MacDonald & Gizis (2018) note that the instability should coincide with local extremal changes in the luminosity, seen here as a local maximum in the derivative. Baraffe & Chabrier (2018) also showed model isochrones that illustrate this idea, with each isochrone curve exhibiting abrupt luminosity shifts at slightly different masses corresponding to the convective transition. These shifts can be difficult to discern, and they may not be totally evident in the models (see Figure 15). Our analysis illuminated these features as undulations in the sequence derivative. The exact location and width of the features we measured likely reflects the distribution of ages and metallicities of the MUSS field stars, and thus represents the solar neighborhood result for the empirical fully convective boundary mass of M-dwarf stars.

Well-defined empirical relationships for stellar properties are critical for many astrophysics investigations. The new mass–radius relationship presented in this article (Section 3) is one such example, and can be readily used to empirically assess masses and/or radii in the low-mass star regime. This and the other empirical relationships used in characterizing the MUSS will surely be supplanted as calibration data sets expand; however, the methods we have demonstrated (Section 4.1) will facilitate their future incorporation in Bayesian model fitting.

To further illustrate how: imagine a hypothetical joint Bayesian regression analysis of stellar radial velocity data, and an exoplanet transit light curve for some system. The host star mass and radius are among the parameters of such a fit. Without any additional information, our formulation for an empirical mass–radius relationship (e.g., Equation (10)) can be multiplied to the model likelihood to give stellar mass and radius solutions (and thus the planet properties that depend on them) entirely consistent with the empirical relationship. This could properly and fully incorporate all the sources of uncertainty, as well as their correlations. Additional constraints (e.g., stellar angular diameters) could provide improved precision that would naturally propagate through to all of the parameter posteriors.

Often in the literature, the approach to this problem does not fully account for the sources of known uncertainty, or ignores their correlations, potentially biasing results or mischaracterizing the actual error distributions. In particular, ignoring the correlations can lead to overestimated errors. These effects could lead to potentially missing important astrophysical trends, or to overconfidence in possible results. The solution often is to acquire more data, but better statistical characterization can more accurately direct significant astrophysical findings.

Consistently applying our methods further enables a clear understanding of potential underlying systematics, and provides a platform for ensemble studies of stellar properties, as we have demonstrated with the MUSS stellar sequences (Section 4.3.2). This GP regression is flexible and provides a powerful tool for examining the stellar sequences. Features along these curves can be difficult to discern, and are often elided through standard parametric regressions assuming a functional form. However, by employing the GP derivative predictions, we have revealed features that may otherwise have been missed. Using this analysis, we were able to define the plateau in the effective temperature–mass relation to be at ∼0.26 M_⊙. Most importantly, the derivatives in L_bol–${ \mathcal M }$ and T_eff–${ \mathcal M }$ showed undulations that were predicted by models of stellar interiors as a consequence of the transition between partly and fully convective interiors. We thus empirically defined that boundary to be at ∼0.34 M_⊙ with the solar neighborhood MUSS stars.

The exact measures of these features reflect average sequence properties for solar neighborhood M dwarfs, likely a consequence of age and metallicity distributions, but present a clear empirical result for the characteristics of nearby stars. These methods provide a powerful new way to examine empirical stellar sequences, for comparison to model predictions. Our results contribute to the empirical evidence of the convective transition in low-mass stellar interiors (e.g., Jao et al. 2018), and a first empirical measure of the corresponding transition mass.

Larger stellar samples with narrower metallicity distributions and accurately/consistently defined properties, as we have employed here with the MUSS stars, will enable more precise examinations of how ensemble stellar sequences can be used to probe the physics of stellar interiors.

We summarize our key results as follows:

• 1.  
  We have defined a new M-dwarf mass–radius relation that gives 3.1% uncertainties at fixed mass across a range of 0.09–0.7 M_⊙; see Section 3.
• 2.  
  We have developed a self-consistent Bayesian framework to simultaneously incorporate empirical relationships in Bayesian model fitting. Although general in nature, we illustrated its application in stellar parameter estimation, flexibly incorporating observational constraints; see Section 4.1.
• 3.  
  Combining multiple photometric empirical relations and simultaneously accounting for the known uncertainty across relations, we applied these methods to the MUSS field stars, consistently defining their ${ \mathcal M }$, ${ \mathcal R }$, T_eff, and L_bol; see Section 4.2. These stars are of interest to the community and these properties have already been used by the analyses of Pineda et al. (2021) and Youngblood et al. (2021).
• 4.  
  The MUSS results define stellar sequences, which we fit using GP regression models, now publicly available; see Section 4.3.2.
• 5.  
  We empirically assessed the boundary mass between partly and fully convective M dwarfs to be ${{ \mathcal M }}_{b}=0.337{\pm }_{0.026}^{0.013}$ M_⊙; see Section 4.3.3.
• 6.  
  A skew normal distribution provides a useful way to approximate published results reporting asymmetric error bars that preserves the central confidence interval and the centroid, for use as prior information in additional studies; see Appendix B.

The authors would like to thank Zachary Berta-Thompson, Elisabeth Newton, Girish Duvvuri, Aurora Kesseli, and Eileen Gonzales for useful discussions and commentary during the preparation of this work. The authors would also like to thank the anonymous referee for thoughtful comments and constructive feedback.

Support for Program numbers HST-GO 14640, 14633, and 15071 were provided by NASA through grants from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.

A.Y. acknowledges support by an appointment to the NASA Postdoctoral Program at Goddard Space Flight Center, administered by USRA through a contract with NASA.

This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

In Section 4, we discussed our methods for determining the stellar properties of the MUSS field stars. These methods relied on combining multiple empirical relations and available data for constraining the properties, usually photometry, and whenever available, measured bolometric fluxes or interferometric angular diameters. We show the specific data used in the parameter estimates of all field stars in Table 5. A few of these stars warranted additional commentary, which we provide below.

Kepler-138. For Kepler-138, there have been some inconsistent parameter estimations published in the literature. Almenara et al. (2018) discuss the history of these estimates and possible implications for the planetary system properties. Our new values are consistent with those used by Kipping et al. (2014) from Pineda et al. (2013), and show agreement with the analysis of Berger et al. (2020). For comparison with transit data, our inferred stellar density for Kepler-138 is $4.95{\pm }_{0.54}^{0.61}$ g cm⁻³, whereas Almenara et al. (2018) report a value of $3.92{\pm }_{0.66}^{0.81}$ g cm⁻³ from light-curve modeling, which is consistent within uncertainties. Their lower density value, however, led them to infer a larger stellar radius than we report, and in spectral energy distribution fitting thus infer a distance of 74.3 ± 5.8 pc, which is larger than the Gaia DR2 parallax inferred distance of 66.99 ± 0.11 pc. The smaller radius (higher density) value would move that estimate toward agreement. We additionally note that the spectrophotometric distance from Pineda et al. (2013) of 66.5 ± 7.3 pc is surprisingly close to the Gaia value, lending additional validation to those methods and to studies of Kepler-138 that have since used their parameter estimates.

TOI-1235. TOI-1235 (TYC 4384-1735-1) has recently been shown to host a transiting exoplanet near the gap in the exoplanet radius distribution (Bluhm et al. 2020; Cloutier & Menou 2020). Our analysis of the stellar properties of this star ($M=0.617{\pm }_{0.017}^{0.016}{M}_{\odot }$, R = 0.604 ± 0.027) give slightly smaller values in mass and radius than that used in those exoplanet papers, although generally within uncertainties. Comparing these values, the Cloutier & Menou (2020) work relied on Benedict et al. (2016), which gives systematically higher masses than our methods (see discussion within Mann et al. 2019). Both Cloutier & Menou (2020) and Bluhm et al. (2020) rely on spectroscopic effective temperatures tied to model atmospheres, and determine T_effs that are cooler than our more empirical approach. Given that our bolometric luminosity estimates are in strong agreement, this difference in assumed temperatures implies larger stellar radii. In general, our stellar parameters are in better agreement with those used by Bluhm et al. (2020) than Cloutier & Menou (2020), with our implied stellar density $3.95{\pm }_{0.44}^{0.51}$ g cm⁻³ in good agreement with the density from transit light-curve modeling ($3.74{\pm }_{0.31}^{0.30}$ g cm⁻³; Bluhm et al. 2020). Our smaller stellar radius implies a smaller planet, and systematically shifts the planet away from the gap in the exoplanet radius distribution.

TRAPPIST-1. Our properties for TRAPPIST-1 may be of particular interest. The bolometric flux is taken from Gonzales et al. (2019). Although not listed in that work explicitly, the value and uncertainties quoted in Table 5 were provided by them (E. Gonzales 2021, private communication). Our stellar parameter values are empirical (M = 94.1 ± 2.6 M_Jup, R = 1.09 ± 0.06 R_Jup), and agree remarkably well with their evolutionary-model-dependent approach (M = 90 ± 8 M_Jup, R = 1.16 ± 0.03 R_Jup).

Probabilistic modeling of astronomical phenomena has become a widespread practice across a variety of subfields, enabled by increased computational ability and efficient Monte Carlo methods for sampling likelihood functions. An outcome of this work is often probability distributions describing parameters of interest, usually quoted as a "best" value and 1σ error bars corresponding to the middle ∼68% confidence interval. These properties are then typically assumed to describe a normal distribution with mean at the "best" value and variance σ². In most instances, the probability distribution for the parameter of interest can indeed be approximated by a normal distribution and these results can be readily incorporated into subsequent studies as prior information. However, with typically nonlinear models and transformations between random variables of interest, results often do not yield symmetric, normally distributed posterior distributions. The 68% confidence interval is then denoted with unequal values, σ₋ and σ₊, below and above the "best" estimate, respectively. These results can then become problematic to use in new analyses. A common and unjustified practice is to average the uncertainties above and below separately when combining these data with other estimates (Barlow 2003). Another typical approach is to approximate the probability distribution as Gaussian using an estimate for the mean and averaging σ₋ and σ₊ to approximate the distribution standard deviation, giving an easy-to-use informative prior on the physical quantity of interest (D'Agostini 2004). This is reasonable when the posterior distribution is well-constrained, such that the log of the distribution around its peak can be approximated by a parabola. Furthermore, central limit theorem arguments make the normal distribution a good limiting function for use in many cases. As we show below, this practice, while giving reasonable results, can shift the central 68% confidence interval, as well as the mode of the distribution, potentially biasing the analysis, when an informative prior is incorporated within Bayesian inference modeling.

As an illustration of this effect, we consider two clearly asymmetric posterior distributions: the first one is from the literature, the parameter a₅ corresponding to the coefficient of the highest-order term in the fifth-order polynomial fit for the mass–luminosity relationship of Mann et al. (2019); the second is the marginalized posterior for the fractional scatter, ${\sigma }_{{ \mathcal R }}$, for the mass–radius relation defined in this work (Section 3). In Figure 18, we plot the histograms of these distributions for a₅ (left) and ${\sigma }_{{ \mathcal R }}$ (right) in gray, with two normal approximations as dashed lines. The line "Normal-1" plots a Gaussian distribution centered at the mean and standard deviation matching that of the posterior samples. The line "Normal-2" indicates a similar Gaussian distribution, but now centered at the peak of the posterior with standard deviation given by the average of σ₊ and σ₋ defined by the central 68% confidence interval of the posterior samples relative to the peak.

Zoom In Zoom Out Reset image size

Because the distributions are asymmetric, neither "Normal" approximation is ideal for conveying the probability density function of these variables. "Normal-1" shifts the peak of the distributions away from the most-probable value, whereas "Normal-2" retains that peak but places much greater probability mass away from the actual skew in the posterior distribution. Both approximations bias subsequent results away from the most-likely confidence interval defined by the actual posterior samples. If the posterior distribution for a given quantity of interest is available, it can be sampled directly, avoiding such issues. However, in practice this is rarely the case, and published values generally give only a "best" value with a 1σ confidence interval.

We consider here a method to estimate a useful probability density function from these "best" values that is Gaussian-like in shape but preserves the asymmetry of the quoted result. In effect, we draw upon the class of three parameter probability density functions that are perturbations to a standard normal distribution. In particular, we employ the skew normal distribution,

$$\begin{eqnarray}&&{ \mathcal S }(x)=\displaystyle \frac{2}{\sigma }\phi \left(\displaystyle \frac{x-\mu }{\sigma }\right){\rm{\Phi }}\left(\alpha \displaystyle \frac{x-\mu }{\sigma }\right),\end{eqnarray} \tag{ 21 }$$

where ϕ is the standard normal distribution, Φ is the cumulative distribution function for ϕ, μ is a location parameter, σ sets the scale, and α is the shape parameter. In the limit of α → 0, the distribution matches that of a Gaussian and becomes increasingly asymmetric for increasing ∣α∣, with the limit of ∣α∣ → ∞ yielding a half-normal distribution. The sign of α is also indicative of the direction of skew. For example, for α > 0, the skew normal distribution exhibits more probability mass at values greater than the mean (i.e., to the right), while the opposite is true for α < 0.

In Figure 18, we show that a skew normal approximation (solid line) does a remarkably good job of matching the shape of the posterior distributions for a₅ and ${\sigma }_{{ \mathcal R }}$: the kernel density estimate (KDE, dotted) and the skew normal curve for both parameters largely overlap, matching the envelope of the histogram.

The skew normal approximation is determined by numerically solving the system of equations for the three variables μ, σ, and α, such that the resulting skew normal distribution yields a central 68% confidence interval and median (or mode) that match the desired input "best" value with its asymmetric error bars. Using the cumulative function for the skew normal distribution, the system can be expressed as

$$\begin{eqnarray}&&{\int }_{-\infty }^{{x}_{u}}{ \mathcal S }(x;\mu ,\sigma ,\alpha )\,{dx}-0.84=0,\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{\int }_{-\infty }^{{x}_{l}}{ \mathcal S }(x;\mu ,\sigma ,\alpha )\,{dx}-0.16=0,\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{\int }_{-\infty }^{{x}_{m}}{ \mathcal S }(x;\mu ,\sigma ,\alpha )\,{dx}-0.50=0,\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\alpha \phi \left(\alpha \,\displaystyle \frac{{x}_{p}-\mu }{\sigma }\right)-{\rm{\Phi }}\left(\alpha \displaystyle \frac{{x}_{p}-\mu }{\sigma }\right)\displaystyle \frac{{x}_{p}-\mu }{\sigma }=0,\end{eqnarray} \tag{ 25 }$$

where we use only one of either Equation (24) or (25), depending on whether the central value in question is the median or mode, respectively, x_m is that input median value, x_p is the desired peak of the output skew normal distribution, x_u is the upper value on the interval bound (e.g., x_m + σ₊), and x_l is the lower value on the interval bound (e.g., x_m − σ₋).

We provide code for determining this numerical approximation in an online repository (see footnote 5). This approach has the advantages of formally matching the central confidence interval for the quantity of interest and giving a probability density function skewed in the appropriate direction, and it can be applied to any result publishing asymmetric error bars as long as it is clear what defines the central quantity (either a median or mode). The skew normal distribution is already implemented in commonly used software, including Python's scipy.

This approach remains an approximation and should not replace directly sampling the posterior when available, although there may be use cases for which having a well-defined functional form may be advantageous. In practice, while we only demonstrate its applicability for two example posteriors, this method should work well for distributions with Gaussian tails and only slight deviations from symmetric. A skew normal distribution constructed in this way will preserve the range of the central confidence interval and likely match the proper median (or mode) of the parameter in question. By definition, this will provide better approximations to asymmetric distributions that are assumed to be Gaussian-like rather than assuming a symmetric Gaussian distribution. Sample distributions with long tails that depart more significantly from a normal distribution will consequently not be approximated well with a skew normal distribution, although a similar approach could be taken with an exponential Gaussian distribution or similar three-parameter probability density functions.

We note throughout this article that ignoring the correlations in the uncertainties enlarges the resulting error estimates when using empirical fit relationships, for example, from published mass–radius relations for M dwarfs. Although the correlations are not often published, typical polynomial fits of any kind, when applied to data, will yield parameters that have correlated uncertainties. In general, the severity of ignoring the correlations when using such relationships will depend on the details of the polynomial relation and its parameter uncertainties. Nevertheless, that process will always yield errors larger than implied by the underlying data. This is because, by ignoring the correlations, one allows the fit parameter uncertainties to include parameter space that is actually strongly disfavored by the regression analysis. Here, we provide a specific example using the results from Section 3, for the magnitude of this effect.

If we consider a 0.6 M_⊙ star, its implied radius using our fractional scatter ${ \mathcal M }$–${ \mathcal R }$ relation (Section 3) is $0.588{\pm }_{0.019}^{0.018}$ R_⊙. If we ignore the correlations of the fit parameter uncertainties, the implied result becomes $0.588{\pm }_{0.020}^{0.021}$ R_⊙. That corresponds to a ∼10% increase in the standard deviation. Similarly, for a 0.1 M_⊙ star, the corresponding results would be 0.129 ± 0.005 R_⊙ and 0.129 ± 0.006 R_⊙ for including and excluding the uncertainty correlations, respectively. That corresponds to a ∼20% increase in the standard deviation. In these cases, the difference is not especially dramatic, because the uncertainties are dominated by the intrinsic scatter term. However, when just considering the contribution from the polynomial parameter uncertainties, the implied errors are at least twice as large if the correlations are ignored.

Whether these effects are very significant depends on the science question. Nevertheless, it is straightforward to account for uncertainty correlations, and doing so enables more accurate results and reliable probabilistic calculations.

In Section 4.3.1, we compared our new stellar parameter determinations with typical values found in the literature. We collected the masses, radii, and effective temperatures for that comparison from the following list of references: Bonfils et al. (2005), Demory et al. (2009), Houdebine (2010), Berta et al. (2011), Kraus et al. (2011), von Braun et al. (2011), Boyajian et al. (2012), Rojas-Ayala et al. (2012), von Braun et al. (2012), Anglada-Escudé et al. (2013), Lépine et al. (2013), Kordopatis et al. (2013), Pecaut & Mamajek (2013), Pineda et al. (2013), Rajpurohit et al. (2013), Santos et al. (2013), Tuomi & Anglada-Escudé (2013), Gaidos et al. (2014), Neves et al. (2014), von Braun et al. (2014), Berta-Thompson et al. (2015), Mann et al. (2015), Terrien et al. (2015), Awiphan et al. (2016), Gillon et al. (2016), Houdebine et al. (2016), Kopytova et al. (2016), Astudillo-Defru et al. (2017), Dressing et al. (2017), Newton et al. (2017), Houdebine et al. (2017), Berger et al. (2018), Damasso et al. (2018), Van Grootel et al. (2018), Houdebine et al. (2019), Rabus et al. (2019), and Jeffers et al. (2020).