Chemical Properties of the Local Disk and Halo. II.  Abundances of  3745 M dwarfs and Subdwarfs from Improved Model Fitting of Low-Resolution Spectra

Neda Hejazi
Department of Physics and Astronomy, Georgia State University, Atlanta, GA 30303, USA

Sébastien Lépine
Department of Physics and Astronomy, Georgia State University, Atlanta, GA 30303, USA

Thomas Nordlander
Research School of Astronomy & Astrophysics, Australian National University, Canberra, ACT 2611, Australia
The ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions, Canberra, ACT 2611, Australia

We present a  model-fit pipeline to determine the stellar parameters of M-type dwarfs, which is an improvement upon our previous work described in Hejazi et al. 2020.  We apply this pipeline to analyze the low-resolution (R$\sim$2000) spectra of 3745 M dwarfs/subdwarfs, collected at the MDM Observatory, Lick Observatory, Kitt Peak National Observatory, and Cerro Tololo Interamerican Observatory. We examine the variation of the inferred parameter values in the HR diagram constructed from their Gaia Early Data Release 3 (EDR3) parallaxes and optical magnitudes. We also study the distribution of our stars in the abundance diagram of [$\alpha$/Fe] versus [M/H]  and inspect the variation of their metallicity class, effective temperature, and surface gravity as well as their Galactic velocity components U, V, and W in this diagram. In addition, the analyses of the stars’ projected motions in the two-dimensional UV, VW, and UW planes and the variation of their chemical parameters in these planes, and also their distribution in the abundance-velocity diagrams are important parts of this study. The precision of our model-fit pipeline is confirmed by the clear stratification of effective temperature and chemical parameters in the HR diagram, the similarity of  the stars’ distribution in the  [$\alpha$/Fe] versus [M/H] diagram and in the metallicity-velocity planes with those from other studies, revealing substructure in the abundance-velocity diagrams, and  chemical homogeneity between the components of a set of binary systems.

§ INTRODUCTION

M-type dwarfs dominate the stellar population of our Galaxy, including $\sim$ 70$\%$ of all stars by number (Reid & Gizis 1997; Bochanski et al. 2010) and hold enormous potential for being used as tools to probe  the structure and chemo-dynamical evolution of the Milky Way, provided kinematics and chemical composition can be reliably measured for very large numbers of them. Metal-deficient, M subdwarfs, in particular, are kinematically associated with the old Galactic halo and thick disk populations and carry the nucleosynthetic signatures of the early evolution of the Galaxy.  M dwarfs/subdwarf are more reliable probes of chemical composition in old populations because, due to the complexity of  the evolutionary sequence of more massive stars,  various enrichment processes may change their chemical pattern. For example, “atomic diffusion” can perturb the atmospheric elemental abundances of of more massive dwarfs with shallow convection zones, i.e.,  the higher the mass, the larger the effect (Gao et al. 2018). In comparison to more massive FGK dwarfs, atomic diffusion has a negligible impact on M-type dwarfs  because of their lower mass and deeper convection zone. The Mixing (“dredge-up”) process can also vary the surface abundance values to some extent in  giants, causing an offset in detailed chemical abundances  between dwarfs and giants  even with a common birthplace. Being unevolved main-sequence stars, there are therefore no ambiguities in the surface composition  of M dwarfs caused by the dredge-up process. Consequently,  the atmospheric abundances of M-type dwarfs  are pristine indicators of  the chemical properties of the progenitor molecular clouds where these stars formed (unless they merge or exchange mass with exoplanets, comets, or their companion in binary systems).

M dwarfs have also been shown to be excellent targets to measure the isotopic abundance ratios of titanium (an $\alpha$ element) as a robust test of Galactic chemical evolution (GCE) models (Chavez & Lambert 2009).  The accurate measurement of Ti  isotopic ratios in stars and the interstellar medium (ISM)  provides important information about the Galaxy's enrichment history, and  in particular, the production of Ti in different types of supernovae (that yield different isotopic ratios). However, because the thermal and macroscopic broadening of atomic lines in stars is comparable to the isotopic splitting, the isotopic effects can essentially  be observed in molecular bands, rather than in atomic lines (although the small  broadening of atomic lines in the ISM makes it possible to measure isotopic effects even in atomic lines in such environments),  and for this reason, isotopic analyses can be  conceived using the spectra of cool stars, which are dominated by molecular lines. Particularly, the prevailing TiO molecular bands in M dwarf spectra allow the measurement of relative abundances for the most stable isotopes of Ti, resulting in valuable clues on the nucleosynthesis of $\alpha$ elements.

However, M-type dwarfs have not been used to their fullest potential in Galactic archaeological studies within the last decade, as compared to numerous  analyses using FGK dwarfs and  giants. Due to their intrinsic faintness, the acquisition of high-resolution, high signal-to-noise spectra of M dwarfs and M subdwarfs demands large telescopes with long exposure times. High-resolution spectroscopy of M dwarfs has therefore been limited to small samples of very nearby stars, mostly from the Galactic disk (e.g., Mann et al. 2013; Rajpurohit et al. 2014, 2018; Veyette et al. 2017; Woolf 2020). On the other hand,  low-resolution spectra of M dwarfs/subdwarfs  have been acquired  in large numbers using more economically efficient observations (e.g., West et al. 2011; Lepine et al. 2013; Zhang 2021). To exploit such a wealth of data, we have attempted to develop an automated pipeline to infer stellar parameters, i.e., effective temperature T$_\textrm{\footnotesize{eff}}$, overall metallicity [M/H], alpha-element enhancement [$\alpha$/Fe], and surface gravity log g for  M dwarfs/subdwarfs by fitting low-resolution spectra to the BT-Settl model grid calculated from the PHOENIX model atmospheres (see Hejazi et al. 2020, hereafter Paper I, for a detailed description). The resulting  values  enable us to analyze the correlation between these parameters, in particular between [M/H] and [$\alpha$/Fe] that has proven powerful in detecting and identifying distinct stellar populations (e.g., Adibekyan, 2012, 2013; Hayden et al. 2015 and references therein). The relationship between the stellar chemical parameters and  photometric and kinematic properties of stars is  also of great importance in distinguishing different  populations (e.g., Lepine et al. 2007;  2013; Hawkins et al. 2015; Yan et al. 2019, also see Section 9 in Paper I for more useful references in this regard).

As it turns out,  different spectroscopic  studies utilize different methods using spectra with different  resolutions and over different wavelength regimes, and  a comparison of the results unfortunately shows significant systematic offsets in  abundance measurements. To assemble a precise  abundance dataset, it is thus critical to obtain  chemical abundances  of all stars in  a sample using the same technique over the same wavelength region with observed spectra of  comparable  resolutions. Although our low-resolution spectroscopic analyses  cannot result in  highly accurate values of   chemical parameters,  since we apply the same technique over a certain spectral range to all stars, the resulting  measurements are precise enough to separate  different stellar population using photometric, chemical and kinematic diagrams.

The work reported here  is an improvement upon the method described in Paper I with a more extensive investigation of the model-fit process and its complications. We revisit the spectra from the sample of Paper I, which focused on high proper-motion stars ($\mu{\geq}0.4{\arcsec}{\textrm{yr}}^{-1}$),  but also expand the analysis to an additional 2045  spectra from nearby (d$<$25 pc) stars with lower proper motions, down to  $\mu=0.03{\arcsec}{\textrm{yr}}^{-1}$. In Section 2, we  describe the spectroscopic sample used in this study. The selected model grid is detailed in section 3. The correlation and degeneracy effect between physical parameters that may cause systematic uncertainties in the  best-fit parameters values are also analyzed in Sections 3. We outline our revised model-fit technique in Section 4 and the error analysis in Section 5.  The  resulting values from our pipeline and their distributions in several parametric and photometric diagrams are presented in Sections 6 and 7. The effect of surface gravity variation and the residual correlation between wavelength datapoints on the best-fit parameter estimates are also addressed in these sections. We analyze the precision of the model-fit chemical parameter values  using a set of common proper-motion pairs in Section 8. The distributions of our stars in the kinematic planes as well as in the abundance-velocity diagrams are examined in Section 9. Finally, we summarize our work in Section 10.

§ SPECTROSCOPIC SAMPLE

Our dataset consists of spectra from 3745 M dwarfs/subdwarfs collected between 2002 and 2012 in multiple observing runs at the MDM observatory, the Kitt Peak National Observatory (KPNO),  the Lick Observatory, and the Cerro Tololo Interamerican Observatory (CTIO). A summary of the observations and data reduction is presented in Paper I (Section 2). For the present analysis, spectral variances were re-calculated from the original spectra using  the read noise and Poisson statistics of the integrated photon count at each spectral pixel, evaluated before sky subtraction and spectrophotometric calibration.

All stars are classified using the template-fit method outlined in Paper I (section 3)  into twelve metallicity classes (from 1 for very metal-rich stars to 12 for very metal-poor stars) and into standard spectral subtypes. The sample consists of stars with spectral subtype between M0.0 and M6.5  and with 2MASS magnitudes K$\leq$15, excluding  objects with unreliable parallaxes and/or photometric magnitude measurements from Gaia DR2/EDR3 (Gaia collaboration et al. 2018; 2021). All the objects used here are either stars that were originally selected for follow-up observations because of their high proper motions ($\mu{\geq}0.4{\arcsec}{\textrm{yr}}^{-1}$) as measured in the SUPERBLINK proper motion survey (Lépine & Shara 2005), or because their pre-Gaia photometric distances suggested they were within 25 parsecs of the Sun (Lépine 2005). All the  brighter objects  with J$<$9  (1564 stars) were previously presented and analyzed in Lépine et al. (2013), and the subset of stars with large proper motions ($\mu{\geq}0.4{\arcsec}{\textrm{yr}}^{-1}$, 1544 stars) was presented and analyzed in Paper I.

Data for the 3745 stars are presented in Tables 1 and 2. Table 1 lists right ascensions, declinations, parallaxes, and proper motions obtained from Gaia EDR3 for 3545 of the stars; entries could not be found in EDR3 for the remaining 200 objects, and for those we are listing their data from the Gaia DR2 catalog. The Gaia flag in the last column of Table 1  shows the  Data Release from which the data is drawn. Table 2 presents the optical and infrared magnitudes of the stars as listed in the Gaia catalogs (DR2 or EDR3 as appropriate) and from  the 2MASS survey.

§ SYNTHETIC SPECTRA

§.§ Grid Selection

We utilize the latest version of BT-Settl synthetic spectra generated by the PHOENIX model atmospheres, as described in Paper I (Section 5.1). Our grid selection is almost the same as explained in Paper I (Section 5.2), limited to T$_\textrm{\footnotesize{eff}}$ from 2700 to 4000 K (consistent with our  selected M0.0-M6.5 spectral type range) in steps of 100 K, log g from 4.5 to 5.5 dex in steps of 0.5 dex, [M/H] from $-$2.5 to +0.5 dex in steps of 0.5 dex, and with a range of [$\alpha$/Fe] values: from $-$0.2 to +0.4 dex for [M/H]$\geqslant$$-$0.5 and from 0 to +0.6 dex
for [M/H]<$-$0.5 dex, in steps of 0.2 dex[Micro-turbulence velocity, $\xi$, is another parameter that can affect spectral absorption lines, but is assumed to be important  in the radiative part of the atmosphere when calculating  BT-Settl synthetic spectra. The turbulent convective velocity as determined by the mixing length theory is added to micro-turbulence velocity in quadrature, but in practice, only increases it by a few percent in M dwarfs, and mainly in deeper layers that do not significantly contribute  to the line core formation. This parameter  is estimated for M dwarfs to a first approximation assuming a simple   function of  effective temperature: $\xi$=(T$_\textrm{\footnotesize{eff}}$ - 1250 K)/7200 (km/s).].

We perform  a four-dimension spline interpolation at every step of 50 K in T$_\textrm{\footnotesize{eff}}$, 0.05 dex in [M/H], 0.025 dex in  [$\alpha$/Fe], and 0.05 in log g, giving  a total of 869,211 grid points in the 4D parameter space. In the present work, we introduce an expanded  set of grid points in surface gravity, i.e., log g from 4.5 to 5.5 dex in steps of  0.05 dex, whereas Paper I examined a more limited set of gravity grid points  from 4.8 to 5.2 dex with a step size of 0.1 dex. Moreover, to reduce the computational time of   model-fit runs,  we first apply a Gaussian broadening kernel and then resample all original and interpolated synthetic spectra across  a wavelength step of 0.1 Å. Our analysis shows that the resulting model-fit parameters of our observed spectra (which have a spectral resolution of $\simeq$2-3 Å) do not change when using these down-sampled synthetic spectra[We randomly selected 300 stars from our sample and applied the pipeline to this subset using both  original and down-sampled synthetic spectra and found no significant difference in the two sets of results.], as compared to parameters inferred from model spectra sampled on the original, high-resolution wavelength grid (0.02 Å).

§.§ Spectral Structure and Synthetic Fitting

In the presence of low temperatures and high pressures,  molecules such as  TiO, VO, CO, CaH, FeH, MgH, H$_\textrm{\footnotesize{2}}$, and H$_\textrm{\footnotesize{2}}$O are the main opacity source in  the atmospheres of M dwarfs (Allard & Hauschildt 1995). The modeling of stellar atmospheres therefore requires  accurate opacity data  for millions of molecular transitions and spectral lines. These molecules have significant impacts on  the detailed structure of M dwarf spectra, influencing  both individual lines and  the overall spectral shape through the cumulative effect of large numbers of overlapping lines at various wavelength ranges (Valenti et al. 1998).

The forests of molecular lines hinder the identification of  continuum levels even at high spectral resolution, which  complicates the accurate determination of atomic line strengths and also introduce obstacles in the techniques that are readily applied to spectra of more massive dwarfs, i.e., F, G, and K  dwarfs, for deriving their physical parameters. The most accurate approach to determine the stellar parameters  of M dwarfs is to normalize both the observed spectrum  of interest and  all the synthetic spectra relative to their corresponding continuum, and compare the resulting rectified spectra to find the best-fit model. Although the continuum of synthetic spectra  of M dwarfs can, in principle, be calculated using the number density of key species contributing to the continuous opacity such as H, H$^\textrm{\footnotesize{$-$}}$, and H$_\textrm{\footnotesize{2}}$, the dominant molecular bands make it difficult to properly define the continuum level of an observed M dwarf spectrum, particularly in the optical regime where the TiO molecular lines cause substantial flux depressions against the continuum background. This is a significant issue because the comparison between a continuum-normalized observed spectrum and a set of continuum-normalized synthetic spectra,   with regard to the real depth of spectral lines and features (measured from the continuum) that  indeed depend on  atmospheric parameters, is expected to be the most reliable way to search for the best-fit model. Unfortunately, this approach is not applicable  to M dwarfs, which poses challenges in estimating their stellar parameters.

There are two common methods used to perform  synthetic fitting over broad spectral ranges when the continuum cannot be correctly determined.  In one method, both the observed spectrum and  the synthetic spectra are normalized to a pseudo-continuum profile separately, and these pseudo-rectified spectra are then used  in the minimization procedure to obtain the best-fit model (e.g., Lee et al. 2011a; Kuznetsov et al. 2019).  In another approach, each synthetic spectrum in the  model grid, or alternatively the observed spectrum, is flux-renormalized by a polynomial that, in either case, can be determined from the ratio of the observed spectrum to each of the synthetic spectra (e.g., Lépine et al. 2013; Hejazi et al. 2020), or  equivalently, from the loss function as described in  Zhang et al. (2021). The flux-renormalized  synthetic spectra (or observed spectrum) will then be compared to the unchanged observed spectrum (or  synthetic spectra)  to find the best matched model. In both methods, the observed spectrum is compared to the synthetic models in terms of  “apparent line depths”,  that are measured from  pseudo-continuum levels,  and  “apparent spectral features”, effectively overlooking potentially important physical information embedded in the real continuum.

§.§ Spectral Sensitivities

Prior to describing our model-fit technique, we  aim to investigate the sensitivity of synthetic spectra to their physical parameters and examine  how their spectral morphology changes when each parameter is varied. This  can help one understand how two synthetic  spectra with different parameters may be  degenerate, which can lead to possible systematic errors in resulting model-fit parameters.  To this end, we illustrate  the change in the  apparent spectral shape  of  synthetic spectra that are degraded to the typical resolution of our observed spectra.  For ease of comparison, we normalize these spectra relative to the flux at an arbitrary wavelength, e.g., the longest wavelength in the plots mentioned below (8725 Å), setting the flux at this wavelength to unity.

The shape, depth, and broadness of prominent spectral lines are found to vary  with physical parameters in a complicated way. Figures 1-4 show how the spectral shape changes as one physical parameter (effective temperature, metallicity, $\alpha$-element enhancement, or surface gravity) is varied while the other parameters are kept fixed.  The position of some key elements and the regions over which certain molecular bands have an important effect are also presented in these figures. In Figure 1, synthetic spectra are shown for different values of temperature,  i.e., T$_\textrm{\footnotesize{eff}}$=3000, 3200, 3400, 3600, 3800, and 4000 K; all models have  [M/H]=0.0 dex,  [$\alpha/$Fe]=+0.2 dex, and  log g=5.0 dex. As T$_\textrm{\footnotesize{eff}}$ decreases,  the star becomes significantly redder and the overall flux drops. In cooler stars,  there is less  thermal energy to populate higher atomic states, so only low-excitation atomic lines remain prominent. Furthermore, with decreasing effective temperature,  more atoms contribute to  form molecules and grains, which increases the effect of molecular bands on individual atomic spectral lines. While some atomic lines such as Ca I ($\sim$6104, 6124, 6164, 6441, 6463, and 6496 Å) Ca II ($\sim$8500, 8544, and 8664 Å), Li I ($\sim$6710 Å), Fe I ($\sim$8329 and 8390 Å), Ti I ($\sim$8438 Å), and V I ($\sim$8692 Å) are weaker at low temperatures, the alkali atomic lines, i.e.,  the doublet  K I ($\sim$7667 and 7701 Å) and Na I ($\sim$8186 and  8197 Å), are  prominent even in the spectra of very cool M dwarfs. These alkali atomic lines are therefore of importance in the determination of atmospheric parameters, especially at low spectral resolution.

The dependence of spectral morphology on metallicity is more complex. Figure 2 depicts the model spectra with [M/H]=$-$2.0, $-$1.5, $-$1.0, $-$0.5, 0.0, and +0.5 dex; all these models are associated with T$_\textrm{\footnotesize{eff}}$=3400 K, [$\alpha/$Fe]=+0.2 dex, log g=5.0 dex. As metallicity increases, the molecular bands become progressively stronger, increasingly depressing the flux level.  It is important to note that with increasing metallicity, the TiO bands (and more likely other oxide metals) grow stronger, but the CaH (and also other hydrides) bands are not affected. This is because  when metallicity increases, the gas pressure decreases and thus the concentration of hydride molecules remains nearly the same (Bessell 1991). The effect of metallicity on  real and apparent line strengths differs from one line to another, which is  due to  the correlation of various factors with stellar metal content. In general, the spectral shape dramatically changes at low metallicities;  the  weaker the molecular bands are, the more prominent individual atomic lines become. Further, due to different levels of flux depression, the spectral features in metal-poor spectra significantly differ from those of metal-rich ones, resulting in significant  changes in broadband colors even for stars with similar effective temperatures. This is the reason why our metal-poor stars ([M/H]$\leq$$-$1.0 dex) are normally well separated from the metal-rich ones ([M/H]$\geq$$-$0.5 dex)  using the inferred metallicities from the synthetic fitting (see section 6).

Spectral sensitivity to effective temperature: synthetic spectra with [M/H]=0.0 dex,  [$\alpha/$Fe]=+0.2 dex, and  log g=5.0 dex, but different  values of T$_\textrm{\footnotesize{eff}}$= 3000, 3200, 3400, 3600, 3800, and 4000 K. The spectra are normalized at 8725 Å.

Spectral sensitivity to metallicity: synthetic spectra with T$_\textrm{\footnotesize{eff}}$=3400 K, [$\alpha/$Fe]=+0.2 dex, log g=5.0 dex, and [M/H]=$-$2.0, $-$1.5, $-$1.0, $-$0.5, 0.0, and +0.5 dex. The spectra are normalized at 8725 Å.

Spectral sensitivity to $\alpha$-element enhancement: synthetic spectra with T$_\textrm{\footnotesize{eff}}$=3400 K, [M/H]=0.0 dex,  log g=5.0 dex, and [$\alpha/$Fe]=0.0, +0.1, +0.2, +0.3, and +0.4 dex. The spectra are normalized at 8725 Å.

Spectral sensitivity to surface gravity: synthetic spectra with T$_\textrm{\footnotesize{eff}}$=3400 K, [M/H]=0.0 dex, [$\alpha/$Fe]=+0.2 dex, and log g=4.5, 4.8, 5.0, 5.2, and 5.5 dex. The spectra are normalized at 8725 Å.

Spectral degeneracy from effective temperature and metallicity: synthetic spectra with [$\alpha/$Fe]=0.2 dex, log g=5.0, and  T$_\textrm{\footnotesize{eff}}$=3000, 3100, 3200, and 3300 K and [M/H]=$-$0.5, $-$0.25, 0.0, and +0.25 dex, respectively. The spectra are normalized at 8725 Å.

Spectral degeneracy from effective temperature and $\alpha$-element enhancement: synthetic spectra with [M/H]=0.0 dex, log g=5.0,  T$_\textrm{\footnotesize{eff}}$=3100, 3200, and 3300 K, and [$\alpha$/Fe]= 0.0, +0.2, and +0.4 dex, respectively. The spectra are normalized at 8725 Å.

Spectral degeneracy from effective temperature and surface gravity: synthetic spectra with [M/H]=$-$0.5 dex, [$\alpha$/Fe]=+0.2 dex,    T$_\textrm{\footnotesize{eff}}$=3000, 3100, 3200,  and 3300 K, and log g=5.1, 4.9, 4.7, and 4.5 dex, respectively. The spectra are normalized at 8725 Å.

Spectral degeneracy from metallicity and $\alpha$-element enhancement: synthetic spectra with T$_\textrm{\footnotesize{eff}}$=3400 K and log g=5.0 dex, but five different combinations of [M/H] and [$\alpha$/Fe]: ([M/H]=$-$0.2 dex, [$\alpha$/Fe]=+0.0 dex), ([M/H]=$-$0.5 dex, [$\alpha$/Fe]=+0.0 dex), ([M/H]=$-$0.5 dex, [$\alpha$/Fe]=+0.1 dex), ([M/H]=$-$0.5 dex, [$\alpha$/Fe]=+0.2 dex), and ([M/H]=$-$0.5 dex, [$\alpha$/Fe]=+0.3 dex). The spectra are normalized at 8725 Å.

Spectral degeneracy from metallicity and surface gravity: synthetic spectra with  T$_\textrm{\footnotesize{eff}}$=3400 K,  [$\alpha$/Fe]=0.0 dex, [M/H]=+0.2, +0.2, and +0.4 dex, and log g=5.0, 5.5, and 5.5 dex, respectively. The spectra are normalized at 8725 Å.

Spectral degeneracy from $\alpha$-element enhancement and surface gravity: synthetic spectra with  T$_\textrm{\footnotesize{eff}}$=3400 K,  [MH]=0.0 dex,  [$\alpha$/Fe]=+0.0, +0.0, and +0.1 dex, and log g=5.0, 5.4, and 5.4 dex, respectively. The spectra are normalized at 8725 Å.

The influence of $\alpha$-element enhancement on spectral shape is shown in Figure 3, which displays the spectra with [$\alpha/$Fe]=0.0, +0.1, +0.2, +0.3, and +0.4 dex but all having the same T$_\textrm{\footnotesize{eff}}$=3400 K, [M/H]=0.0 dex, log g=5.0 dex.  As can be seen from the figure, higher values of [$\alpha/$Fe] generally cause a higher degree of flux depression from molecular bands. The atomic lines are remarkably influenced by the background molecular opacities, but in different ways. For example, the $\alpha$-element lines, i.e., Ca I, Ca II,  and Ti, are significantly blended with the ambient molecular bands at higher [$\alpha/$Fe] values. On the other hand, there is an opposite trend for the other atomic lines such as the doublet Na I that grow stronger when [$\alpha/$Fe] decreases. Decreasing the abundances of  $\alpha$ elements may lead to  increasing the abundances of some other elements  to keep  the overall metallicity constant for the plotted model spectra, though many factors have to be taken into account to completely explain the effect of  [$\alpha/$Fe] on spectral lines and features.

Figure 4 demonstrates the change of spectral structure due to the variation of surface gravity, showing the spectra with different values of log g=4.5, 4.8, 5.0, 5.2, and 5.5 dex but the same T$_\textrm{\footnotesize{eff}}$=3400 K, [M/H]=0.0 dex, and [$\alpha/$Fe]=+0.2 dex. Clearly, surface gravity has a weaker effect on overall flux and spectral shape, even over a wide  range of 1.0 dex, as compared to the other parameters.  Nevertheless, this parameter plays an important role in the determination of the other stellar parameters using  synthetic model fitting (Section 6). The sensitivity of spectral lines to  log g varies from one element to another.  For instance, the doublet K I and Na I become stronger as log g increases, which can be explained by the variation of  the population of atomic states with pressure.  Based on numerical tests, we find that for a typical  model with a low value of  log g, the majority of  Na  and K atoms are  ionized while Ca and Ti are mainly  neutral.  Increasing log g thus increases the population of neutral Na I and K I atoms, but has no effect on the population of neutral Ca I and Ti I atoms.

All these examinations imply that each parameter has a unique effect on the flux level and the full structure of synthesized spectra over a wide spectral range. Some combinations of the parameters
shown in Figures 1-4 may not, however, happen in real stars. For example, since  M dwarfs don't age,  a change in log g for a fixed mass (or T$_\textrm{\footnotesize{eff}}$) must be due to changes in the abundance scale. In other words, stars with the same T$_\textrm{\footnotesize{eff}}$ and the same chemical parameters will have nearly the same log g. This means while other parameters are constant, the variation of log g will be extremely narrow in reality, as compared to the broad range  (1.0 dex)  presented  in Figure 4.

In effect,  the spectral sensitivity to a particular physical parameter  depends  on the value of the other parameters.  The shape of spectral lines and features changes in a complex  way when more than one parameter vary, and for this reason,  to fully understand the exact dependence of spectral structure on stellar parameters, extensive spectroscopic analyses  combined  with high-resolution optical and infrared spectra (which is beyond the scope of the present investigation) may ultimately be required.

§.§ Spectral Degeneracies and Systematics

Apparent similarities and degeneracies between spectra  have been shown to be a serious problem in  M dwarf synthetic fitting   (e.g., Passegger et al. 2016 and 2018; Rains et al. 2021). Despite the complicated dependence of spectral shapes on physical parameters, there are
clear similarities in the overall morphology of model spectra generated using different
parameter values, which are expected to cause  significant  levels of  degeneracy in model fitting,
in particular when the observed/synthetic spectrum is flux-renormalized to deal with the absence of a clearly defined continuum, as is usually the case. The perplexing correlations between parameters makes a full degeneracy analysis difficult, and in this section, we will  illustrate the cases arising from the simultaneous variation of two parameters.  A complete exploration  of  spectral degeneracy when more than two parameters are involved remains very challenging and is beyond the scope of the present paper. Our model-fit results indicate systematic relationships  between physical parameters  that likely have physical origins (for example, a reduced metallicity is expected to naturally yield a higher  $\alpha$-element enhancement), other correlations are likely due to issues in model atmospheres or  systematic errors in inferred parameter values resulting from degeneracies in synthetic spectra, or even a combination of all these factors. Our purpose here is to attempt  to find  the underlying reason for these  trends, however, one should keep in mind that the broad conclusions below  are only valid for modeling of low-resolution spectra, because at high resolutions, the more detailed spectral structures can break the degeneracies that are plaguing  spectral modeling at low resolutions.

Effective Temperature and Metallicity: The correlation between effective temperature and metallicity is essential in  synthetic fitting and resulting parameter values.  As shown in Section 6, our observed spectra resist being fit with high-temperature (T$_\textrm{\footnotesize{eff}}$$\gtrsim$3600 K), high-metallicity  ([M/H]$\gtrsim$0.2 dex) model spectra, no matter which initial values, minimization formalism, or  iterative process are chosen in the model-fit pipeline.  And this is despite the fact that many of our lower temperature stars are fitted with high-metallicity model spectra, which means that at least some of our hotter stars should be high-metallicity objects as well. This indicates that  high-temperature, high-metallicity synthetic models may not be  in agreement with observations, which makes a perceptible systematic trend  in our estimated parameters. In the atmosphere of M dwarfs, the number of molecules decreases with increasing  effective temperature, which in turn reduces the influence of molecular bands on the spectral flux. On the other hand, as metallicity increases, the abundance of heavy elements rises, which makes the effect of molecular lines stronger. Apparently, model-fits of low-resolution spectra are not able to correctly disentangle these two effects for  high-temperature, high metallicity M dwarfs, and most likely tend to fit high-metallicity stars with models of slightly lower effective temperatures as a response to seeing stronger molecular bands. This problem is compounded by the fact that high-temperature M dwarfs have weak molecular bands (more like K dwarfs) which provides reduced leverage in the model-fitting compared with cooler M dwarfs whose molecular bands are very strong.

An examination of synthetic spectra also shows some sources of  degeneracies that arise from the coincident variation of T$_\textrm{\footnotesize{eff}}$ and [M/H]. As [M/H] rises,  molecular band opacities increase, which depress the continuum and, in turn, decrease apparent line depths (as measured from the local pseudo-continuum).  These changes can be somewhat nullified by increasing T$_\textrm{\footnotesize{eff}}$, which has the opposing effect of reducing molecular opacities and increasing apparent line depths. This effect is illustrated in Figure 5. There is, for example, a clear resemblance between any two models with a difference of 100 K in T$_\textrm{\footnotesize{eff}}$ and  0.25 dex in [M/H] over many wavelength regions. However, the significant discrepancies  in some line strengths and features  between models with more than 100 K apart in T$_\textrm{\footnotesize{eff}}$ break this degeneracy.  We find a variation of T$_\textrm{\footnotesize{eff}}$ within $\sim$100 K with a simultaneous  variation of [M/H] within $\sim$0.25-0.30 dex may lead to degeneracy,  at both  high and low metallicities. Such degeneracies cannot, however, explain the above-described systematic trend that only occurs for the high-temperature, high-metallicity  stars.

Effective Temperature and $\alpha$-element Enhancement:  As detailed in Section 6, our model-fit  results  also show a systematic trend between T$_\textrm{\footnotesize{eff}}$ and  [$\alpha$/Fe] mainly for  the near-solar-metallicity stars ($-$0.5$\leq$[M/H]$\leq$+0.5 dex), again regardless the initial values, the  $\chi$$^\textrm{\footnotesize{2}}$ formalism, and the extent of the examination model grid used in the  minimization processes.  Our high-temperature stars tend to fit with low values of [$\alpha$/Fe], while  low-temperature stars are better fitted with high values  of [$\alpha$/Fe], which is very unlikely to represent a real correlation since our stars are selected from a volume-limited sample, i.e., they should all have similar distributions of [$\alpha$/Fe] regardless of temperature. Instead, this effect likely points to issues with the synthetic spectra. Most of the optical region in M dwarf spectra is strongly influenced by the  molecular lines of  TiO that is composed of two $\alpha$ elements (i.e., Ti and O). As effective temperature increases, the number of TiO molecules decreases, which should not be confused with a decrease in $\alpha$-element abundances. However, as one can see from Figure 3, the decrease in [$\alpha$/Fe] also reduces the effect of molecular bands on the overall flux, and this can be the reason why higher-temperature stars are better matched with models of lower  [$\alpha$/Fe] values. This trend therefore introduces systematic errors in modeling the atmosphere of M dwarfs which cannot properly differentiate  between the effect of the TiO molecular bands and $\alpha$-element abundances.

The two parameters T$_\textrm{\footnotesize{eff}}$ and [$\alpha$/Fe] can also cause spectral degeneracies to some extent. An increase in [$\alpha$/Fe] does not affect the resulting spectrum in the exactly  opposite way as an increase in T$_\textrm{\footnotesize{eff}}$. Despite this fact, there are still similarities in low-resolution spectra over some important wavelength regions between any two  spectra with a difference of 100 K in T$_\textrm{\footnotesize{eff}}$ and 0.2 dex in  [$\alpha$/Fe],   as shown in Figure 6. However, the significant differences in spectral morphology between models with effective temperatures that differ by more than 100 K restrict the degeneracy to a relatively narrow range of T$_\textrm{\footnotesize{eff}}$ values. All in all, a variation of T$_\textrm{\footnotesize{eff}}$ within $\sim$100 K along with a variation  of  [$\alpha$/Fe] within $\sim$0.2 dex   can produce relatively degenerate synthetic spectra, in both high- and low-metallicity regimes, but such a temperature change  cannot be responsible for the above-described systematic trend observed in our resulting model-fit parameters.

Effective Temperature and Surface Gravity: The relation between effective temperature and surface gravity is of importance in characterizing stellar atmospheres. As presented in Section 6, our synthetic fitting shows a tendency for the higher-temperature stars to be better matched with lower values of log g, and inversely, the lower-temperature stars are better matched with higher values of log g, primarily for stars of  near-solar metallicities, which is  consistent with the prediction of stellar  evolutionary theory (Section 4).  In addition, these two  parameters can induce degeneracies when varied together simultaneously, as presented in Figure 7. For instance, any two spectra with a difference of 100 K in T$_\textrm{\footnotesize{eff}}$ and 0.2 dex in log g are quite similar, giving rise to   degeneracy effects in the synthetic model fitting. However, the considerable discrepancies in some atomic lines and features between the models with more than 100 K apart  in T$_\textrm{\footnotesize{eff}}$ do not result in systematic errors in the inferred model-fit parameters.   We note a variation of T$_\textrm{\footnotesize{eff}}$  within $\sim$100 K together with  a variation of log g within $\sim$0.2 dex  may yield degenerate spectra,  at both high and low metallicities,   although such degeneracies may not be extensive enough to cause the  observed relation between these two parameters in our results.

Metallicity and $\alpha$-element Enhancement:  There is a general correlation between the values of metallicity and $\alpha$-element enhancement derived from our model-fit pipeline: [$\alpha$/Fe] increases as  one moves from metal-rich M dwarfs to metal-poor M subdwarfs, which is  in line with the trends predicted by  Galactic chemical evolutionary  models (Section 6). But because $\alpha$ elements have a dominant effect on M dwarf spectra, these two chemical parameters affect the spectral shape in roughly the same way and an increase in [M/H] may be counteracted by a decrease in [$\alpha$/Fe] and vice versa. However, the change of spectral shape due to  a large variation in [M/H] (e.g., by  $\gtrsim$ 0.5 dex) cannot be fully compensated by adjusting [$\alpha$/Fe]  over our selected grid range. As a result, the degeneracy can only happen when these two parameters vary over small ranges, as shown in Figure 8. Considering the model with  [M/H]=$-$0.2  and [$\alpha$/Fe]=0.0 dex (red) as the reference, we begin with the model having  [M/H]=$-$0.5  and [$\alpha$/Fe]=0.0 dex (blue) and increase  [$\alpha$/Fe] by 0.1 (black) and then by 0.2 dex (green). One can see how the spectral shape becomes more similar   to the reference model by varying [$\alpha$/Fe] over many spectral regions. However, increasing [$\alpha$/Fe] by $\gtrsim$ 0.3 dex (yellow)  makes the spectral shape more  dissimilar to  the reference model, and this is where the degeneracy breaks. On closer examination, we find that a variation of [M/H] within $\sim$0.3-0.35 dex and a variation of  [$\alpha$/Fe] within  $\sim$0.2 dex at the same time may create degenerate synthetic spectra in both low- and high-metallicity regimes. While such degeneracies are too limited to be responsible to the  global  trends in the [M/H]-[$\alpha$/Fe] diagram, some local trends associated with other parameters may be attributed to these spectral degeneracies (Section 6). Apart from this degeneracy, these two parameters are tightly correlated together over larger ranges, as shown in Section 5.

Metallicity and Surface Gravity:  The correlation between metallicity and surface gravity is critical  in  synthetic-model fitting. Our model-fit results show a systematic relation between these two parameters when surface gravity is allowed to vary, in spite of the selected initial values, the minimization procedure, or the extent of the examination grid, as described more at length in Section 6.  There appears to be some degeneracies between spectra with a small difference in metallicity  but with a  relatively large difference in surface gravity (as compared to the parameter ranges in our model gird), as illustrated  in Figure 9.  A  decrease in metallicity by 0.2 dex  normally raises the flux level and increases the apparent line depths, which may be counteracted with a decrease in  surface gravity by 0.5 dex to some extent over a significant part of the selected wavelength coverage. However, the spectral morphology begins to become  inconsistent as the difference in  metallicity and surface gravity increases, which breaks the degeneracy. We find a metallicity change  within $\sim$0.2 dex and a gravity change  within $\sim$0.5 dex can create degenerate synthetic spectra, resulting in systematic trends in the best-fit parameter values, mostly for stars of  near-solar metallicities. This can contribute to the observed relationship between [M/H] and log g in our results for  [M/H]$\geq$$-$0.5 dex, when log g is treated as a free parameter; i.e., the higher-metallicity stars tend to have higher values of surface gravity and the lower-metallicity stars tend to have lower values of surface gravity. However, the effect is not as important for metal-poor spectra, which implies that degeneracies from these two parameters are less significant  in the low-metallicity regime.  This is the reason why  our results show no considerable systematic relation between metallicity and surface gravity for metal-poor M subdwarfs.

$\alpha$-element Enhancement and Surface Gravity:
As will be shown in section 6 below, there is   a  trend between our [$\alpha$/Fe] best-fit estimates and the inferred values of log g from empirical photometric relations (which are kept fixed in the model fitting), particularly for  high-metallicity stars.  This trend follows  in very much the same way as the trend between [$\alpha$/Fe]  and T$_\textrm{\footnotesize{eff}}$, which indeed reflects the tight correlation between T$_\textrm{\footnotesize{eff}}$ and log g. On the other hand,  there is a weak degeneracy effect between [$\alpha$/Fe] and log g when both are varied together. Figure 10 compares  the  spectral shapes of three models with distinct values of  log g and [$\alpha$/Fe].  There is some level of similarity in the spectral morphology of two   spectra with a difference of 0.1 dex in [$\alpha$/Fe] and 0.4 dex in log g. The effect is, however, not as strong as the degeneracy between [M/H] and log g, but can be  responsible for the systematic shifting of some metal-rich stars to higher values of   [$\alpha$/Fe]  and log g (Section 6).  We also find that this degeneracy begins to  break once the difference in  [$\alpha$/Fe] and  log g increases. Our examination shows a weaker degeneracy effect for metal-poor spectra, and as a consequence, there is no noticeable trend between the model-fit values of these two parameters in the low-metallicity regime.

As mentioned above, a thorough investigation of the correlation between parameters needs more careful analyses and examinations.  Detailed spectral shapes are influenced by all stellar parameters together, making the understanding of degeneracy effects difficult. Overall,  when  observed spectra provide limited information  because of low spectral resolutions, instrumental noise, and renormalization of the continuum, these spectra can be easily fitted to a range of synthetic spectra with different stellar parameter values, which introduces significant uncertainties in the spectral modeling. Ideally, one would want to know the effective ranges over which the different physical parameters may produce degenerate synthetic spectra, and thus determine the level of uncertainty for each physical parameter for any specific star. This is a challenging proposition in general, but the rough estimates provided above in a few specific cases provide a general guideline to what one should expect.

The executive summary is  effective temperature has the largest effect on the shape and flux level of spectra, and for this reason, it cannot make significant degeneracies along with other parameters beyond a variation of $\sim$100 K, and this is the reason why this parameter is usually well determined by our spectroscopic model fitting. Any trend or relationship  related to T$_\textrm{\footnotesize{eff}}$  is most likely due to a physical origin or deficiencies in the model atmospheres.  On the other hand, the three other parameters, i.e.,  [M/H], [$\alpha$/Fe], and  log g, produce variations in the synthetic spectra that are strongly correlated, and this causes significant degeneracies at low spectral resolutions. Most notably, there may be a critical degeneracy effect in the model fitting (Section 6.3) when the change of spectral morphology due to a relatively large variation of a parameter (e.g.,  log g) can be compensated  by
a relatively small variation of another parameter (e.g., [M/H] or [$\alpha$/Fe]).

Photometric surface gravity versus photometric effective temperature for the 3745 stars (blue dots)  along with the isochrones for a  stellar age of 10 Gyr, and  chemical parameters [M/H]=+0.0 dex and [$\alpha$/Fe]=+0.0, +0.2, and +0.4 dex (solid, dashed and dotted black lines, respectively),  [M/H]=$-$1.0 dex and [$\alpha$/Fe]=+0.2, +0.4, and +0.6 dex (solid, dashed and dotted red lines, respectively), and [M/H]=$-$2.0 dex and [$\alpha$/Fe]=+0.2, +0.4, and +0.6 dex (solid, dashed and dotted green lines, respectively) from the Dartmouth stellar evolutionary models.

Photometric surface gravity versus photometric effective temperature for the 3745 stars (the same as shown in Figure 11), color-mapped according to their metallicity classes.

§ SYNTHETIC-FIT METHOD

In comparison to our previous work, we extent the fitting spectral region to 6090-8700 Å, which is covered by the majority of our observed spectra. However, there is a  number of  stars whose spectra do not fully encompass this wavelength region and may have a starting wavelength  longer than 6090 Å or an ending wavelength  shorter than 8700 Å.  Nevertheless, these stars all minimally cover the spectral range between 6300 Å  and 8600 Å, which includes the most important features sensitive to stellar parameters  (see Paper I). The effect of missing spectral portions  on the resulting best-fit stellar parameters is therefore negligible. We exclude the  region 6550-6580  Å, consisting of the H$_\textrm{\footnotesize{$\alpha$}}$ emission line of hydrogen that is  detected  in some of  our stars, but cannot be reproduced by model atmospheres. The atomic line Li I ($\sim$ 6710 Å) is not well-modeled for a number of synthetic spectra, particularly in the  metal-poor regime, and consequently, we remove the region 6690-6726 Å from the fitting region. In addition,  the wavelength ranges contaminated by telluric absorption bands, i.e., 6855-6885 Å, 7590-7655 Å, 8125-8170 Å, and 8216-8270 Å are also masked. It should be noted that these regions are consistent with the most prominent parts in the transmission spectrum (R$\sim$4000) of the Earth's atmosphere over the KPNO  (Hinkle, Wallace &  Livingston 2003). We also develop a routine to automatically reject highly noisy and problematic regions as well as occasional spectral artifacts from each observed spectrum.

§.§ Minimization Process

We carry out  a number of $\chi$$^\textrm{\footnotesize{2}}$ minimizations, each of which compares one observed spectrum with a number of synthetic spectra associated with  known values of the parameters T$_\textrm{\footnotesize{eff}}$, [M/H], [$\alpha$/Fe], and  log g; these four parameters are referred to as the "primary" parameters. In each minimization, one or more of the primary parameters are allowed to vary (“free”), while the others are kept fixed (“constant”). Radial velocity shift, $v_{r}$, and  convolution factor, C[Synthetic spectra  are convolved by a Gaussian profile to match the spectral resolution of observed spectra.], both referred to as "secondary" parameters, may also be varied or held fixed.  All minimizations follow the same sequence, as summarized below:

* Free parameters  and their corresponding ranges  are determined, and a  specified value is assigned to each constant parameter.

* Synthetic spectra are convolved based on the convolution factor C.

* Observed wavelengths are shifted according to the radial velocity shift $v_{r}$ and the synthetic spectra are interpolated at the shifted wavelengths.

* The observed spectrum is divided by each interpolated synthetic spectrum and the ratio is fitted by a polynomial of  high order  (order of 8) with three successive rejection runs of 10$\sigma_{r}$, 8$\sigma_{r}$, and 6$\sigma_{r}$, where $\sigma_{r}$  is the root mean square error (RMSE) between the quotient and the fit.

* The  interpolated synthetic spectra are multiplied by their corresponding polynomials evaluated at each shifted observed wavelength to renormalize synthetic spectra relative to the observed spectrum.

* The flux-renormalized, interpolated model spectra at the shifted wavelengths, F$_{\textrm{mod}}$($\lambda_{\textrm{i,shift}}$), are compared to the observed spectrum at the unshifted wavelengths, F$_{\textrm{obs}}$($\lambda_{\textrm{i,unshift}}$), by a $\chi$$^\textrm{\footnotesize{2}}$  as defined by

\begin{equation}
\chi^\textrm{\footnotesize{2}} =   {
\sum \limits_{i} \frac{(\mathrm{F}_{\mathrm{obs}}(\lambda_{\mathrm{i,unshift}})  - \mathrm{F}_{\mathrm{mod}}(\lambda_{\mathrm{i,shift}}))^{2}} {\mathrm{E}_{\mathrm{obs}}(\lambda_{\mathrm{i,unshift}})^{2}}}
\end{equation}

where E$_{\textrm{obs}}$($\lambda_{\textrm{i,unshift}}$) is the empirical uncertainty (statistical and instrumental error) of the observed flux at the unshifted wavelengths.

* The best-fit model spectrum minimizes the $\chi$$^\textrm{\footnotesize{2}}$ value.

Our technique comprises several sequential minimizations that  follow the above-described routine, but with different sets of free and constant  parameters and different ranges within which the selected free parameters vary. Resulting parameters from each minimization are used in the next minimization, progressively improving  fit solution.

§.§ Initial Parameter Values

Prior to running the pipeline, we provide an  initial estimate for each primary parameter [M/H], T$_\textrm{\footnotesize{eff}}$, and   log g based on the star's estimated  metallicity class (MC, 1<MC<12 from standard spectral classification) and  available photometric calibrations for M dwarfs. We also provide an initial value for [$\alpha$/Fe] by applying a simple relation between [M/H] and [$\alpha$/Fe] as follows:  we obtain  the  initial  metallicity estimates  based on  the measured metallicity class of our  stars  and the relationship between metallicity class and [M/H] derived in Paper I (Figure 34). Motivated by our previous results and other studies of nearby FGK dwarfs (e.g., Adibekyan et al. 2012, 2013; Recio-Blanco et al. 2014), we define a first-order relationship between [M/H] and [$\alpha$/Fe] to estimate the initial values of [$\alpha$/Fe]:

\begin{equation}
\begin{aligned}
& \mathrm{[\alpha/Fe]}=+0.5 \text{ for }   \mathrm{-2.5 \leq [M/H] \leq -2.0}\\
& \mathrm{[\alpha/Fe]}=+0.4 \text{ for }   \mathrm{-1.9 \leq [M/H] \leq -1.5} \\
& \mathrm{[\alpha/Fe]}=+0.3 \text{ for }   \mathrm{-1.4 \leq [M/H] \leq -1.0} \\
& \mathrm{[\alpha/Fe]}=+0.2 \text{ for }   \mathrm{-0.9 \leq [M/H] \leq -0.5} \\
& \mathrm{[\alpha/Fe]}=+0.1 \text{ for }   \mathrm{-0.4 \leq [M/H] \leq -0.1} \\
& \mathrm{[\alpha/Fe]}=+0.0 \text{ for }   \mathrm{+0.0 \leq [M/H] \leq +0.5}.
\end{aligned}
\end{equation}

We calculate the initial values of T$_\textrm{\footnotesize{eff}}$  using the photometric calibrations  of Mann et al. (2015) who developed a number of metallicity-independent relationships between  T$_\textrm{\footnotesize{eff}}$ and colors associated with the SDSS (r,z), Johnson-Cousins (V,I$_\textrm{\footnotesize{c}}$), 2MASS (J,H), and Gaia (G$_\textrm{\footnotesize{BP}}$,G$_\textrm{\footnotesize{RP}}$)  bandpasses. We choose the relation with the least scatter ($\sigma$=49 K), which is based on the 2MASS and Gaia colors J$-$H and G$_\textrm{\footnotesize{BP}}$$-$G$_\textrm{\footnotesize{RP}}$:

\begin{equation}
\mathrm{T_{\footnotesize{eff}}} = 3500{\{}\sum \limits_{i=0}^{{i=4}}  a_{i}(\mathit{J-H})^{i} + \sum \limits_{i=5}^{{i=6}}  a_{i}(\mathit{G_{\footnotesize{BP}}}- \mathit{G_{\footnotesize{RP}}})^{i}{\}}
\end{equation}

where ${a}_{0}$=3.1720, ${a}_{1}$=$-$2.4750, ${a}_{2}$=1.0820, ${a}_{3}$=$-$0.2231, ${a}_{4}$=0.0174, ${a}_{5}$=0.0878, and ${a}_{6}$=$-$0.0436. We also obtain an initial value for the  radius ($\sigma$=2.89%) and the mass  ($\sigma$=2.0%)  of the stars using the  metallicity-independent, photometric relations from Mann et al. (2015) and Mann et al. (2019), respectively, as follows:

\begin{equation}
log(\frac{\mathrm{R_{*}}}{\mathrm{R_{\sun}}}) = b_{1}+b_{2}\mathit{M_{K}}+b_{3}\mathit{M_{K}}^2
\end{equation}

\begin{equation}
log(\frac{\mathrm{M_{*}}}{\mathrm{M_{\sun}}}) =
\sum \limits_{i=0}^{i=5}  c_{i}(\mathit{M_{K}} - zp)^{i}
\end{equation}

where M$_{\footnotesize{K}}$ is the absolute magnitude in 2MASS K bandpass,  ${b}_{1}$=1.9515, ${b}_{2}$=$-$0.3520, ${b}_{3}$=0.0168, ${c}_{0}$=$-$0.6420, ${c}_{1}$=$-$0.2080, ${c}_{2}$=$-$8.43$\times$10$^{-4}$, ${c}_{3}$=7.87$\times$10$^{-3}$, ${c}_{4}$=1.42$\times$10$^{-4}$, ${c}_{5}$=$-$2.13$\times$10$^{-4}$,  and $zp$=7.5. We then determine the  surface gravity of the stars  by the equation (Prsa et al. 2016):

\begin{equation}
{log(\emph{g}})=4.438+log(\frac{\mathrm{M_{*}}}{\mathrm{M_{\sun}}}) -2log(\frac{\mathrm{R_{*}}}{\mathrm{R_{\sun}}}).
\end{equation}

The photometric calibration of T$_\textrm{\footnotesize{eff}}$ (Eq. 3)  was established before the launch of the Gaia mission; the pre-launched Gaia bandpasses were somewhat off from their later derived values, and the estimated temperatures from this relation are unlikely to be reliable. However, since we use these estimates only  as initial values to start the iterative model-fit pipeline, the photometric values derived above should be sufficiently accurate for our purposes. As shown in Section 6, our determined effective temperatures show strong correlation with color and absolute magnitude in the HR diagram, indicating their high precision, even as only approximate initial values of T$_\textrm{\footnotesize{eff}}$ are used.

Figure 11 shows the photometric surface gravity versus photometric effective temperature of the 3745 stars in our sample. Overplotted are isochrones for a stellar age of 10 Gyr  and metallicities [M/H]=$-$2.0, $-$1.0, and +0.0 dex,  each with three different values of [$\alpha$/Fe], all from the Dartmouth stellar evolutionary model grid (Dotter et al. 2008). The solar-metallicity isochrone with [$\alpha$/Fe]=0 appears to be offset towards higher values of  T$_\textrm{\footnotesize{eff}}$ and log g, while the solar-metallicity isochrones with higher values of $\alpha$-element enhancement, i.e.,  [$\alpha$/Fe]=+0.2 and +0.4, are better matched with the average near-solar metallicity, disk stars (see Figure 31). However, it is difficult to properly interpret  the location of isochrones relative to the stars' distribution  because  the photometric temperatures and surface gravities, and also isochrones  may be inconsistent to some extent with the true parameter values. Figure 12 displays the same plot as shown in Figure 11, but color-mapped based on the metallicity class of the stars. Although the empirical photometric relations have not been tested for  low-metallicity stars ([M/H]$<$$-$0.6 dex) that make up about 6$\%$ of our sample, the  metal-poor stars (MC$\geq$6) are appropriately located in the low-metallicity region of the diagram. Therefore, we believe the photometric estimates provide acceptable initial values for the low-metallicity stars to start the pipeline.

§.§ Fitting Pipelines

Our main goal is to develop a synthetic-fit pipeline that can be accomplished  in a reasonable amount of time. The performing of  a minimization process over 860,000 grid points is extremely time-consuming and exploring all the grid points is a very inefficient process. To this end, our pipeline carries out multiple grid searches  over limited ranges of parameters, beginning with the derived initial values as described above. This procedure significantly reduces the required number of computing processors (CPU or central processing unit) and  the operation time, providing an efficient technique for large, low-resolution surveys.

We examine two different methods and analyze the resulting best-fit parameter values and their distributions in parameter spaces as well as in photometric and kinematic diagrams. In the first method, that we name the “normal method”,  the residual correlation is not taken into account and we normally include all  the wavelength data points in the fitting.  In the second, that we call the “reduced-correlation method”, the residual correlation between neighboring data points is considered and the minimizations  are performed over a subset of wavelength points that are distanced from each other to reduce the effect of this correlation. Each method is  also  investigated in two ways;  one approach keeps  surface gravity constant  and one treats surface gravity as a free parameter. The results from these methods are then compared to each other and  the influence of the residual correlation and surface gravity on derived best-fit parameter values is addressed separately.

§.§.§ Normal Method

Beginning with the above-described initial parameter values  and rounding each value to the closest grid point of the corresponding parameter, if necessary, and also setting an initial convolution factor C equal to the average value of our previous results (Paper I), the model-fit pipeline carries out  ten  successive minimization steps each of which  follows the same routine as outlined in Section 4.1. Best-fit parameter values inferred at each step are used in the next step as revised initial parameters. The first five steps examine different values of the parameters on a coarser model grid, while the last five steps examine the values on a finer model grid. This pipeline, when surface gravity is allowed to vary, is outlined as below:

* Vary the radial velocity shift $v_{r}$ from $-$500 to 500 km/s, in steps of 10 km/s.

* Vary  the convolution factor C from 0.1 to 8, in steps of 0.1.

* Vary T$_\textrm{\footnotesize{eff}}$ within $\pm$200 K, in steps of  100 K, around its corresponding initial value.

* Simultaneously vary  [M/H] within $\pm$0.5 dex, in steps of 0.1 dex,   [$\alpha$/Fe] within $\pm$0.3 dex, in steps of 0.05 dex,  and log g within $\pm$0.2 dex, in steps of 0.1 dex, around their respective initial values.

* Setting all parameters derived from previous steps as new initial values, iteratively repeat all the steps from 1 to 4  until the parameter values converge.

* Setting all parameters derived from the step 5 as new initial values, vary $v_{r}$ from $-$500 to 500 km/s, in steps of 5 km/s.

* Vary  C from 0.05 to 8, in steps of 0.05.

* Vary T$_\textrm{\footnotesize{eff}}$ within $\pm$100 K, in steps of  50 K, around its corresponding initial value.

* Simultaneously vary  [M/H] within $\pm$0.25 dex, in steps of 0.05 dex,   [$\alpha$/Fe] within $\pm$0.1 dex, in steps of 0.025 dex,  and log g within $\pm$0.1 dex, in steps of 0.05 dex, around their respective initial values.

* Setting all parameters derived from the steps 6-9 as new initial values, iteratively repeat these steps  until the parameter values converge.

In the case when surface gravity is regarded as a fixed parameter, the pipeline executes the same steps as above, except for the steps 4 and 9, in which  log g is not allowed to change and its values is kept constant and equal to its initial value.

As it happens, surface gravity plays an important role in low-resolution spectral fitting, although it has the least effect on the spectral shape as compared to the other parameters, notably effective temperature and metallicity. Generally, the accurate determination of  surface gravity, even in high-resolution analyses, has proved to be challenging. Based on our careful grid search,  when this parameter is allowed to vary over a large range of values, it tends to diverge and find best-fit values near the limits of the range for some stars, which do not represent physically realistic solutions knowing that our objects are dwarf stars. The other parameters are accordingly adjusted to this value of log g and unlikely fitted to the correct models. For this reason, we minimize this effect by limiting the range of values, within which log g is free to change in each minimization step, i.e., at most $\pm$0.2 dex.  Nevertheless, there is still a considerable  variation of surface gravity during the model fitting, which has a significant  influence on our final  results and the distributions of stars in different  diagrams (Section 6).  The best treatment for this problem, particularly in low-resolution spectroscopy,  is to use the acceptable values of surface gravity determined using independent methods (e.g., photometry) and keep this parameter fixed in the entire model-fit pipeline. Unfortunately, reliable estimates of  log g are not available for our stars, and this effect remains as a source of uncertainty  for our low-resolution analysis.

§.§.§ Reduced-Correlation Method

Due to the finite resolution of astronomical spectrographs, the instrumental line-spread function (LSF, that describes the distribution of a monochromatic light source at the local plane) is usually sampled over several pixels, resulting in a broadening of observed spectral features. Consequently, adjacent pixels or wavelength datapoints do not represent the independent samples of the true spectrum, and a discrepancy between a synthesized and observed spectral feature may yield a correlated residual that spreads over multiple datapoints (Czekala et al. 2015, hereafter C15). As emphasized in C15 (and pointed out by the first author through our private communication), the Poisson errors (or noises) from astrophysical or instrumental effects inserted on the observed spectrum are likely uncorrelated, but it is the models that do not have enough degrees of freedom to fit the data accurately. A minor mismatch between the model and observed spectra can create correlated fit residuals that propagate over a characteristic scale comparable to the instrumental broadening.

To account for this effect, C15 have developed a formulation of $\chi$$^\textrm{\footnotesize{2}}$ minimization, including  a nontrivial covariance matrix with off-diagonal terms that measures the correlation between any pair of wavelength residuals, which is iterated using Markov Chain Monte
Carlo (MCMC) simulations (Starfish[http://iancze.github.io/Starfish]). Although this framework can perfectly apply
for detailed analyses of a small number of stars at a time,  the large size of our spectroscopic dataset, combined with the time-consuming and computationally expensive MCMC methods, made this procedure too prohibitive to run on available computer resources. As an alternative, we  introduce a significantly quicker  approach to reduce the residual correlation between adjacent datapoints, as described below.

Figure 13 compares the observed spectra (red) with their corresponding  flux-renormalized best-fit models (blue), at unshifted wavelengths,  for four typical  stars  spanning a wide range of metallicities with parameters: T$_\textrm{\footnotesize{eff}}$=3300, 3500, 3500, and 3400 K, [M/H]=+0.3, 0.0, $-$0.6, and $-$1.1 dex, [$\alpha$/Fe]=+0.2, +0.05, +0.275, and +0.325 dex, and  log g=5.35, 5.0, 5.3, and 5.35 dex, respectively. The respective residuals R$_{\textrm{i}}$  between these spectra and their best-fit models, as described by:

\begin{equation}
\mathrm{R}_{\mathrm{i}}=   {
{\mathrm{F}_{\mathrm{obs}}(\lambda_{\mathrm{i,unshift}})  - \mathrm{F}_{\mathrm{mod}}(\lambda_{\mathrm{i,shift}})}},
\end{equation}

which is evaluated at each  datapoint, are also shown in the lower panels. The autocorrelation length scale (ACL) of  the residual array  R$_{\textrm{i}}$ can be determined using the  autocorrelation function (ACF), i.e., from the trend in the autocorrelation value as a function of the pixel lag number. The ACF of the residual arrays from the four stars in Figure 13 are shown in Figure 14.  Here we calculate their residual ACF relative to their best-fit model from the “normal method” when surface gravity is a free parameter[In general, there is a smaller distance between the observed spectra and their best-fit models derived from the pipeline in which log g is a free parameter, which leads to apparently lower values of $\chi$$^\textrm{\footnotesize{2}}$, as compared to those inferred from the pipeline in which log g is constant. However,  there may be a strong degeneracy effect underlying the resulting parameter values from the former method for some stars, which is difficult to perceive just by comparing the observed spectrum and and its best-fit model. Nevertheless, we find no considerable difference in the  estimated ACL using the two approaches.]. Since the stars were observed by similar instrumental setups, all the plots show more and less similar decreasing trends with pixel lag number.  The  ACL  is defined by the lag numbers  where the ACF of the residuals are reduced by about a third ($\sim$20-40%) of the value at lag number=0. For our observed spectra, we find that ACL$\sim$10 pixels in the vast majority of cases. In  other words,  the residual at each point is significantly propagated over $\sim$10 neighboring  datapoints, which means that the residuals of any 10 adjacent datapoints are highly correlated. There is a relationship between  the ACL and the spectral resolution of the related spectrograph: the higher resolution, the less instrumental broadening, and the shorter ACL.

In the “reduced-correlation method”, we therefore perform each minimization routine over every 10-th wavelength datapoint, i.e.,  selecting a subset of one tenth of the total datapoints that show little to no correlation and thus efficiently sample our low-resolution spectra. Starting from the first datapoint, we select every 10-th wavelength, and apply the model-fit pipeline to these selected wavelengths to derive the parameter values. We then repeat the same process, but starting from the second  datapoint and derive a second set of parameter values. This approach is repeated ten times up to the tenth datapoint and each time a new set of parameter values are inferred. We lastly take an average of the ten values obtained from each step and report this average as our final best-fit estimate.

Since we run this procedure on each spectrum after running the “normal method” described in Section 4.3.1, we use the final best-fit values from the “normal method”  as initial values in this “reduced-correlation method”.  We  again perform the “reduced-correlation method” for two cases of surface gravity, i.e., when it is a free or a constant parameter, and the parameter estimates  from the “normal method” using each  gravity-modeling approach are then used as initial values in the corresponding approach of the "reduced-correlation method”. We also treat  the   secondary parameters  as fixed parameters in the entire process, equal to the well-determined values obtained from the “normal method”. The pipeline, in which surface gravity is a variable parameter, is thus shortened to only three steps, as follows:

* Vary T$_\textrm{\footnotesize{eff}}$ within $\pm$100 K, in steps of  50 K, around its corresponding initial value.

* Simultaneously vary  [M/H] within $\pm$0.25 dex, in steps of 0.05 dex,   [$\alpha$/Fe] within $\pm$0.1 dex, in steps of 0.025 dex,  and log g within $\pm$0.1 dex, in steps of 0.05 dex, around their respective initial values.

* Setting all parameters derived from  the steps 1 and 2 as new initial values, iteratively repeat these steps  until the parameter values converge.

In the case when surface gravity is a fixed parameter, the pipeline performs the three steps above, except for the steps 2 and 3 in which log g is kept  constant equal to its initial value.

Comparison between the observed spectra (red) and their respective  flux-adjusted, best-fit models (blue), at unshifted wavelengths,  for four stars with metallicities  ranging from +0.3 dex to $-$1.10 dex. The associated residuals between the spectra and their best-fit models are also shown in the lower panels. The parameters of these stars are reported in the legends.

Residual Autocorrelation versus lag numbers for the same stars as those shown in Figure 13.

$\chi$$^\textrm{\footnotesize{2}}$ maps in two-dimensional parameter planes for the star PM J01305+2008  with T$_\textrm{\footnotesize{eff}}$=3000 K, [M/H]=$-$0.05 dex, [$\alpha$/Fe]=+0.25 dex, and  log g=5.0 dex. The best-fit grid point is shown by a white dot. The number of wavelength datapoints used in the $\chi$$^\textrm{\footnotesize{2}}$ calculation is 682.

$\chi$$^\textrm{\footnotesize{2}}$ maps in two-dimensional parameter planes for the star PM J03440+2730 with T$_\textrm{\footnotesize{eff}}$=3400 K, [M/H]=$-$0.9 dex, [$\alpha$/Fe]=+0.425 dex, and  log g=5.45 dex. The best-fit grid point is shown by a white dot. The number of wavelength datapoints used in the $\chi$$^\textrm{\footnotesize{2}}$ calculation is 1063.

$\chi$$^\textrm{\footnotesize{2}}$ maps in two-dimensional parameter planes for the star PM J03352+2055  with T$_\textrm{\footnotesize{eff}}$=3400 K, [M/H]=$-$1.1 dex, [$\alpha$/Fe]=+0.325 dex, and  log g=5.35 dex. The best-fit grid point is shown by a white dot. The number of wavelength datapoints used in the $\chi$$^\textrm{\footnotesize{2}}$ calculation is 734.

$\chi$$^\textrm{\footnotesize{2}}$ maps in two-dimensional parameter planes for the star PM J06069+1706  with T$_\textrm{\footnotesize{eff}}$=3300 K, [M/H]=$-$1.45 dex, [$\alpha$/Fe]=+0.375 dex, and  log g=5.45 dex. The best-fit grid point is shown by a white dot. The number of wavelength datapoints used in the $\chi$$^\textrm{\footnotesize{2}}$ calculation is 606.

§ ERROR ANALYSIS

The uncertainty in our inferred parameters originates from distinct sources as follows:

* Uncertainties in observed spectra, e.g.,  instrumental errors and statistical noise from photon count. The observational uncertainties  are taken into account as an inverse weight in the $\chi$$^\textrm{\footnotesize{2}}$  formalism (Eq.1).

* Uncertainties in the interpolated synthetic spectra. Since the corresponding  spectra calculated from  the atmosphere code are not available yet (as their computations would be exceedingly time consuming), no comparison can be made to quantify the difference between the interpolated and calculated spectra. Nevertheless, we assume that these interpolated model spectra represent a fair approximation of the calculated ones.

* Uncertainties in the flux-renormalization of synthetic spectra relative to the observed spectrum of interest. This is difficult to quantify using the present $\chi$$^\textrm{\footnotesize{2}}$ method. However, some other methods such as  Starfish (see above) provide a routine through MCMC simulations (C15) that can account for any uncertainties induced by the modification of model spectra. Unfortunately, such methods would be extremely computationally intensive for our large M dwarf/subdwarf sample, and  imperfections in the flux-renormalization of synthesized models thus remain as an underlying source of error in the inferred parameter.

* Errors from the residual correlation between wavelength datapoints. Such an error factor  can be significantly reduced by using the “reduced-correlation method” described in Section 4.3.2, although this does not make a considerable difference in the distribution of the derived parameters, which is likely due to the low spectral resolution of our spectra.

* Errors  from  the degeneracy between physical parameters. As explained  in Section 3.4 (and also shown in Section 6), the degeneracy effect from log g and each of chemical parameters, i.e., [M/H] and  [$\alpha$/Fe],  can cause significant errors in the resulting parameters. The best solution is to keep surface gravity fixed and equal to the values inferred from independent techniques. There is also a weaker degeneracy between [M/H] and  [$\alpha$/Fe], which may cause errors in the derived parameter values.

* Uncertainties due to possibly undetected unresolved binaries. These binary systems have a composite spectrum, which can make significant errors in the derived stellar parameters.

* Uncertainties associated with  limitations of models  to reproduce all observed spectral details. These are caused by the incompleteness of available opacity data, and by simplifying assumptions  in modeling stellar atmosphere, e.g.,  the local thermodynamic equilibrium (LTE) that overlooks  non-LTE (NLTE) effects, which may cause significant abundance errors (Bergemann & Nordlander 2014; Nordlander et al 2017) or  simplifications used in the Mixing Length Theory (Kippenhahn & Weigert 1990) connected with stellar convection. Such uncertainties are much more difficult to quantify and our error analysis does not include the systematic uncertainties due to the disagreement between models and observations. Nevertheless, we find systematic trends between parameters that are unlikely to have physical origins and may be attributed to the inadequacy of model atmospheres.

Figures 15-18 show the  $\chi$$^\textrm{\footnotesize{2}}$ maps in six two-dimensional parameter spaces for four typical stars in different metallicity regimes. The parameter values are estimated using the “normal method” in which surface gravity is considered as a free parameter. To create these two-dimensional maps, we vary two parameters that are shown in X and Y axes, while keeping the two other  parameters fixed and equal to their best-fit estimates (otherwise, the  $\chi$$^\textrm{\footnotesize{2}}$ would be a multivalued function). The best-fit grid point is presented by a white dot  in each plot. The number of wavelength  datapoints ($\sim$the degrees of freedom) used in the  $\chi$$^\textrm{\footnotesize{2}}$ calculation for each spectrum (i.e., after removing the regions including  the H$_\textrm{\footnotesize{$\alpha$}}$ emission line of hydrogen, telluric absorption lines, and  spectral artifacts as well as highly noisy ranges)  is shown in the respective captions. The best-fit grid points are located in the regions with lowest $\chi$$^\textrm{\footnotesize{2}}$ values, which shows the consistency of our model-fit pipeline, though in some cases the best-fit parameter is off-center from the $\chi$$^\textrm{\footnotesize{2}}$ distribution, which suggests small systematic errors in the minimization algorithm, likely due to the limited resolution of the model grid. We also note that the reduced $\chi$$^\textrm{\footnotesize{2}}$ ($\chi$$^\textrm{\footnotesize{2}}$/degrees of freedom) is often near or even less than 1, indicating a potential issue with the flux uncertainties (E$_{\textrm{obs}}$) that, in addition to Poisson noises,  incorporate terms propagated from uncertainties in the flux calibration. These errors are partially,  but not fully, canceled by flux-renormalization (Section 4.1), and the flux uncertainty values are therefore likely to be somewhat overestimated.

There is a relatively tight correlation between [M/H] and [$\alpha$/Fe], directed  from high-[M/H] and low-[$\alpha$/Fe] to low-[M/H] and high-[$\alpha$/Fe] values. As explained in Section 3, these two chemical parameters change  the spectral morphology in roughly the same way, but at different levels. This anti-correlation between [M/H] and [$\alpha$/Fe], however, suggests that the combined parameter [$\alpha$/Fe]+[M/H] should be significantly better constrained than either one individually. This anti-correlation  also appears as  tilted distributions  of stars with nearly  the same best-fit values of log g (even with different effective temperatures) in the [$\alpha$/Fe] versus [M/H] diagram, no matter which gravity-modeling approach is used (Section 6.3). As can be seen from these figures, T$_\textrm{\footnotesize{eff}}$ is much less correlated with [M/H] and [$\alpha$/Fe], while the correlation between log g and  the two chemical parameters is  prominent.

The joint error distribution of  T$_\textrm{\footnotesize{eff}}$, [M/H], [$\alpha$/Fe], and log g can be estimated using the $\chi$$^\textrm{\footnotesize{2}}$ values of a model grid that spans over ranges  adequately wider  than the expected parameter errors[Critical values of the $\chi$$^\textrm{\footnotesize{2}}$ distribution with various degrees of freedom are tabulated and available in the literature. However, these values do not apply to our analysis, yielding too small parameter errors.]. Focusing on the “normal method”,  we  find the “joint confidence grid”,  that contains the  models with $\chi$$^\textrm{\footnotesize{2}}$ values within 20% deviation from the minimum value, is large enough for calculating the errors of the derived parameters. The deviation from the minimum value of $\chi$$^\textrm{\footnotesize{2}}$ by 20% typically spans parameter ranges around best-fit values within T$_\textrm{\footnotesize{eff}}$$\sim$$\pm$150-160 K, [M/H]$\sim$$\pm$0.4-0.45 dex, [$\alpha$/Fe]$\sim$$\pm$0.2-0.25 dex, and  log g$\sim$$\pm$0.25 dex (if using the variable-gravity approach) that are considerably larger than the average errors of our results from Paper I. In addition, the variation of the inferred parameter values from the present work in the abundance diagram of [$\alpha$/Fe] versus [M/H] (Section 6.3) or abundance-velocity  spaces (Section 9) are quite comparable to those of other studies, which indicates that the expected parameter errors must be  less than the above values. In Section 9, we also identify substructure such as streams in abundance-velocity spaces that would not be possible if the resulting chemical parameters had large uncertainties. The 20% deviation is a conservative choice,  in comparison with  the 5%  departure chosen in the high-resolution analysis of Rajpurohit et al. (2014) and  the 10%  departure chosen in our previous low-resolution analysis of Paper I.  The parameter errors  are then determined by the standard deviation of the individual confidence ranges that are  obtained by projecting the  four-dimensional joint confidence grid on  the corresponding one-dimensional  parameter spaces.

In the present study, the inferred parameter values from the “reduced-correlation method” are only used to compare with those from the “normal method”  and  to investigate how the results change when using independent wavelength datapoints.  The standard deviation of the ten measurements from the mean value for each parameter is calculated to show the scatter of different estimates from different sets of nearly uncorrelated datapoints.

§ RESULTS: NORMAL METHOD

We apply the “normal method” to the 3745 M dwarfs/subdwarfs for two approaches when surface gravity is kept constant (as specified by “const-grav”) and when surface gravity is allowed to vary (as specified by “var-grav”). In the following sections, we examine the distribution of the  inferred parameter values using different diagrams and address the effect of surface gravity variation on the model-fit values of the other parameters.

§.§ Effective Temperature Distribution

Figure 19 compares the photometric temperatures with  the spectroscopic model-fit values inferred from the two gravity-modeling  approaches, which are randomized within  the step size of the temperature gird. Both plots present a  similar dispersion around the 1:1 line (black line), though there is a slightly better consistency between the  photometric temperatures and the model-fit values  from the constant-gravity approach (top panel) compared those from the variable-gravity approach  (bottom panel). These panels indicate  an average offset (downward shift relative to the 1:1 line) of $\sim$120-140 K for stars with T$_\textrm{\footnotesize{eff}}$$\gtrsim$3700 K and  $\sim$50 K  for stars with T$_\textrm{\footnotesize{eff}}$$\lesssim$3100-3200 K. The largest discrepancies occur in metal-poor stars with a typical offset (upward shift relative to the 1:1 line) of $\sim$170 K between two temperature measurements, which suggests that the photometric temperatures for low-metallicity stars may not be as reliable; this is likely because temperature-color relationships have been largely calibrated and tested using metal-rich stars. We thus assume that the spectroscopic model-fit values are more reliable in this case.

§.§ Gaia HR Diagram

We aim to examine the variation of parameter values in the  HR diagram using Gaia magnitudes and parallaxes. However, the position of stars in  such a diagram is affected by the precision of Gaia astrometric measurements which is quantified by a quality flag, i.e., the  Renormalized Unit Weight Error (RUWE)[The detailed description of RUWE can be found in the Gaia  DPAC public documentation: https://gea.esac.esa.int/archive/documentation/GDR2/Gaia_ar chive/chap_datamodel/sec_dm_main_tables/ssec_dm_ruwe.h tml]. The RUWE is a measure of the astrometric noise and is expected to be  around 1.0 for a well-behaved measurement of a single star  where the single-star model provides a good fit to the astrometric observations, despite the primary color or magnitude. On the other hand, a RUWE value significantly greater than 1.0 (say, >1.4) indicates that the source is non-single or otherwise problematic for the astrometric solution. The RUWE is set to null for sources with only a two-parameter solution, since this value would be difficult to interpret in such cases. Some studies have used this parameter as an indication of binarity. For example, Ziegler et al. (2020) showed  that 86% of their planet candidate sample with RUWE>1.4 had companions undetected by the Gaia survey, but observed in high-resolution follow-up. The effect and application of Gaia RUWE values can be found in a number of papers (e.g., Ziegler et al. 2020, Belokurov et al. 2020; Wood et al. 2021; Kervella et al. 2022).

Figure 20 shows the RUWE values of our 3745 stars versus their distances  (in  logarithmic scale). We find high values of RUWE for a sizable fraction of our sample. In total 1066 stars, i.e., 28% of the sample, have  RUWE>1.4, of which 459 stars have 1.4<RUWE$\leq$2, 373 stars have 2<RUWE$\leq$10, and 234 stars have RUWE>10.  As seen from the figure, the high-RUWE stars are mostly  located at lower distances ($\lesssim$50 pc), presenting a high level of binarity. In general, unresolved binaries are easier to detect from astrometric deviations when the stars are relatively nearby. Our sample consists of a large number of nearby stars (within 25 pc) compared with typical  samples from Gaia that usually include stars at  larger distances, and consequently,  unresolved systems might not show up as well in the astrometric solution. Moreover, our spectroscopic survey selected nearby stars based on (pre-Gaia) photometric distances, and   binary stars  are more likely to be selected in a “nearby” sample based on a photometric distance cut because unresolved binaries are overluminous at a given color and thus appear to be closer. Figure 21 displays the Gaia HR diagram,  M$_{\footnotesize{G}}$ (absolute G magnitude) versus  G$_\textrm{\footnotesize{BP}}$-G$_\textrm{\footnotesize{RP}}$ (top panel) and M$_{\footnotesize{G}}$ versus  G -G$_\textrm{\footnotesize{RP}}$ (bottom panel), of  two subsets: 2679 stars with RUWE<1.4 (as shown in blue) and 1066 stars with RUWE>1.4 (as shown in black). Clearly, the majority of high-RUWE stars reside towards the upper end of the distribution, which is likely due to being over-luminous as a result of their  binarity.

Figure 22 presents the Gaia HR diagram, M$_{\footnotesize{G}}$ versus  G$_\textrm{\footnotesize{BP}}$-G$_\textrm{\footnotesize{RP}}$,  of the full set of 3745 stars (top panels) and a cleaner subset of  2679 stars with RUWE<1.4  (bottom panels).  All these plots are color-mapped based on the  effective temperatures derived from both  procedures when  log g is  fixed (left panels) and when  log g is variable (right panels); in this case we find that there is no significant difference between the inferred temperatures from these two  gravity-modeling approaches.  Obviously, the effective temperature of stars is consistent  with their Gaia color and absolute magnitude, and this close agreement is a strong indication that our model-fit pipeline produces relatively high precision values for T$_\textrm{\footnotesize{eff}}$ in most cases. The absence of  high-temperature, low-metallicity stars  is a result of sample selection effect;  among metal-poor, low-mass stars, those with T$_\textrm{\footnotesize{eff}}$>3600 K look significantly bluer and tend to be classified as K dwarfs, and are thus excluded. In addition, low-temperature, metal-poor M subdwarfs are often  too faint to be observed using the selected telescopes and spectrographs (see Tables 1 and 2 in Paper I for more details), and for this reason,  are not present  in our sample. Both effects explain why metal-poor M subdwarfs are found over a much narrower range of T$_\textrm{\footnotesize{eff}}$ values (3000$\lesssim$T$_\textrm{\footnotesize{eff}}$$\lesssim$3600 K), while the range of T$_\textrm{\footnotesize{eff}}$ values for the more metal-rich M dwarfs is much broader  (2800$\lesssim$T$_\textrm{\footnotesize{eff}}$$\lesssim$4000 K). Figure 23 depicts the HR diagram, but with another form, i.e.,  M$_{\footnotesize{G}}$ versus  G-G$_\textrm{\footnotesize{RP}}$, of the same samples, color-mapped by the same effective temperatures, as shown in the respective panels of Figure 22. The impact of RUWE appears to be more important in the HR diagram with the color  G-G$_\textrm{\footnotesize{RP}}$ than the color  G$_\textrm{\footnotesize{BP}}$-G$_\textrm{\footnotesize{RP}}$.

[Fixed surface gravity]

[Variable surface gravity]
Comparison between the photometric temperatures and the inferred model-fit values for the 3745 stars, when surface gravity is considered as a fixed parameter  (top panel) and when surface gravity is considered as a free parameter (bottom panel), using the normal method.

RUWE versus distance of the 3745 stars in  logarithmic scale.

HR diagram,  M$_\textrm{\footnotesize{G}}$ versus G$_\textrm{\footnotesize{BP}}$ - G$_\textrm{\footnotesize{RP}}$ (top panel) and M$_\textrm{\footnotesize{G}}$ versus G - G$_\textrm{\footnotesize{RP}}$ (bottom panel), of the 3745 stars divided into two subsets: 2679 stars with RUWE<1.4 (blue) and 1066 stars with RUWE>1.4 (black).

[Fixed surface gravity]

[Variable surface gravity]

[Fixed surface gravity]

[Variable surface gravity]
HR diagram,  M$_\textrm{\footnotesize{G}}$ versus G$_\textrm{\footnotesize{BP}}$ - G$_\textrm{\footnotesize{RP}}$,  of the 3745 stars (top panels) and a subset of 2679 stars with RUWE<1.4 (bottom panels),  when surface gravity is set as a fixed parameter (left panels) and when surface gravity is set as a free parameter (right panels), color-mapped based on  effective temperatures, using the normal method.

[Fixed surface gravity]

[Variable surface gravity]

[Fixed surface gravity]

[Variable surface gravity]
HR diagram,  M$_\textrm{\footnotesize{G}}$ versus G - G$_\textrm{\footnotesize{RP}}$,  of the 3745 stars (top panels) and a subset of 2679 stars with RUWE<1.4 (bottom panels),  when surface gravity is set as a fixed parameter (left panels) and when surface gravity is set as a free parameter (right panels), color-mapped based on  effective temperatures, using the normal method.

[Fixed surface gravity]

[Variable surface gravity]

[Fixed surface gravity]

[Variable surface gravity]
HR diagram,  M$_\textrm{\footnotesize{G}}$ versus G$_\textrm{\footnotesize{BP}}$ - G$_\textrm{\footnotesize{RP}}$, of the 3745 stars (top panels)  and a subset of 2679 stars with RUWE<1.4 (bottom panels),  when surface gravity is set as a fixed parameter (left panels) and when surface gravity is set as a free parameter (right panels), color-coded based on five groups with different  metallicity ranges, using the normal method. The number of stars in each metallicity group, N, is panel “a”:  N$_\textrm{{Black}}$=105,  N$_\textrm{{Yellow}}$=2849, N$_\textrm{{Blue}}$=593, N$_\textrm{{Red}}$=85, and N$_\textrm{{Green}}$=113, panel “b”:  N$_\textrm{{Black}}$=784,  N$_\textrm{{Yellow}}$=2278, N$_\textrm{{Blue}}$=477, N$_\textrm{{Red}}$=107, and N$_\textrm{{Green}}$=99,    panel “c”:  N$_\textrm{{Black}}$=82,  N$_\textrm{{Yellow}}$=1943, N$_\textrm{{Blue}}$=461, N$_\textrm{{Red}}$=83, and N$_\textrm{{Green}}$=110,     panel “d”:  N$_\textrm{{Black}}$=502,  N$_\textrm{{Yellow}}$=1599, N$_\textrm{{Blue}}$=378, N$_\textrm{{Red}}$=103, and N$_\textrm{{Green}}$=97.

[Fixed surface gravity]

[Variable surface gravity]

[Fixed surface gravity]

[Variable surface gravity]
HR diagram,  M$_\textrm{\footnotesize{G}}$ versus G - G$_\textrm{\footnotesize{RP}}$, of the 3745 stars (top panels)  and a subset of 2679 stars with RUWE<1.4 (bottom panels),  when surface gravity is set as a fixed parameter (left panels) and when surface gravity is set as a free parameter (right panels), color-coded based on five groups with different  metallicity ranges, using the normal method. The number of stars in each metallicity group is the same as that shown in Figure 23 associated with the respective panels.

[Fixed surface gravity]

[Variable surface gravity]

[Fixed surface gravity]

[Variable surface gravity]
HR diagram for  the stars with T$_\textrm{\footnotesize{eff}}$$\leq$3550 K, divided into two groups with [M/H]$\geq$$-$0.5 dex (top panels) and  [M/H]<$-$0.5 dex  (bottom panels), when the parameters are determined with fixed surface gravity  (a: 2714 stars and c: 162 stars) and with variable surface gravity (b: 2701 stars and d: 143 stars), color-mapped with associated values of the combined chemical parameter [M/H]+[$\alpha$/Fe], using the normal method. The values of RUWE are not taken into account in this figure.

Comparison between the inferred model-fit parameters, T$_\textrm{\footnotesize{eff}}$, [M/H], and  [$\alpha$/Fe]  when surface gravity is a fixed parameter  and their corresponding values when surface gravity is a free parameter, for the 2829 stars with T$_\textrm{\footnotesize{eff}}$ $\leq$3550 K in both gravity cases, color-mapped based on the surface gravity values from the former approach  (top panels) and from the latter approach (middle panels), using the normal method. The bottom panels show the histograms of the differences between corresponding values inferred from the two gravity-modeling approaches.

[Fixed surface gravity]

[Variable surface gravity]
Histogram of estimated metallicities for the stars with  T$_\textrm{\footnotesize{eff}}$$\leq$3550 K and [M/H]$\geq$$-$0.5 dex from the approach when surface gravity is constant  (2714 stars, top panel) and when surface gravity is varying (2701 stars, bottom panel), using the normal method.

[Fixed surface gravity]

[Variable surface gravity]

[Fixed surface gravity]

[Variable surface gravity]

[Fixed surface gravity]

[Variable surface gravity]
[$\alpha$/Fe] vs. [M/H] diagram for the stars with T$_\textrm{\footnotesize{eff}}$$\leq$3550 K and $-$2.5<[M/H]<$+$0.5 dex when surface gravity is a fixed parameter (left panels, 2868 stars) and when surface gravity is a free parameter (right panels, 2668 stars), color mapped based on metallicity class (top panels), T$_\textrm{\footnotesize{eff}}$ (middle panels), and log g (bottom panels), using the normal method.

[Fixed surface gravity]

[Variable surface gravity]

[Fixed surface gravity]

[Variable surface gravity]

[Fixed surface gravity]

[Variable surface gravity]
[$\alpha$/Fe] vs. [M/H] diagram for the stars with T$_\textrm{\footnotesize{eff}}$$\leq$3550 K and $-$0.5$\leq$[M/H]<$+$0.5 dex when surface gravity is a fixed parameter (left panels, 2698 stars) and when surface gravity is a free parameter (right panels, 2511 stars), color mapped based on metallicity class (top panels), T$_\textrm{\footnotesize{eff}}$ (middle panels), and log g (bottom panels), using the normal method.

Top panel: photometric surface gravity versus inferred model-fit effective temperature for 2877 stars  with T$_\textrm{\footnotesize{eff}}$$\leq$3550 K, color-mapped based on the metallicity of  stars. All best-fit values are derived using the normal method when surface gravity is fixed and equal to the photometric values. Bottom panel: the same plot as shown in the top panel but for stars with [M/H]$\geq$$-$0.5 dex along with solar-metallicity isochrones related to  three different values of  [$\alpha$/Fe]=+0.0, +0.2, and +0.4 dex (solid, dashed and dotted black lines, respectively) from the Dartmouth stellar evolutionary models.

[Fixed surface gravity]

[Variable surface gravity]
HR diagram for the stars with  T$_\textrm{\footnotesize{eff}}$$\leq$3550 K from the approach when surface gravity is  constant  (2877 stars, top panel) and when surface gravity is varying (2851 stars, bottom panel), color-mapped based on surface gravity values, using the normal method.

Figures 24 and 25 show the same HR diagram panels as in Figures 22 and 23, respectively, but color-coded according to five groups  of stars with different  metallicity ranges, as noted in the legends. Results from the constant-gravity and variable-gravity methods mostly differ in the range of best-fit values, with the variable-gravity yielding a wider range of metallicities, notably at the metal-rich end. On the other hand, both methods show a clear stratification in the color-absolute magnitude diagram. This signifies that our pipeline can separate different metallicity populations with high precision, whether surface gravity is a fixed or a free parameter in the fitting process.  While the two gravity approaches do not  make a considerable  difference  in the resulting model-fit effective temperatures, the derived metallicities, particularly in the near-solar metallicity regime,  are  significantly dependent on which approach is used. As shown in Sections 3 and 5, surface gravity  is more correlated with  metallicity  than effective temperature, and as a result, the best-fit values of metallicity are  remarkably affected by the variation of   surface gravity (Sections 6.3 and 6.4).

Both Figures 24 and 25 also indicate that, in either gravity approach, only very few  observed  spectra are matched with high-temperature, high-metallicity models (T$_\textrm{\footnotesize{eff}}$>3550 K and [M/H]>+0.1), and our high-temperature stars (whose best-fit models resemble the atmosphere of  late K dwarfs with a weaker level of molecular opacity, rather than M dwarfs) are mostly concentrated in the metallicity range $-$0.5$\leq$[M/H]$\leq$+0.1 dex. This systematic concentration  is found to be  independent of the selected initial values, the $\chi$$^\textrm{\footnotesize{2}}$ formalism (e.g., whether the observed flux errors are included in the $\chi$$^\textrm{\footnotesize{2}}$ expression or not), the number of iterative processes,  and the extent of the examination grid. Therefore, the lack of high-temperature, high-metallicity stars in our sample  is likely due to the inconsistency between models and observed spectra, rather than issues in the model-fit pipeline.

Figure 26 presents the same HR diagrams  as  shown in the top panels of Figure 24, but for stars with T$_\textrm{\footnotesize{eff}}$$\leq$3550 K (that are believed to  have more accurate parameters),  for two metallicity ranges, i.e., [M/H]$\geq$$-$0.5 dex (top panels) and [M/H]<$-$0.5 dex  (bottom panels), and for two gravity approaches, i.e.,  when surface gravity is a fixed parameter (left panels) and when surface gravity is a free parameter (right panels). All these plots are color-mapped based on the combined parameter [M/H]+[$\alpha$/Fe]. As can be seen, the distribution of  [M/H]+[$\alpha$/Fe] values is well stratified in both near-solar metallicity and low-metallicity regimes, which again indicates a high precision in   the best-fit values of  chemical-parameters, no matter which surface gravity approach is used.

§.§ Constant versus Variable Gravity

In this section we examine how the variation of surface gravity  in the model-fitting affects the values of the other parameters. As explained  in Section 6.2, the inferred parameter values of high-temperature stars (T$_\textrm{\footnotesize{eff}}$>3550 K), whose best-fit models are closer to spectra of late-K dwarfs than they are to M dwarfs (i.e., they have relatively  weak molecular bands), are most likely to suffer from large systematic uncertainties as a result. Because of this, we exclude those stars from our present analysis.

Figure 27 compares the best-fit values of T$_\textrm{\footnotesize{eff}}$, [M/H], and [$\alpha$/Fe]  obtained from the two gravity-modeling  approaches, as specified by “const-grav” and “var-grav”, for the 2829 stars with T$_\textrm{\footnotesize{eff}}$ $\leq$ 3550. To ease the interpretation, the best-fit values for each parameter are randomized within the corresponding step size of the model grid (otherwise they would fall along the evenly spaced grid points). The  color map in the top panels is based on log g (const-grav) and in the middle panels is based on log g (var-grav).  The bottom panels show the histograms of the differences between the values of the corresponding parameters obtained from the two gravity-modeling approaches.  Although there is a significant change  in the best-fit values of  log g when being treated as a free parameter  relative to those constant values (i.e., photometric surface gravities), there appears to be only a small (downward) shift in the distribution of temperatures with respect to the line 1:1, as shown in the left panels “a” or “d”.  The histogram plotted in the bottom left panel “g” is also slightly off relative to the zero value towards higher values by $\sim$10-20 K. This is the reason why there is no significant difference in the variation of T$_\textrm{\footnotesize{eff}}$ in the HR diagram from the two gravity-modeling approaches (Figures 22 and 23), which is consistent with the weak degeneracy effect between effective temperature and the other parameters (Section 3.4).

Clearly, the variation of surface gravity has a  more important effect on the best-fit values of metallicity,  as illustrated in the middle panels of Figure 27.  As shown in Section 3.4, the simultaneous variation of log g and [M/H] can make a substantial degeneracy effect in the inferred model-fit  values, particularly for near-solar metallicity spectra. Accordingly,  once surface gravity becomes a free parameter,  the  surface gravity and metallicity  estimates  of a significant number of stars with near-solar metallicities shift to higher values, while  there is no noticeable change   in these two parameters for metal-poor stars.  This can clearly be seen from the middle panels (“b” and “e”) that show a relatively large fraction of stars which  are shifted  below  the 1:1 line (black line). The histogram in the panel “h” also shows an overall shift of the distribution towards higher values by $\sim$0.12-0.13 dex relative to the values obtained from the fixed-gravity model-fit. During the model-fit process, the values of log g and [M/H] competitively increase in each iteration, and for some stars, may even drift to the upper limits of the  model grid  (i.e., log g=5.5 dex and/or  [M/H]=+0.5 dex), which are unlikely to represent realistic values of these parameters for the stars in our sample. There is also an overall  increase in [$\alpha$/Fe], mainly  associated with near-solar metallicity stars,  which is also owing to the degeneracy effect between log g and [$\alpha$/Fe] (Section 3.4).  This can be observed in  the  right panels  “c” or “f”, that show how the overall distribution of stars  are systematically off below the 1:1 line,  mostly for [$\alpha$/Fe]$\lesssim$+0.4 dex. This effect  yields  a shift in the overall histogram distribution shown in  the panel “i” towards higher values by $\sim$0.03-0.04 dex relative to values obtained with the fixed-gravity method.

The end effect is quite dramatic: in Figure 28, we show the metallicity distribution of stars with [M/H]$\geq$$-$0.5 dex  and T$_\textrm{\footnotesize{eff}}$$\leq$3550 K when surface gravity is kept constant  (2714 stars, top panel) and when surface gravity is allowed to vary (2701 stars, bottom panel). There is a clear shift of metallicities towards higher values when  surface gravity becomes a free parameter in the model fit to synthetic spectra. The constant gravity approach results in  just six stars with metallicities at  the upper limit of the model grid, i.e., [M/H]=+0.5 dex, and  98 stars  with +0.3$\leq$[M/H]<+0.5 dex. On the other hand,  the varying gravity approach yields 181 stars having metallicities at  [M/H]=+0.5 dex and 582 stars with +0.3$\leq$[M/H]<+0.5 dex. The large fraction of stars with  high metallicity values  ([M/H]$\gtrsim$+0.3 dex) in our volume-limited  sample derived from the latter method  is not consistent with the metallicity distribution of more massive (FGK) disk stars  in the solar neighborhood. On the other hand, the metallicity distribution obtained from the constant gravity method  is similar to those of previous studies for nearby dwarf stars (e.g., Rocha-Pinto & Maciel 1998; Haywood 2001; Nordström et al. 2004, hereafter N04; Du et al. 2004; Hejazi et al. 2015; An 2019). This implies that the high-metallicity fits  from the varying gravity approach  most likely suffer from large systematic errors and  do not seem to represent the realistic measurements of metal content.

In Figure 29, we examine the distribution of our stars in the [$\alpha$/Fe]  versus [M/H] diagram, color-mapped by metallicity class (top panels), effective temperature (middle panels), and surface gravity (bottom panels), for the two  gravity-modeling approaches, i.e., when log g is kept fixed (2868 stars, left panels) and when gravity is free to change (2668 stars, right panels).  This time, we also exclude all stars with best-fits values of [M/H]=+0.5 dex,  because these stars have best-fit values in one parameter that are at the edge of the model grid and we assume that all best-fit parameters (notably gravity) may be biased one way or another. The values of [M/H] and [$\alpha$/Fe] for the remaining stars are randomized within the associated step sizes of the model grid, again for ease of interpretation. In order to better present the detailed distribution of the near-solar metallicity (disk) stars,  Figure 30 shows the same panels as those in Figure 29, but zoomed-in over the metallicity range $-$0.5$\leq$[M/H]<$+$0.5 dex.

There is a global trend in the distribution of stars for both gravity-modeling methods: the higher-metallicity stars, on average, have lower values of $\alpha$-element enhancement, while  lower-metallicity stars, on average, tend to have higher values of $\alpha$-element enhancement, which is consistent with current models of  the Milky Way's enrichment history  (e.g., Croswell 1995; Pagel 1997; Chiappini et al. 1999, 2001; Chiappini 2001). The time scales of chemical enrichment, particularly by type II supernovae (SNe II) and type Ia supernovae (SNe Ia),  are critical to establish the evolutionary models for elemental abundance ratios such as [$\alpha$/Fe];  $\alpha$-elements are mostly produced by SNe II over time scales of 3 to 30 Myr, while iron is primarily  generated by SNe Ia  over larger time-scales ranging from 30 Myr to the Hubble time (Pipino & Matteucci 2009). During  the chemical  enrichment  of  the interstellar medium (ISM)   by core-collapse  SNe II,  heavy elements are  highly  rich in  $\alpha$-elements. However, after the “time delay”  ($\sim$ a few billion years, e.g., Heringer et al. 2019), when white dwarfs in binary systems   grow up to  the Chandrasekhar mass by accretion from their companions, which leads to SNe Ia explosions,  the ejected materials start to enrich the ISM with iron, and accordingly  [$\alpha$/Fe]  begins to drop.  This figure  also shows how  stars get redistributed when log g becomes a free parameter, shifting a considerable number of  stars towards  higher values of [M/H], log g and/or higher values of  [$\alpha$/Fe].

There is a very clear correlation between the location of a star in this diagram and the star's metallicity class (MC) estimated from spectral classification, as seen from the color-coding on the top panels. The trend is nearly  independent of the gravity-modeling approach used in the synthetic fitting. Regions with similar values of MC are tilted from higher-[M/H] and lower-[$\alpha$/Fe] to lower-[M/H] and higher-[$\alpha$/Fe] values, which indicates that the metallicity class depends on  both metallicity and $\alpha$-element enhancement. The spectral metallicity classification based on the TiO and CaH molecular bands (Paper I) is thus determined by both chemical parameters, and is not merely a measure  of “metallicity” in the general sense (i.e., of [M/H]). We find that the distribution  of stars with nearly the same values of log g, for both gravity-modeling approaches, are also tilted in nearly the same direction as the distribution of stars with similar MC values, which shows a tight correlation between [M/H] and [$\alpha$/Fe] (Sections 3 and 5).

More perplexing is the presence of a systematic trend between T$_\textrm{\footnotesize{eff}}$ and [$\alpha$/Fe], primarily for stars with [M/H]$\geq$$-$0.5 dex, that occurs while using both gravity-modeling  approaches, as seen from the middle panels in Figures 29 and 30. The higher-temperature stars tend to  have lower values of [$\alpha$/Fe] and lower-temperature stars have higher values of  [$\alpha$/Fe]. This correlation cannot have any physical basis because all stars are drawn from the same local population, and it is not possible that their chemical abundances be somehow correlated with their effective temperature. This resulting trend is independent of the initial values, $\chi$$^\textrm{\footnotesize{2}}$ expression (i.e., whether weighted with flux errors or not), and the extent of the searching grid in the model-fit pipeline. The range of values over which there is a potential degeneracy in the synthetic spectra while T$_\textrm{\footnotesize{eff}}$  and [$\alpha$/Fe] are varied (see Section 3) are not sufficiently extensive to cause a noticeable systematic trend in the diagram. Moreover, even if these ranges were large enough, the resulting trend would be in the opposite direction, i.e., higher-temperature stars would shift toward higher values of [$\alpha$/Fe], and inversely,  lower-temperature stars would shift towards lower values of  [$\alpha$/Fe]. Since  this trend is independent on the gravity-modeling approach, the correlation between surface gravity and   either T$_\textrm{\footnotesize{eff}}$ or [$\alpha$/Fe] does not have an important  role in the correlation between these two parameters. As a result, we believe that this trend is more likely due to some issues in the model atmospheres themselves (Section 3.4).

The bottom panels in Figures 29 and 30 show the same distribution as the upper panels, but now color-mapped according to the best-fit surface gravity of stars, for both gravity-modeling approaches. There is a clear trend between surface gravity and $\alpha$-element enhancement using the constant-gravity approach  (panel “e”) and between surface gravity and metallicity using the variable-gravity approach (panel “f”), in the  near-solar metallicity regime. The former trend can be attributed to the correlation between surface gravity, that are equal to the photometric values,  and effective temperature, as shown in Figure 31.   The surface gravity values seem to follow quite the same trend in [$\alpha$/Fe] as effective temperatures do (panel “c”), which   indicates the tight correlation between   log g and T$_\textrm{\footnotesize{eff}}$;  stars with   lower values of  log g have higher values of T$_\textrm{\footnotesize{eff}}$  (as predicted by theoretical isochrones, Figures 11 and 31), and accordingly,  lower values of [$\alpha$/Fe].  However, the latter trend is  most likely due to systematic errors caused by the strong degeneracy  between the two parameters log g and [M/H]. When gravity is allowed to vary during model-fitting, a significant number of stars  shift towards high values of [M/H]  and  log g. In addition,  some metal-rich  stars move toward the region with higher values of [$\alpha$/Fe], which is mostly due to  the degeneracy between log g and [$\alpha$/Fe]. This shift is clearly inconsistent with previous analyses of FGK-dwarf distributions  in the [$\alpha$/Fe]  versus [M/H] diagram (e.g., Lee et al. 2011b; Adibekyan et al. 2012, 2013; Recio-Blanco et al. 2014); metal-rich dwarfs, that belong to the Galactic disk and are among the young stellar population, are expected to have low values of  [$\alpha$/Fe], as compared to the old halo stars that  are believed to be rich in $\alpha$ elements.

We find that our model-fit values from the pipeline in which surface gravity is kept fixed  are generally  in  good agreement with those of previous analyses. More importantly, the  distribution of  the chemical parameter estimates  inferred from the constant-gravity approach in the abundance diagram of [$\alpha$/Fe]  versus [M/H]  is  comparable with that from  the latest (third) release of the Galactic Archaeology with HERMES (GALAH) survey  (Buder et al.  2021, Figure 5). The GALAH DR3 provides accurate atmospheric parameters and individual chemical abundances of 30 elements (11 of which are based on NLTE computations) from high-resolution optical spectroscopy for  mostly nearby  $\sim$500000 stars[The sample includes 65% dwarfs, 34% giants, and 1% as unclassified objects.] in the Galactic disk. This indicates that our best-fit  estimates  are sufficiently precise to show trends compatible with those from other studies using much more extensive samples.

As shown in the top panel of Figure 31, the parameter values derived from the constant-gravity approach can show a clear contrast in log g between metal-rich and metal-poor stars.  The metal-poor stars have larger values of log g, as compared to the metal-rich stars with the same T$_\textrm{\footnotesize{eff}}$. For a given effective temperature, low-metallicity  M subdwarfs have been shown to be as much as five times smaller than their solar-metallicity counterparts (Kesseli et al. 2019). Consequently,  the surface gravity of these metal-poor dwarfs is expected to be larger than that of metal-rich stars, for a specific  temperature. However, we find neither  a trend in log g between metal-poor and metal-rich stars nor a relationship between log g and T$_\textrm{\footnotesize{eff}}$ using parameter values obtained from the variable-gravity approach, as these values have large systematic uncertainties. The bottom panel of Figure 31 shows the same plot as presented in the top panel, but for near-solar metallicity stars ([M/H]$\geq$$-$0.5 dex) with a more detailed metallicity stratification. Overplotted are the solar-metallicity isochrones associated with a stellar age of 10 Gyr  and three different values of [$\alpha$/Fe]=+0.0, +0.2, and +0.4 dex from the Dartmouth stellar evolutionary models. The isochrones with higher values of [$\alpha$/Fe], i.e., +0.2 and +0.4  seem to be in better agreement with the model-fit metallicities, as compared to the isochrone with [$\alpha$/Fe]=+0.0 dex[Since the number of metal-poor stars in our sample is not statistically large, the comparison between the isochrones  and the distribution of stars with specific metallicity values in low-metallicity regime is not appropriate.]. However, since the photometric surface gravities, model-fit parameter values, and isochrones may not be completely accurate, no clear conclusion can be made for the location of the isochrones in the diagram. Nevertheless, the values of T$_\textrm{\footnotesize{eff}}$  show a decreasing trend with increasing log g, which  is consistent with that of the  theoretical  isochrones, regardless of their location relative to the stars' distribution.

Figure 32 compares HR diagrams of the stars in our subset, color-mapped according to the values of surface gravities derived from the constant-gravity approach (or photometric values, 2877 stars, top panel) and from the variable-gravity approach (2851 stars, bottom panel). The variation of photometric surface gravities in the color-absolute magnitude diagram is consistent with the relationship between these values and effective temperatures, as predicted by theoretical evolutionary models (Figures 11 and 31). Furthermore,  metal-poor M subdwarfs tend to have higher values of  surface gravity, which is expected from their smaller size, as compared to their metal-rich counterparts. However, the values of log g obtained from the variable-gravity method, that are affected by large systematic uncertainties, show no acceptable trend in  log g.

§ RESULTS: REDUCED-CORRELATION METHOD

To see whether a different fitting method may change the results, we apply the “reduced-correlation method” to the 3745 M dwarfs/subdwarfs  again for two gravity-modeling approaches, i.e., when surface gravity is kept constant and when surface gravity is allowed to vary. We analyze  the distribution of the  best-fit estimates using  some key diagrams and investigate   the influence of surface gravity on the values of the  other parameters, similar to our examination of resulting values from  the “normal method” (Section 6). We also compare the  results  from the two methods and discuss how the reducing of the residual correlation between wavelength datapoints can affect the derived parameter values.

§.§ Effective Temperature Distribution

Figure 33 depicts the same comparison as shown in Figure 19, but using the inferred temperature estimates from the “reduced-correlation method”.  The scatter  around the 1:1 lines (black) shown in each panel is very similar to the one presented in the corresponding panel in Figure 19. This indicates that the “reduced-correlation method” does not make  a major difference  in the best-fit temperature values relative to those obtained from the “normal method”.

§.§ Gaia HR Diagram

The top  and the bottom panels of Figure 34 are the same as the top panels of Figures 22 and 23, respectively, but color-mapped using the temperature estimates from the “reduced-correlation method”. Similarly, the top  and the bottom panels of Figure 35 are the same as the top panels of Figures 24 and 25, respectively, but color-coded based on  the metallicity estimates from the “reduced-correlation method”. The RUWE values are not taken into consideration in Figures 34 and 35. Despite some slight differences, the variation of temperatures and metallicities in these color-absolute magnitude diagrams are much the same as what can be seen from the corresponding plots in the top panels of  Figures 22-25. The stratification in T$_\textrm{\footnotesize{eff}}$ and [M/H] is still prominent, which again suggests that our model-fit method is relatively reliable. However, there  seems to be again  the same heavy concentration of high-temperature  stars (T$_\textrm{\footnotesize{eff}}$$\gtrsim$3600 K) in the narrow metallicity range $-$0.5$\lesssim$[M/H]$\lesssim$+0.1 dex, which is likely  owing to issues in the models of M dwarf atmospheres.

We do not include  any further plots that are  equivalent to those of Figures 26-32 when using the  best-fit values derived from the “reduced-correlation method”  because such plots are all remarkably similar to the corresponding ones using the  best-fit values inferred from the “normal method”, and  consequently,  would provide no more information about the  distributions of resulting parameter values.

§.§ Normal versus Reduced-Correlation Method

Figure 36 compares the best-fit estimates of T$_\textrm{\footnotesize{eff}}$, [M/H], and [$\alpha$/Fe] obtained from the normal and “reduced-correlation method”  for stars with T$_\textrm{\footnotesize{eff}}$$\leq$3550 K using the constant gravity approach (top panels).  All parameter values are randomized within the corresponding bin size of our model grid. These plots are color-mapped based on the constant (photometric) surface  gravity values.  The bottom panels show the histograms of the differences between the values of respective parameters derived from the two methods. Regardless of very slight offsets, the scatter of stars around the 1:1 lines is almost symmetric. The histogram panels indicate nearly normal distributions with an average around zero value.  We calculate the standard deviation of these differences: $\sim$44 K  in T$_\textrm{\footnotesize{eff}}$,  0.08 dex in [M/H], and 0.04 dex in [$\alpha$/Fe].

Figure 37 shows the same comparison as presented in Figure 36, but using the variable gravity approach. The plots are color-mapped based on the surface gravity values inferred from the “normal method”. We find a little larger shift from the 1:1 line for effective temperature (panel “a”) and metallicity (panel “b”), as compared to the equivalent panels in Figure 36. On the other hand, there is no noticeable offset relative to the 1:1 line for $\alpha$-element enhancement (panel “c”). These are also clear from the histograms in the bottom panels; the distributions show a $\sim$20-30 K offset from the zero value in the panel “d” and a $\sim$0.02-0.03 dex  offset from the zero value in the panel “e” towards higher values, while there is no  shift in the panel “f”. We determine the standard deviation of these differences around the average  values: $\sim$46 K  in T$_\textrm{\footnotesize{eff}}$, 0.1 dex in [M/H], and and 0.04 dex in [$\alpha$/Fe].

In general, there is a deviation within  $\sim$45 K in T$_\textrm{\footnotesize{eff}}$, within $\sim$0.1 dex in [M/H], and within $\sim$0.04 dex in [$\alpha$/Fe] between the values obtained from the normal and “reduced-correlation method”. As a result, the “reduced-correlation method”, which decreases the correlation of the residuals between spectral datapoints due to instrumental limitations, does not make a significant difference in the model-fit parameter values and their resulting distributions, which is likely due to the low spectral resolution of our spectra. The uncertainties caused by these correlations are thus unlikely  to be of importance in the fitting of synthetic spectra to our observed spectra (or in general, to low-resolution M dwarf spectra). This suggests that systematic errors in the model-fit method are unlikely to be due to the fitting technique itself, but more likely due to issues in the models and/or degeneracies in the synthetic spectra discussed in previous sections.

[Fixed surface gravity]

[Variable surface gravity]
Comparison between the photometric temperatures and the inferred model-fit values for the 3745 stars, when surface gravity is considered as a fixed parameter  (top panel) and when surface gravity is considered as a free parameter (bottom panel), using the reduced-correlation method.

[Fixed surface gravity]

[Variable surface gravity]

[Fixed surface gravity]

[Variable surface gravity]
HR diagram,  M$_\textrm{\footnotesize{G}}$ versus G$_\textrm{\footnotesize{BP}}$ - G$_\textrm{\footnotesize{RP}}$ (top panels) and  M$_\textrm{\footnotesize{G}}$ versus G - G$_\textrm{\footnotesize{RP}}$ (bottom panels), of the 3745 stars,  when surface gravity is set a fixed parameter (left panels) and when surface gravity is set a free parameter (right panels), color-mapped based on  effective temperatures, using the reduced-correlation method. The values of RUWE are not taken into account in this figure.

[Fixed surface gravity]

[Variable surface gravity]

[Fixed surface gravity]

[Variable surface gravity]
HR diagram,  M$_\textrm{\footnotesize{G}}$ versus G$_\textrm{\footnotesize{BP}}$ - G$_\textrm{\footnotesize{RP}}$ (top panels) and  M$_\textrm{\footnotesize{G}}$ versus G - G$_\textrm{\footnotesize{RP}}$ (bottom panels), of the 3745 stars,  when surface gravity is set a fixed parameter (left panels) and when surface gravity is set a free parameter (right panels), color-coded based on five groups with different  metallicity ranges, using the reduced-correlation method. The number of stars in each metallicity group, N, is panels “a” and “c”:  N$_\textrm{{Black}}$=78,  N$_\textrm{{Yellow}}$=2822, N$_\textrm{{Blue}}$=644, N$_\textrm{{Red}}$=98, and N$_\textrm{{Green}}$=103, panel “b” and “d”:  N$_\textrm{{Black}}$=757,  N$_\textrm{{Yellow}}$=2294, N$_\textrm{{Blue}}$=494, N$_\textrm{{Red}}$=108, and N$_\textrm{{Green}}$=92. The values of RUWE are not taken into account in this figure.

Comparison between inferred model-fit parameters, T$_\textrm{\footnotesize{eff}}$, [M/H], and  [$\alpha$/Fe] using the reduced-correlation method and corresponding values using the normal method for 2836 stars with T$_\textrm{\footnotesize{eff}}$$\leq$3550 K (top panels). Surface gravity is treated as a constant parameter in both methods. The bottom panels show the histograms of the differences between  corresponding values inferred from the two methods.

Comparison between inferred model-fit chemical  parameters, T$_\textrm{\footnotesize{eff}}$, [M/H], and  [$\alpha$/Fe] using the reduced-correlation method and corresponding values using the normal method for 2788 stars with T$_\textrm{\footnotesize{eff}}$$\leq$3550 K (top panels). Surface gravity is treated as a free parameter in both methods. The bottom panels show the histograms of the differences between corresponding values inferred from the two methods.

[Fixed surface gravity]

[Fixed surface gravity]

[Fixed surface gravity]

[Variable surface gravity]

[Variable surface gravity]

[Variable surface gravity]
Comparison between the inferred model-fit values of  [M/H], [$\alpha$/Fe], and the combined parameter [M/H]+[$\alpha$/Fe] of the primaries and  those of their respective companions, for binary systems whose both components have T$_\textrm{\footnotesize{eff}}$$\leq$3550 K, when surface gravity is a constant parameter (top panels, 74 binaries) and when surface gravity is a free parameter (bottom panels, 71 binaries) in the model fitting process, using the normal method.

[Fixed surface gravity]

[Fixed surface gravity]

[Fixed surface gravity]

[Variable surface gravity]

[Variable surface gravity]

[Variable surface gravity]
Comparison between the inferred model-fit values of  [M/H], [$\alpha$/Fe], and the combined parameter [M/H]+[$\alpha$/Fe] of the primaries and those of their respective companions, for binary systems whose both components have T$_\textrm{\footnotesize{eff}}$$\leq$3550 K, when surface gravity is a constant parameter (top panels, 68 binaries) and when surface gravity is a free parameter (bottom panels, 68 binaries) in the model fitting process, using the reduced-correlation method.

Kinematic distributions of the stars in three different directions of the sky:  the Galactic poles (Group UV, left panels),  the Galactic center and anti-center (Group VW, middle panels), and  the solar apex and anti-apex (Group UW, right panels).   These panels show the projected motion in the Galactic UV, VW, and UW planes for a subset of stars having  T$_\textrm{\footnotesize{eff}}$$\leq$3550 K and  $-$2.5<[M/H]<+0.5 dex, that are color-mapped according to their metallicities,  inferred  from the  normal method when surface gravity is regarded as a fixed parameter: 1013 stars (Group UV), 803 stars (Group VW), and 1051 stars (Group UW).

[$\alpha$/Fe] vs. [M/H] diagram for the stars with T$_\textrm{\footnotesize{eff}}$$\leq$3550 K and $-$2.5<[M/H]<+0.5 dex, in Group UV+UW (top panels), Group UV+VW (middle panels), and Group UW+VW (bottom panels). The parameter values are inferred from  the normal method using the constant-gravity approach: 2065 stars in Group UV+UW, 1817 stars in Group UV+VW, and 1849 stars in  Group UW+VW. The stars are color-mapped based on U, V, and W velocity components for  Groups UV+UW, UV+VW, and UW+VW, respectively.

Velocity component  U (left panels),  V (middle panels), and W (right panels) versus [M/H] (top panels) and [M/H]+[$\alpha$/Fe] (bottom panels) of the stars  with T$_\textrm{\footnotesize{eff}}$$\leq$3550 K and $-$2.5<[M/H]<+0.5 dex in the respective combined Group UV+UW,  UV+VW, and UW+VW, as  those shown in Figure 41, but with narrower ranges of velocities mainly associated with the disk stars:  $-$200$\leq$U$\leq$200 km/s (2028 stars in Group UV+UW), $-$40$\leq$V$\leq$260 km/s (1747 stars in Group UV+VW) , and $-$100$\leq$W$\leq$100 km/s (1780 stars in Group UW+VW). The parameter values are inferred from the  normal method in which surface gravity is  a fixed parameter.

Velocity component V versus [M/H]+[$\alpha$/Fe] of the stars  with T$_\textrm{\footnotesize{eff}}$$\leq$3550 K and $-$2.5<[M/H]<+0.5 dex in the combined Group UV+VW, as shown in the bottom middle panel of Figure 42.  Three stellar streams are marked in red.

§ PRECISION ANALYSIS FROM COMMON PROPER-MOTION PAIRS

Wide binaries are believed to be formed in different processes (see e.g., Hawkins et al. 2019 for more details). Pairs with smaller  separations, i.e., from  a hundred to a few thousand AU, are more likely to form during the turbulent  fragmentation of a molecular core (e.g. Offner et al. 2010; Lee et al. 2017). Wider pairs with separations between 0.01 to 1 pc  are thought to be born  by the dynamical scattering of  one component in  triple systems (Reipurth & Mikkola 2012), the dissipation of stellar clusters (e.g. Kouwenhoven et al. 2010; Moeckel & Clarke 2011), or the pairing of  adjacent pre-stellar  cores  that move slowly enough relative to each other (e.g. Tokovinin 2017).  In most cases, it can be assumed that the two components in wide binary systems are formed from the same molecular cloud and have nearly the same chemical composition, and this is generally borne out by most observations of co-moving pairs. Occasional observed discrepancies in the elemental abundances of the components in some co-moving pairs have been attributed to  the general dynamical coldness of the thin disk or resonances induced by the Galactic tidal field (e.g., Simpson et al. 2019). It has been argued that these pairs of stars are “coincidentally” on similar orbits because of such dynamical interactions. Detailed investigations of chemical abundance distributions of stellar groups or pairs  are of great importance to understand the circumstances under which the members exhibit significant abundance differences and to find out whether these stars truly originate from a single parent cloud or have coincidentally converged to a similar orbit.