{'fraction_non_alphanumeric': 0.032124326293704535, 'fraction_numerical': 0.0246366069072955, 'mean_word_length': 4.524186272468596, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Near-infrared Metallicities, Radial Velocities and Spectral Types for 447 Nearby M Dwarfs

We present metallicities, radial velocities and near-infrared spectral types for 447 M dwarfs determined from moderate resolution (R~2000) near-infrared (NIR) spectra obtained with IRTF/SpeX. These M dwarfs are targets of the MEarth Survey, a transiting planet survey searching for super Earths around mid-to-late M dwarfs within 33pc. We present NIR spectral types and new IRTF spectral templates in the Z, J, H and K-bands, created using M dwarfs with near-solar metallicities. We developed two spectroscopic distance calibrations that use NIR spectral type or an index based on the curvature of the K-band continuum. Our distance calibration has a scatter of 14%. We searched 27 NIR spectral lines and 10 spectral indices for metallicity sensitive features, taking into account correlated noise in our estimates of the errors on these parameters. We calibrated our relation using 36 M dwarfs in common proper pairs with an F, G or K-type star of known metallicity. We validated the physical association of these pairs using proper motions, radial velocities and spectroscopic distance estimates. Our resulting metallicity calibration uses the sodium doublet at 2.2um as the sole indicator for metallicity. It has an accuracy of 0.12dex inferred from the scatter between the metallicities of the primaries and the estimated metallicities of the secondaries. Our relation is valid for NIR spectral types from M1V to M5V and for -1.0<[Fe/H]<+0.35dex. We present a new metallicity relation using J-H and J-K colors that directly relates two observables: the distance from the M dwarf main sequence and equivalent width of the sodium line at 2.2um. We measured radial velocities by modeling telluric features to determine the absolute wavelength calibration of our spectra, and used M dwarf binaries, observations at different epochs, and comparison to precisely measured radial velocities to demonstrate 4.4km/s accuracy.


INTRODUCTION
MEarth is a transiting planet survey looking for super Earths around nearby mid to late M dwarfs. As part of our efforts to characterize the local M dwarf population, the MEarth team and collaborators are gathering a diverse data set on these low mass stars. These unique data have already begun to bear fruit. Charbonneau et al. (2009) reported the discovery of a super Earth transiting the mid M dwarf GJ 1214. Irwin et al. (2011a) took advantage of our long-baseline photometry to measure rotation periods as long as 120 days for 41 M dwarfs and investigated their angular momentum evolution, finding that strong winds may be needed to explain the population of slowly rotating field M dwarfs. Irwin et al. (2011b) presented a long period M dwarf-M dwarf eclipsing binary and measured the masses of the two components and the sum of their radii. They find the radii to be inflated by 4% relative to theoretical predictions, reflecting a well-known problem with stellar models at the bottom of the main sequence (e.g. Lopez-Morales 2007; Boyajian et al. 2012).
Interest in M dwarfs is fueled by prospects for testing theories of planet formation. Creating a planetary system around a small star is one of the simplest ways to test the effect of initial conditions: the disk out of which planets form is less massive around an M dwarf than around a more massive star. Core ac-cretion and gravitational instability models predict different rates of occurrence of planets around low-mass stars, with the formation of giant planets through core accretion being hampered by the low disk surface density and long orbital time scale in M dwarf protoplanetary disks (Laughlin et al. 2004). Recent results from Kepler showed that giant planets are less likely to be found around K and early M stars than around F and G stars, lending support to the core accretion model (Borucki et al. 2011;Fressin et al. 2013). A similar finding was reported for M dwarfs targeted by radial velocity surveys (Johnson et al. 2007;Cumming et al. 2008). The high metallicity of solar-type stars that host close-in giant planets was confirmed over a decade (e.g. , but smaller planets have been found around stars of a range of metallicities (Buchhave et al. 2012). Efforts have been made to extend these relations to the lowest stellar masses (e.g. Johnson & Apps 2009;Schlaufman & Laughlin 2010;Rojas-Ayala et al. 2012), but have been limited by the small number of planets currently known around M dwarfs.
M dwarfs present a unique opportunity for the detection and characterization of habitable Earth-sized planets. Mid to late M dwarfs are favorable targets for transiting planet searches (Nutzman & Charbonneau 2008). Their low luminosity puts the habitable zone at smaller orbital radii, making transits more likely and more frequent: for an M4 dwarf, the period of a habitable planet is two weeks, compared to one year for a solar-type star. Because the transit depth is set by the planetto-star radius ratio, smaller planets are more readily detectable around these stars. The small radius of an M dwarf is also favorable for follow-up studies of an orbiting planet's atmosphere with transmission or occultation techniques and nearby mid M dwarfs are bright enough in the NIR for precise spec-troscopic studies (e.g. Bean et al. 2011;Crossfield et al. 2011;Berta et al. 2012).
In contrast to solar type stars, the physical parameters of M dwarfs are not in general well understood and present a major hurdle for studying transiting planets orbiting M dwarfs. Few M dwarfs are bright enough for direct measurement of their radii (e.g. Berger et al. 2006;Boyajian et al. 2012), and discrepancies between the observed radii and theoretical predictions persist (see Torres 2013, for a review). Empirical calibrations provide an inroad. For example, Muirhead et al. (2012a) and Muirhead et al. (2012b) exploited the K-band metallicity and temperature relations of Rojas-Ayala et al. (2012, hereafter R12) to estimate new planet properties for the Kepler Objects of Interest (KOIs) orbiting the coolest Kepler stars and discovered the planetary system with the smallest planets currently known, the Kepler-42 system (née KOI-961). Johnson et al. (2012) combined existing photometric relations to estimate the stellar properties of KOI-254, one of the few M dwarfs known to host a hot Jupiter. Ballard et al. (2013) used M dwarfs with interferometric radii as a proxy to constrain the radius and effective temperature of Kepler-61b.
Several studies have used M dwarf model atmospheres matched to high resolution spectra to determine stellar parameters. Woolf & Wallerstein (2005) estimated M dwarf temperatures and surface gravities from photometry, then, fixing these parameters, inferred the metallicity from the equivalent widths (EWs) of metal lines. Updating and modifying the spectral synthesis method of Valenti et al. (1998), Bean et al. (2006a) used TiO and atomic lines in combination with NextGen PHOENIX model atmospheres (Hauschildt et al. 1999) to measure the physical properties of M dwarfs. Most recently,Önehag et al. (2012) matched model spectra from MARCS (Gustafsson et al. 2008) to observations of FeH molecular features in the infrared and found metallicities higher than those inferred by Bean et al. (2006b). The MARCS model atmospheres do not include dust formation and are not applicable to M dwarfs later than M6V. However, uncertain sources of opacity in the model atmospheres complicate direct interpretation of observed spectra throughout the M spectral class.
An effective technique for quantitatively studying the metallicities of M dwarfs makes use of cool stars in common proper motion (CPM) pairs with an F, G or K-type star, where the primary has a measured metallicity. Assuming the two are coeval, one can infer the metallicity of the low-mass companion and subsequently use a sample of CPM pairs to confirm or empirically calibrate tracers of M dwarf metallicity. Gizis & Reid (1997) applied this idea to the M subdwarf population, using observations of late-type companions to F and G subdwarfs of known metallicity to confirm the metallicity relation of Gizis (1997), which used optical spectral indices to infer the metallicity of M subdwarfs. Bonfils et al. (2005) pioneered the empirical calibration of M dwarf metallicities using CPM pairs. The authors found that a metal-rich M dwarf has a redder V − K color at a given absolute K magnitude, due to increased line blanketing by molecular species, particularly TiO and VO. The calibration is valid for 4 < M K < 7.5, 2.5 < V − K < 6 and −1.5 < [Fe/H] < +0.2 dex. Bonfils et al. (2005) reported a standard deviation of 0.2 dex. Johnson & Apps (2009), finding the calibration of Bonfils et al. (2005) to systematically underestimate the metallicities of metal-rich stars, updated the relation by considering the offset from the mean main se-quence (MS), assuming the mean MS defined an isometallicity contour with [Fe/H] = −0.05 dex. Their calibration sample used six metal-rich calibrators. Schlaufman & Laughlin (2010) found that the previous works had systematic errors at low and high metallicities and further updated the photometric relation. They used a larger calibration sample comprised only of M dwarfs with precise V magnitudes in CPM pairs with an F, G or K-star, where the primary's metallicity had been determined from high resolution spectroscopy. They also updated the determination of the mean MS, finding that it corresponded to an isometallicity contour with [Fe/H] = −0.14 dex. However, external information was still required to determine the mean MS. The standard deviation of their fit was 0.15 dex. Neves et al. (2012) tested the photometric calibrations of Bonfils et al. (2005), Johnson & Apps (2009), and Schlaufman & Laughlin (2010) on a new sample of FGK-M CPM pairs that had precise V -band photometry. With their sample of 23 M dwarfs, they found the Schlaufman & Laughlin (2010) calibration had the lowest residual mean square error (RMSE = 0.19 ± 0.03 dex) and highest correlation coefficient (R 2 ap = 0.41 ± 0.29), performing marginally better than the Bonfils et al. (2005) calibration. They updated the Schlaufman & Laughlin (2010) calibration, though the diagnostic values did not improve by more than the associated errors.
Rojas-Ayala et al. (2010, hereafter R10) took a different approach and used moderate resolution K-band spectra (R ≈ ∆λ/λ ≈ 2700) to measure metallicity. They used the EWs of the Na I doublet and Ca I triplet to measure metallicity and the H 2 O-K2 index to account for the effects of temperature. The calibration was updated in Rojas-Ayala et al. (2012, hereafter R12), who demonstrated that their empirical metallicities gave reasonable results for solar neighborhood M dwarfs. With 18 calibrators, this method yielded RMSE = 0.14 dex and R 2 ap = 0.67 for their [Fe/H] calibration. The lines used in this calibration are isolated across the entire M dwarf spectral sequence and are located near the peak of the M dwarf spectral energy distribution (SED). Parallaxes and accurate magnitudes, which are scarce for M dwarfs, are not required, placing metallicities within reach for many M dwarfs. Terrien et al. (2012) applied the methods of R10 to spectra obtained with the SpeX instrument on the NASA Infrared Telescope Facility (IRTF), using 22 CPM pairs as calibrators. They updated the K-band R10 calibration (RMSE = 0.14 dex, R 2 ap = 0.74) and presented an H-band calibration (RMSE = 0.14 dex, R 2 ap = 0.73). Mann et al. (2013) expanded the sample of calibrators and identified over 100 metal-sensitive features in the NIR and optical. Their calibration sample included 112 FGK-M CPM pairs, selected on the basis of common proper motion and galactic models. They constructed metallicity relations in the optical and in each of the NIR bands out of metallicity sensitive features and a single parameter to account for temperature dependencies. Their [Fe/H] calibrations had standard deviations between 0.11 dex and 0.16 dex and R 2 ap values ranging from 0.68 to 0.86. They also updated the color-color relation of Johnson et al. (2012) and the Kand H-band spectroscopic relations of Terrien et al. (2012) and R12.
We also note the larger context in which constraints on the physical properties of M dwarfs are applicable. For example, Bochanski et al. (2007) used SDSS M dwarfs to test the Besançon galactic model (Robin et al. 2003), comparing ob-served kinematics to the model and comparing the observed metallicities and active fractions of the thin and thick disk. In this study and others using SDSS, optical molecular indices were used as a proxy for metallicity (e.g. Gizis & Reid 1997;Woolf & Wallerstein 2006). The ζ-index, which uses CaH and TiO molecular band heads, is commonly used to identify subdwarfs and extreme subdwarfs (Lepine et al. 2007;Dhital et al. 2012). Theories of star formation must also match the observed luminosity and mass functions of M dwarfs, which are in turn important input into galactic models. Bochanski et al. (2010), again exploiting SDSS, measured the M dwarf luminosity and mass functions. Photometric distance estimates were used in this work, and one of the primary factors complicating these estimates was uncertainty in how metallicity affects absolute magnitude.
In this work, we present our observation and analysis of near infrared (NIR) moderate resolution (R ≈ 2000) spectra of 447 MEarth M dwarfs. Our sample is presented in §2 and in §3 we discuss our observations and data reduction. We account for correlated noise when estimating the error on our measurements, as we discuss in §4. In §5, we present by-eye NIR spectral types for each star and a new spectroscopic distance calibration. Our metallicity measurements, described in §6, are based on the method developed by R12: we use EWs of spectral features in the NIR as empirical tracers of metallicity, using M dwarfs in CPM pairs to calibrate our relationship. We present a color-color metallicity calibration in §7. In §8, we discuss our method for measuring radial velocities, which uses telluric features to provide the wavelength calibration, and demonstrate 4 km s −1 accuracy. Our data are presented in Table A1 and we include updated parameters for those stars observed by R12 in Table A2. We include radial velocities, spectral types and parallaxes compiled from the literature.

SAMPLE
Our sample consists of 447 M dwarfs targeted by the MEarth transiting planet survey and 46 M dwarfs in CPM pairs with an F, G or K star of known metallicity, a subset of which we used to calibrate our empirical metallicity relation.

MEarth M dwarfs
The MEarth project is photometrically monitoring 2000 of the nearest mid to late M dwarfs in the northern sky with the goal of finding transiting super Earths. Nutzman & Charbonneau (2008) described how the MEarth targets were selected from the Lépine-Shara Proper Motion catalog of northern stars (LSPM-North; Lépine & Shara 2005). For completeness, we summarize their method here. From the subset of stars believed to be within 33 pc (Lépine 2005), using spectroscopic or photometric distance estimates where parallaxes were unavailable, they selected those with V − J > 2.3, J − K S > 0.7, and J − H > 0.15, resulting in a sample of probable nearby M dwarfs. The radius for each probable M dwarf was estimated by first using the absolute K S magnitude-to-mass relation of Delfosse et al. (2000), and inputting this mass into the mass-to-radius relationship from Bayless & Orosz (2006). They subsequently selected all objects with estimated radii below 0.33R ⊙ , driven by the desire to maintain sensitivity to planets with radii equal to twice Earth's.
MEarth is a targeted survey, visiting each object with a cadence of 20-30 minutes on each night over one or more observing seasons. A fraction of the sample has sufficient coverage and quality to estimate their rotation periods, with re-covered periods ranging from 0.1 to 90 days. These will be discussed in a subsequent paper.

Spectroscopy targets
We targeted a subset of the MEarth M dwarfs for NIR spectroscopy. We re-observed 30 stars in common with R12, who focused their efforts on M dwarfs within 8pc, in order to evaluate any systematic differences between our instruments and methods. The IRTF declination limit prevented us from observing stars above +70 • . We divide our targets into four subsamples based on the reason for their selection: • Rotation sample: 181 M dwarfs with preliminary rotation periods measured from MEarth photometry. These show periodic photometric modulation presumed to be due to star spots rotating in and out of view.
• Nearby sample: 257 M dwarfs drawn from the full MEarth sample, for which no clear periodic photometric modulation was detected at the time of selection. This included 131 M dwarfs selected because they have parallaxes available from the literature, 94 M dwarfs with photometric distance estimates, and 32 "photometrically quiet" M dwarfs. The photometrically quiet M dwarfs are those for which phase coverage and photometric noise were sufficient to achieve good sensitivity to rotationally induced photometric modulations, but for which no such modulations were observed.
• Metallicity calibrators: 46 M dwarfs in CPM pairs with an F, G or K, where a metallicity measurement is available for the primary. These are discussed in §6. We used 36 M dwarfs in our final metallicity calibration.
• Potential calibrators: 10 potential calibrators are in CPM pairs with an F, G or K star but do not have a metallicity measurement available for the primary. We did not include these stars in our metallicity calibration.
We present new observations of 447 nearby M dwarfs in Table A1 (the rotation and nearby samples and potential calibrators). Data for our 46 M dwarf metallicity calibrators are presented separately.

OBSERVATIONS
We conducted our observations with the SpeX instrument on the NASA Infrared Telescope Facility (Rayner et al. 2003, IRTF). We used the short cross dispersed (SXD) mode with the 0.3 × 15 ′′ slit. This yielded spectra with R ≈ 2000 covering 0.8 − 2.4µm, with gaps between orders where there is strong atmospheric absorption. Our observations spanned 25 partial nights over 4 semesters. Observing conditions are summarized in Table 1; in moderate clouds, we observed bright targets.
We typically acquired four observations of each object, with two observations at each of two nod positions (A and B), in the sequence ABBA. We used the default A position and nod distance, with the A and B positions falling 3. ′′ 75 from the edge of the slit (a 7. ′′ 5 separation). Most of our targets were observed within half an hour of meridian crossing. For hour angles greater than one, we aligned the slit with the parallactic angle. We observed A0V stars for use as telluric standards within one hour of each science target, at angular separations no more than 15 • , and with airmass differences of no more than 0.1 when possible (see §4). We took flat field spectra (using an internal quartz lamp) and wavelength calibrations (using internal Thorium-Argon lamps) throughout the night, at one hour intervals or after large slews. The typical observation time for a K = 9 target at each nod was 100 seconds (for a total integration time of 400 seconds). Combining four nods yielded a total signal-to-noise ratio (S/N) of 250 per resolution element. We reduced the data with the instrument-specific pipeline Spextool (Cushing et al. 2004), modified to allow greater automation and to use higher S/N flat fields, created by median combining all flat field frames from a given night. Images were first flat-field corrected using the master flat from the given night. After subtracting the A and B images, we used boxcar extraction with an aperture radius equal to the full width at half maximum (FWHM) of the average spatial profile and subtracted the residual sky background. To determine the background sky level in the AB subtracted image, we used a linear fit to the regions beginning 1. ′′ 2 from the edges of the aperture. This step was important near sunrise and sunset and increasingly important in bluer orders, but the K-band was largely unaffected. Each spectrum was wavelength calibrated using the set of Thorium-Argon exposure most closely matching in time.
We combined individual spectra for the same object (typically 4 per object) using the Spextool routine xcombspec. We scaled the raw spectra to the median flux level within a fixed wavelength region and removed low order variations in the spectral shapes. We used the highest S/N region of the H-band for scaling. The modified spectra were combined using the robust weighted mean algorithm, which removed outliers beyond 8σ.
We used xtellcor to perform the telluric corrections (Vacca et al. 2003). We used the Paschen δ line near 1µm in the A0V telluric standard to create a function to describe the instrumental profile and the rotational broadening observed in spectrum. We used xtellcor to convolve this function with a model of Vega and shifted the model to match the star's observed radial velocity. We scaled the line strengths of individual lines to match those observed; for data taken in 2012, we adjusted the scaling by hand. We found this to be a necessary step because even for sub-1% matches to the Vega model, residual hydrogen lines were apparent. The atmospheric absorption spectrum, as observed by the instrument, was found by dividing the observed A0V spectrum by the modified Vega spectrum. We shifted the atmospheric absorption spectrum to match the absorption features in the object spectrum and divided to remove the atmospheric absorption features present. We performed this step separately in each order, using a region dominated by telluric features to shift the spectra.
We performed flux calibration as part of the telluric correction, but variable weather conditions and slit losses made the absolute flux level unreliable. We do not require absolute flux calibration for our project goals.

ESTIMATION OF UNCERTAINTY
Given the high S/N (typically > 200) of our spectra, the uncertainties in quantities measured from our data are dominated by correlated noise, rather than random photoncounting errors. Correlated noise could be introduced by poorly-corrected telluric lines or by unresolved features in the region of the spectra assumed to represent the continuum.
We drew our errors from a multivariate Gaussian with Gaussian weights along the diagonal of the covariance matrix. At each pixel, we simulated Gaussian random noise using the errors returned by the SpeX pipeline, which included photon, residual sky, and read noise and which were propagated through the Spextool pipeline. We multiplied the error realization by a Gaussian centered on that pixel with unit area and full width at half maximum (FWHM) equal to the width of the autocorrelation function. To determine the appropriate FWHM, we autocorrelated each order of several spectra of different S/N and found that a Gaussian with a FWHM of 1.5 pixels approximated the width of the autocorrelation function; we used this FWHM for all stars. We did this for each pixel, resulting in an array of overlapping Gaussians of unit area, one centered on each pixel. We then added the contributions from the Gaussians at each pixel, and took the sum at each pixel to be the error on that pixel. This effectively spread the error associated with a single pixel over the neighboring pixels according to the autocorrelation function.
We then re-measured spectral indices (described below), EWs (described in §6.3) and the radial velocity (as described in §8.1). We repeated this process 50 times and calculated the 1σ confidence intervals, which we took to be the errors on our measurements.
To assess the accuracy of our error estimates, we considered stars that we observed on two separate occasions, which have different observing conditions and S/N. By comparing independent measurements of the same object, we determined whether our error estimates accurately model the observed differences in the measurements. We used EWs, which we measure by numerically integrating within a defined region, as indicators of M dwarf metallicity (our method is described in detail in §6). The line of most interest to us is the Na I line at 2.2µm. The median error on EW Na is 0.17Å, typically 5%, which was achieved with S/N = 300. 92% of our spectra have S/N in the K-band greater than 200 and 67% have errors on EW Na less than 0.2Å. In Figure 1, we compare EW Na for stars that were observed multiple times, finding that our method accurately captures the observed errors. -We compare EW Na measurements for stars for which we have more than one observation. The horizontal axis shows the EW Na of the selected observation and the vertical axis shows the EW Na of the alternate observation, both inÅ. We also include the 1σ confidence intervals from 50 trials.
We also measure 10 spectral indices ( §6.3), including the H 2 O-K2 index, a temperature-sensitive index that measures the curvature of the K-band by considering the flux level in three K-band regions (R12). It is defined as: Angle brackets represent the median flux within the wavelength range indicated, where wavelengths are given in microns. In Figure 2, we compare measurements of the H 2 O-K2 index for objects which were observed multiple times. Our autocorrelation analysis underestimated the true uncertainties. The largest discrepancies arose when airmass differed by more than 0.2 or time of observation differed by more than two hours (these were not typical occurrences amongst our sample). If using the H 2 O-K2 index for metallicity or temperature measurements, we suggest taking particular care to observe a telluric standard immediately before or after each science observation, and as closely matching in airmass as possible, as described in Vacca et al. (2003).

NIR SPECTRAL TYPES
We determined NIR spectral types by eye for each star using the K, H, J and Y -bands. Our NIR spectral types are based on the spectral typing system defined by Kirkpatrick et al. (1991Kirkpatrick et al. ( , 1995Kirkpatrick et al. ( , 1999, hereafter the KHM system. We used a custom spectral typing program to match each science spectrum to a library of spectral type standards created from our data ( §5.1- §5.2). We considered the differences between our NIR spectral types and other spectral type indicators ( §5.3) and calibrated a new spectroscopic distance relation using apparent K S magnitude and either NIR spectral type or the H 2 O-K2 index ( §5.4).

Spectral typing routine
We first estimated the spectral type for each star using the relationship between H 2 O-K2 index and spectral type that was presented in R12. We displayed the object spectrum and two spectral standards: the spectral standard with the estimated spectral type and the spectral standard with the spectral type one subtype later. We indicated the FeH bands identified in Cushing et al. (2005)  index of the selected observation and on the vertical axis, the H 2 O-K2 index of the alternate observation. The errors from 50 trials are smaller than the data points. We indicate the cases of significant airmass discrepancies between the science and telluric spectra as triangles (for ∆AM > 0.2) and diamonds (for 0.2 > ∆AM > 0.1). The two cases with large discrepancies in the H 2 O-K2 index but for which the science and telluric spectra are closely matching in airmass are instances where the science and telluric observations were separated by more than two hours.
FIG. 3.-An example of the output from our spectral typing routine. We included the K, H, J and Y -bands in our program. We show the object spectrum, in this case GJ 1214, in black. We overplot two spectral standards in blue and red. Dashes indicate FeH bands; only the Wing-Ford band head at 0.99µm is apparent in mid M dwarfs. In this case, we selected the blue spectral standard, M4V, as the best match to the object spectrum. The spectral type from Reid et al. (1995) is M4.5V.
Wing-Ford FeH band at 0.99µm is the only band head apparent across the entire M spectral sequence. FeH is known to be sensitive to spectral type (e.g. Schiavon et al. 1997;Cushing et al. 2005). Using a GUI, we checked earlier and later spectral standards as desired, then selected a spectral type for the object. An example is shown in Figure 3. For K7V and M0V, we used the spectral standards from the IRTF library. For the remaining spectral types, we created standards from our observations by median-combining stars of a single spectral type. We were unable to reliably separate M8V and M9V stars and therefore treat them as one spectral category (see §5.1). In practice, we also could not distinguish between K7V and M0V and assigned these a K7/M0V spectral type. We did not consider half-spectral types. We found the difference between late K dwarfs and M0V stars, and similarly between M8V and M9V stars, to be marginal in the NIR. We used a combined M8V/M9V spectral standard in our program. While K7V and M0V spectral standards were included separately in our spectral typing code, in our later analysis stars we considered a joint K7/M0V spectral class. We took a holistic approach to spectral typing due to the metallicity-dependence of many spectral features. We placed more weight on the redder orders and less weight on features known to be sensitive to metal content (such as the sodium line at 2.2µm). Our NIR spectral types are included in Table A1.

IRTF spectral standards
We initially used the M dwarfs in the IRTF spectral library (Rayner et al. 2009) as spectral standards, using the KHM spectral standards except for our M0V (HD19305), M3V (AD Leo/Gl 388) and M6V (CN Leo/Gl 406) spectral standards. However, we noted several differences between the strengths of features in the standard spectra and the typical object spectra. In particular, the M4V spectral standard, Gl 213, is metal poor. This is to be expected: Cushing et al. (2005) identify Gl 213 as a probable low-metallicity object on the basis of its low Fe, Al, Na, and Ca EWs. By comparison with neighboring spectral standards and using our holistic approach to spectral typing, we were nevertheless able to accurately assess the NIR spectral types of solar metallicity stars.
To address the concern of spectral standards with extreme metallicities or other unique features, we created our own standard spectra. We assessed the spectral type of all stars observed through the 2012A semester by eye once, using the IRTF spectral library stars as standards. We then mediancombine stars of a single spectral type that were within 0.2 dex of solar metallicity or, for M5V-M7V stars, within 0.1 dex of solar (see §6 for a description of how we determine metallicities for our stars). There were two stars comprising the M1V spectral standard (with median We included all five M8/9V stars observed through the 2012A semester in the M8/M9V spectral standard. We continued to use the IRTF spectral library standards for K dwarfs and M0V stars. We show our spectral sequence in four IRTF bands, from K7V to M8/9V, in Figs. 4-7. We then re-classified each star by eye using our new standard spectra.

Comparing measures of spectral type
We first compare our by-eye NIR spectral types to those measured with the H 2 O-K2 index, using the relation in R12. These measures agree to within one spectral type; however, our by-eye spectral types are on average half a spectral type later than those measured using the H 2 O-K2 index. We express M subtype numerically as Sp NIR , where positive values are M subtypes (e.g. Sp NIR = 4 corresponds to M4V) and negative values are K subtypes (e.g. Sp NIR = −1 corresponds to K7V and Sp NIR = −2 corresponds to K5V). We find: Over 100 of our objects have optical spectral types from the Palomar/Michigan State University (PMSU) Survey (Reid et al. 1995;Hawley et al. 1996, included for comparison in Table A1). The PMSU survey used the depth of the strongest TiO feature in optical M dwarf spectra as the primary indicator of spectral type, and calibrated their relation against nearly 100 spectral classifications on the KHM system. As in R12, we find a systematic difference between the PMSU spectral types and the NIR spectral types as a function of metallicity, shown in Figure 8 for stars earlier than M5V. For M5V stars, there appears to be no clear trend with metallicity.
For early and mid M dwarfs, the NIR spectral type is typically half a spectral type later than the PMSU spectral type, with more metal poor stars being prone to the largest differences between the PMSU and NIR spectral types. We see the same trend with metallicity as in R12: stars that are metal poor were assigned PMSU spectral types that are earlier than the NIR spectral type we assigned.
We calibrated a metallicity-sensitive function relating NIR spectral type to PMSU spectral type, to facilitate joint use of our data. We found that a linear combination of NIR spectral type and metallicity is sufficient only between NIR spectral types M1V and M3V, while a non-linear combination quali- FIG. 8.-Relation between NIR spectral type, metallicity and PMSU spectral type. The horizontal axis is the NIR spectral type determined by eye in this work. The vertical axis is the spectral type from PMSU (Reid et al. 1995;Hawley et al. 1996), determined from optical spectral features. We represent each bin as a single point, using color to indicate the mean metallicity and size to indicate the number of objects in each bin. In cases where a quarter of the stars fall into a metallicity bin different than the mean, we plot a second data point interior to the first. The area of the interior point relative to the exterior is proportional to the fraction of stars with the second metallicity. Overplotted is our best fitting relation (solid lines). We also include the best fitting linear relation (dashed lines), which extend across the region for which they were calibrated. Contours for our best fits are given by metallicities indicated in the legend and correspond to the colors used for the data points. The metallicity bins used to color data points are: tatively explains the trends seen in our data. Our best fitting non-linear relation is given by: where spectral types are expressed numerically, as described above, and [Fe/H] is given in dex. It is valid over NIR spectral types from M1V-M4V and has a scatter of half a subtype.

Spectroscopic distances
We used NIR spectral type and the H 2 O-K2 index to calibrate a relation with absolute K S magnitude, using 187 M dwarfs with parallaxes and K S magnitudes ( Figure 9). We calculated errors on absolute K S magnitude from the parallax errors, imposing a lower limit of 0.01 magnitudes (this limit was applied to three stars). We performed a linear least squares fit, using the average of the positive and negative errors on the distance to calculate the K S magnitude measurement error. The fit is valid for NIR spectral types M0V-M8V or 0.7 < H 2 O-K2 < 1.06. Expressing the M subtype numerically, our best fits are: To estimate the error on the inferred magnitudes and distances, we remove 5σ outliers and calculate the standard deviation between the measured and inferred absolute magnitudes. Outlier rejection removes four objects for the spectral type relation and three for the H 2 O-K2 relation. The standard deviation is 0.30 magnitudes for the NIR spectral type relation and 0.27 magnitudes for the H 2 O-K2 relation, indicating that most of the scatter is intrinsic, rather than due to binning by spectral type. Using standard Gaussian error propagation, we estimated that the uncertainty in the distances inferred using Equation 5 is approximately 14%. Spectroscopic distance estimates based on the H 2 O-K2 index are included for stars in our sample in Table A1. For binaries where only the total magnitude of two components is available, we assume they contribute equally to the luminosity.

METALLICITY CALIBRATION
We calibrated our metallicity relation using M dwarfs in CPM pairs with FGK stars, where the primary has a measured metallicity ( §6.1). Our method of identifying CPM pairs and additional validation using radial velocities and spectroscopic distance estimates is described in §6.2. We searched the NIR for suitable tracers of metallicity ( §6.3) and looked into potential sources of bias ( §6.4). We tested our calibration using M dwarf-M dwarf binaries and M dwarfs observed at multiple epochs ( §6.5) and compared measurements from R12 to those from this work ( §6.6).

Metallicities of the primary stars
For our potential primary stars, we used FGK stars with metallicities measured by  hereafter VF05), Santos et al. (2004Santos et al. ( , 2005Santos et al. ( , 2011, hereafter San-tos+), Sousa et al. (2006, and Bonfils et al. (2005). We use VF05 metallicities where available. We also considered those stars with metallicities measured from Sozzetti et al. (2009). VF05 andSozzetti et al. (2009) fit an observed spectrum to a grid of model spectra (Kurucz 1992 Sozzetti et al. (2009) was available: LSPM J0315+0103, LSPM J1208+2147N, LSPM J1311+0936 and PM I16277-0104. These M dwarfs are useful in extending the calibration regime to lower metallicities, but the scatter in their measured metallicities was large enough to be of concern, so we did not use these M dwarfs as part of our final calibration sample. However, we did use these four stars to validate the extrapolation of our calibration to [Fe/H] = −1.0 dex.
We used 0.03 dex divided by the square root of the number of spectra analyzed as the error for VF05 metallicity measurements, as described by the authors (typically 1-2 spectra were analyzed in VF05). Errors for metallicities from San-tos+, Sousa+, and Bonfils et al. (2005) were reported individually in the literature. Since the errors were consistent with the scatter we find between VF05 and these measurements, we did not further inflate the error bars.

Identification of calibrators
We used calibrators from previous works (Bonfils et al. 2005;Johnson & Apps 2009;Schlaufman & Laughlin 2010;Terrien et al. 2012), but also identified new calibrators. To locate new FGK-M CPM pairs, we cross-matched the LSPM-North and LSPM-South (Lépine, private communication) catalogs with themselves and with those stars with measured metallicities from VF05, Sousa+, Santos+, or Sozzetti et al. (2009). Our search was subject to the following requirements: the secondary must be within 5 ′ , have colors consistent with an M dwarf (V − J > 2, J − K S > 0.6 and H − K S > 0.1), and have proper motions within 6σ of the primary, where the uncertainties were assumed to be those stated in the LSPM catalogs.
We used a χ 2 statistic to identify CPM pairs. The statistic was constructed from the angular separation (a), the difference in proper motions (∆PM = |P M 1 − P M 2 |), and the difference in distance modulii (∆DM = DM 1 − DM 2 ). For the distance, we used parallaxes where available, and otherwise used the M J versus V − J relation (Lépine 2005) using the V − J estimates from Lépine & Shara (2005): We required χ 2 < 15 for selection of an object as a candidate binary. We note that the M J values estimated from Lépine & Shara (2005) V − J measure were often highly uncertain, because many were derived from photographic estimates of the V magnitude. Thus, the constraints from requiring a common distance modulus are weak in these cases. Additionally, the LSPM catalogs gave the same proper motion value for many very close systems where separate values could not be measured; our analysis assumed that the proper motions were independently measured.
After gathering our observations, we checked that the RV of the primary was in agreement with our measurement of the RV of the secondary and that the distance to the primary was in agreement with our spectroscopic distance estimate for the secondary. We compared the RV and distance measurements for each calibrator and three stars were immediately obvious as outliers. Two have RVs differing by more than 10σ: Gl 806.1B and CE 226. One has a distance differing by 7.5σ: HD 46375B. (This star is noted on SIMBAD as not being a CPM pair, although in MEarth imaging they do appear to move in tandem). LP 731-76, a mid M dwarf, has the same K S magnitude as its primary, an early K dwarf, clearly indicating that these are not associated. We did not include these four stars in our final sample of calibrators. While some of these systems may be physically associated, unresolved hierarchical triples, we consider the purity of our sample to be more important than its completeness.
Two of the remaining calibrators are concerning, but we do not have sufficient cause to exclude them from our sample. LSPM J0045+0015N has a distance estimate of 22pc (compared to 41pc for the primary) and an RV of 16 km s −1 (compared to 32 km s −1 ). 2MASS J03480588+4032226 has a distance estimate of 30pc (compared to 50pc for the primary) and an RV of 0 km s −1 (compared to −10 km s −1 ); the low proper motion of this object means that the evidence for the physical association of the pair from proper motion alone is weakened.
We identified 2MASS J17195815-0553043 as a visual double, and a comparison between the National Geographic Society-Palomar Observatory Sky Survey and 2MASS indicates the pair likely has a common proper motion. The distance estimates and radial velocities of the components also support the pair being physical associated. To estimate the distance to 2MASS J17195815-0553043, we assumed the two components had equal magnitudes such that the sum of their fluxes matched the published value. PM I14574-2124W (Gl 570BC) is a known spectroscopic binary, comprised of 0.6M ⊙ and 0.4M ⊙ components (Forveille et al. 1999). As we demonstrate below, the Na line we use to measure metallicity appears to be only weakly sensitive to temperature over the spectral type range of our calibration, and therefore the EWs should not be strongly influenced by the presence of a binary companion, and this object was not removed from the calibration sample. To be consistent with our treatment of known and unknown spectroscopic binaries, we use the total magnitude of PM I14574-2124W when estimating its distance.
The M dwarf calibrators and our observations are presented in Tables 2 and 3. 46 FGK-M CPM pairs appear in these tables. As previously stated, four of these objects were removed from our final calibration sample because they may not be physically associated. An additional four M dwarfs with measurements of the primary star's metallicity from Sozzetti et al. (2009) were not included in the calibration sample, although we used them to validate our calibration to lower metallicities. Two M0V dwarfs were also not included in our final metallicity calibration, as is discussed in subsequent sections. Our final calibration sample therefore consisted of 36 M dwarfs with NIR spectral types from M1V to M5V, with one M7 dwarf, and metallicities between −0.7 and +0.45 dex. The typical calibrator is an M4 or M5 dwarf and has a metallicity within 0.2 dex of solar.  Table 4. To measure the EW of a feature, we first mitigated the effect of finite pixel sizes by linearly interpolating each spectrum onto a ten-times oversampled wavelength grid with uniform spacing in wavelength. The continuum was estimated by linear interpolation between the median fluxes of the two continuum regions. We then applied the trapezoidal rule to numerically integrate the flux within the feature. We also measured ten spectral indices. We considered three indices quantifying the deformation in the continuum due to water absorption: the H 2 O-K2 index, introduced in §4 (R12), the H 2 O-H index (Terrien et al. 2012) and the H 2 O-J index (Mann et al. 2013). We also measured the flux ratios defined by McLean et al. (2003) and used by Cushing et al. (2005). These ratios quantify absorption in several water, FeH and CO bands. The indices we measured are summarized in Table 5. Finally, we considered three non-linear combinations of parameters. The non-linear combinations we considered were motivated by previous work: Luhman & Rieke (1999) suggested that EW Na /EW CO is temperature-sensitive and R12 used the ratios EW Na / (H 2 O-K2) and EW Ca / (H 2 O-K2) to fit their metallicity relation.
We searched for the combination of three parameters that provide the best fit to metallicity, using the forms: where F n is the EW of one of the 27 spectral features in Table 4, one of the ten indices in Table 5, or one of the three non-linear combinations of parameters described above. We  b Astrometry references. If one reference is provided, it applies to both position and proper motion; if two are provided, the first is for position and the second for proper motion. c Apparent KS magnitudes are from S06. d Errors on the distance estimates are 14%. e Proper motions and distances for primary stars are from Hipparcos (van Leeuwen 2007) except when otherwise noted. f For CE 226 and Gl 376B, the Hipparcos proper motion for the primary was found to be a better match to the observed motion of the secondary from 2MASS to recent epoch MEarth imaging than the proper motion given in Ruiz et al. (2001, for CE 226) or in LSPM-North (for Gl 376B). In these cases, the Hipparcos value has been adopted in the table. g Lépine, private communication. We resolved this object as a binary. An appended "J" indicates a measurement that was derived for the components jointly. We assume the two components contribute equally to the luminosity in order to estimate their spectroscopic distances. h No parallax was available for the primary. Its distance was estimated assuming an absolute KS magnitude of 6, typical for an early K dwarf. i No KS magnitude could be found for Gl 806.1B. We estimated a rough magnitude from 2MASS Atlas images using a 4 pixel aperture radius (this value was chosen to reduce contamination from nearby stars), and applied an aperture correction of 0.04 magnitudes, derived from stars of similar K magnitude elsewhere in the field.    Cushing et al. (2005), but feature and continuum windows were defined based on our observations. b Feature and continuum windows were modified from those defined in Covey et al. (2010). c Atomic feature not identified. used the RMSE as a diagnostic to identify the best potential metallicity relations. There were a multitude of relations with RMSE < 0.14 dex, of which the majority included the EW of the Na line at 2.2µm as the primary indicator of metallicity (sometimes appearing as EW Na /EW CO or EW Na / (H 2 O-K2)) and included a quadratic term. However, we preferred the two-parameter fit [Fe/H] = A (EW Na ) + B (EW Na ) 2 + C because it uses one fewer parameter. Motivated by the clear trend with metallicity seen in EW Na , we also considered functional forms other than a quadratic, including a spline. No other forms tested resulted in a statistically superior fit. We show our result in Figure 11.
In performing our final fit, we did not include the two K7/M0V stars. In Figure 11 these are evident as having an EW Na lower than other calibrators of similar metallicity. As discussed in §6.4, we attempted to find a metallicity relation that was valid through these early spectral types, but did not converge on a suitable result. Our final calibration sample therefore includes 36 M dwarfs with spectral types M1V and later. We address this choice in detail in the following section.
Our best fit is given by: We estimated the error introduced by our limited number of calibrators by bootstrapping. We randomly selected 36 of our calibrators, allowing repeats, and re-fit our metallicity relation. The standard deviation of the difference between the best fitting metallicities of the M dwarf secondaries and the metallicities of the primaries, averaged over 100 bootstrap samples, was 0.12 ± 0.01 dex. The correlation coefficient, R 2 ap is often used to evaluate the goodness of fit. The correlation coefficient indicates how well the fit explains the variance present in the data and is given by: where n is the number of data points and p is the number of parameters. The R 2 ap value for our fit is 0.78±0.07. The bestfitting metallicities for our calibrators are included in Table 3. The errors on metallicity include the errors on EW Na , bootstrap errors and the scatter in our best fit, added in quadrature. We took the bootstrap errors to be the 1σ confidence interval on the resulting metallicities when considering the best fits from 100 bootstrap samples. The intrinsic scatter in the relation (0.12 dex) dominates for all but the lowest metallicity stars.
The scatter in our metallicity relation is similar to those reported by R10, R12, Terrien et al. (2012) and Mann et al. (2013) despite differences in sample size, lending support to the idea that the scatter is astrophysical in origin. We consider potential temperature and surface gravity effects in §6.4. One possibility is variations between the Na abundance and [Fe/H] of the primary solar-type star. We considered whether an M dwarf's EW Na is a better tracer of its primary star's Na abundance than its Fe abundance. 32 of our calibrators have measured abundances for Na from VF05. We related the spectral features and indices in Tables 4 and 5 to the Na abundance of the primary star. We found several suitable tracers; however, none reduced the scatter.
In Table A1, we include the EWs of the Na line at 2.20µm and the Ca line at 2.26µm, the H 2 O-K2 index, our inferred [Fe/H], and their associated errors for each of our targets. The corresponding values for the FGK-M CPM pairs can be found in Table 3. 6.4. Influence of effective temperature and surface gravity on the metallicity calibration We examine the potential influence of differences in the effective temperature and surface gravity on the metallicity calibration presented in §6.3 by computing EW Na for a grid of BT-Settl theoretical spectra for spectral types K5V-M7V, shown in Figure 12 (Allard et al. 2011, the behavior of NIR lines in theoretical spectra are discussed in some detail in R12). The spectral type range corresponds to approximately K5V-M6.5V, depending on the adopted temperature scale (we quote the temperature scale from E. Mamajek, which is available online. 6 ). The BT-Settl theoretical spectra show EW Na varying by 1Å between M0V and M8V stars (Figure 12). We also note that in our K-band SpeX spectral sequence ( Figure  4) the Na line at 2.2µm is broader for the latest spectral types.
We plot in Figure 13 the median EW Na for each NIR spectral type as a function of H 2 O-K2, for two subsamples. Our "nearby sample" ( §2.2) formed the first, and kinematically young stars (V tot < 50km/s) formed the second. We selected the nearby sample to approximate a volume limited sample, which is unlikely to be influenced by selection effects that may exist in the rotation sample. We selected the kinematically young sample in order to isolate stars that are expected to be of similar age and metallicity. We found a similar increase in the median EW Na of mid to late M dwarfs as we noted in the theoretical spectra. This could introduce a systematic error of 0.1 dex in the metallicities of early M dwarfs relative to mid M dwarfs. However, we are uncertain of the origin of this effect, given the differing behavior of our two subsamples and the relative differences in the number of early and late type stars (there are 23 stars with NIR spectral types M0V-M2V and 231 with spectral types M4V-M5V across the two subsamples).
We considered whether an alternative parameterization could account for this potential bias. We show the residuals for our chosen parameterization and three alternatives, including the parameterization from R12, in Figure 14. In Figure  15, we show the effect that the alternative calibrations have on the metallicities of the sample as a whole. With the R12 parameterization, the inferred metallicities of the latest stars decreased by 0.1 dex and metallicities were consistent across The behavior of EW Na in our observed spectra. We plot the median EW Na against the median H 2 O-K2 for each NIR spectral type as shown in Figure 11. The medians for two subsamples are shown. Filled squares include only those stars which are in our nearby sample and open squares include only kinematically young stars. Points are colored by their NIR spectral type, from purple for M0V stars to red for M8V stars, as shown in Figure 11. spectral types. However, the metallicities of M5 were lowered relative to those of M4 dwarfs, the spectral range across which our relation is best calibrated. Furthermore, the fit is unconstrained at the latest spectral types where the choice of the R12 parameterization makes the most difference. Including the EW of magnesium or the H 2 O-K2 as a third parameter in the metallicity calibration improves the fit for the two K7/M0V calibrators and has only a marginal effect at other spectral types. However, only scatter above the best fit plotted in Figure 11 was reduced in this case, while the scatter below our best fit remained.
When the M0V calibrators were not included in the fit, the addition of these extra parameters makes little difference. Therefore, rather than including an additional parameter to fit these two points at the far end of our spectral type range, we simply limit our calibration to a range of spectral types which appear to be well-fit by a relation depending solely on Na.
The insensitivity of NIR spectral types to late K dwarfs may  15.-The median metallicity for two subsamples of stars as a function of NIR spectral type. Filled squares indicate median metallicities for stars without measured rotation periods and open squares indicate the median metallicities for kinematically young stars. Points are colored by their NIR spectral type, from purple for M0V stars to red for M8V stars, as shown in Figure 11. be partially responsible for the behavior seen in our two M0V calibrators. The optical spectral type of PM I07400−0336 places it as K6.5V dwarf (Poveda et al. 2009) and LSPM J1030+5559 has been identified previously as a K7V dwarf (Garcia 1989). However, theoretical models indicate that the EW Na should remain constant between late M and mid K dwarfs (with slight dependence on surface gravity), and -We compare EW Na against EW CO for all stars in our sample. According to Luhman et al. (1998), very young stars would reveal themselves through low EW Na but high EW CO . We have no data in the upper left corner of this plot, indicating it is likely that no very young stars are present in our data. Mann et al. (2013) reported a metallicity calibration that is valid from K5V-M5V.
Surface gravity remains one possible explanation for the K7/M0V discrepancy and has yet to be explored in the context of empirical calibrations. Luhman et al. (1998) demonstrated that in the low surface gravity environments of very young stars, Na may appear abnormally weak. It is therefore possible that an M dwarf with an age of several Myr could be masquerading as a metal-poor object. The CO band head is sensitive to gravity in the opposite manner and is therefore a useful indicator of youth (Luhman et al. 1998). In Figure  16, we plot EW Na against EW CO for all stars in our sample. We found a general positive correlation and spectral dependence, but no object stood apart has having low EW Na but high EW CO . This is not surprising as it is unlikely that we would find a new, bright young star within 25pc.
We considered the potential for other systematics by comparing the difference between our best fitting metallicities and the metallicities of the primaries to the EWs of all other indices. In all cases, we found no significant systematic effects. 6.5. Tests of our metallicity relation As a test of our metallicity calibration, we compared the metallicities we estimated for the components of M dwarf-M dwarf CPM pairs. We have observed 22 such pairs. 11 were placed on the slit together and so share observing conditions, while 11 were observed separately but close in time. In both cases, the two stars were reduced with the same telluric standard. In Figure 17 we show the results of this comparison. The mean metallicity difference between the primary and secondary components is −0.01 dex with a standard deviation of 0.05 dex. This is less than the uncertainty of our metallicity measurement by a significant amount, lending support to the idea that most of the scatter in the metallicity relation is astrophysical in origin, as mentioned in §6.3.
We also compared EW Na measurements for stars that were observed on more than one occasion in Figure 1 (see  §4). We found that our EW Na measurements were consistent even for observations taken in very different conditions and separated in time by months or more. The mean EW Na difference between the observation we elected to keep and the observation we discarded was −0.01 dex with a standard deviation  (Herter et al. 2008). To facilitate joint use of our observations and those from R12, we determined the relationship between TripleSpec and SpeX EWs. We compare our EW Na measurements directly in Figure 18. We used the following relation to convert from TripleSpec to IRTF EW Na : EW Na,N13 = 0.036 + 0.90 (EW Na,R12 ) Similarly for the Ca line at 2.26µm: We also directly compared our metallicity estimates for the 28 stars in common (excluding metallicity calibrators). As seen in Figure 18, the metallicity measurements agreed well for sub-solar metallicities, but for metallicities above solar, the relation in this work gives higher metallicities for late M dwarfs (M5V-M7V). The difference between our inferred metallicity and that from R12 is 0.0±0.07 dex for M1V-M4V stars and 0.08 ± 0.05 for M5V-M7V stars. This difference is consistent with the effects discussed in §6.4, but we note that our relation is most strongly constrained for M4 and M5 dwarfs.
The objects observed by R12 are listed in Table A2. We have included EWs updated using Equations 14 and 15. After applying our EW Na relationship, we can directly compute the metallicities for stars published in R12 using our metallicity calibration. We also present these updated metallicities in Table A2 FIG. 18.-Comparison between our measurements and those from R12. In the top panel, we compare EW Na measured in this work using the SpeX instrument on IRTF to those presented in R12, who used the TripleSpec instrument on Palomar. We show the one-to-one line (dashed line) and our best fit (solid line). In the bottom panel, we compare [Fe/H] estimated in this work directly to that estimated by R12. We over plot the mean metallicity difference for an early subsample (NIR spectral types M1V-M4V) and a late subsample (M5V-M7V). Data are plotted as filled squares if our EW Na measurements agree within the errors and as open squares if the discrepancy is larger than the associated error. In both panels, data are colored by their NIR spectral type, from purple for M0V stars to red for M8V stars, as shown in Figure 11.
We exploited our sample of M dwarfs with spectroscopically determined NIR metallicities to identify which color-color diagrams are metallicity sensitive and to derive an empirical relationship between an M dwarf's NIR color and its metallicity. In Figure 19, we plot JHK S color-color diagrams for the 444 of our targets with the highest quality 2MASS magnitudes (Skrutskie et al. 2006, qual flag=AAA). We also plot the Bessell & Brett (1988) M dwarf main sequence, which coincides with our solar metallicity stars. These diagrams are plotted in the 2MASS photometric system; we used the color transformations updated 7 from Carpenter (2001) to transform the colors from Bessell & Brett (1988) to the 2MASS system.
All color combinations discriminated effectively between low and high metallicity stars. Consistent with Johnson et al. (2012), we found that the J − K S color of an M dwarf is the best single-color diagnostic of its metallicity. We used the vertical (J − K S ) distance from the J − K S , H − K S Bessell & Brett dwarf main sequence (D MS ) as the diagnostic for the metallicity of an M dwarf. We considered using D MS to determine both EW Na and [Fe/H] directly ( Figure 20). We chose to relate D MS to EW Na because the correspondence is linear and because it relates two directly measured quantities.
We determined the relation between EW Na and D MS using those stars with 2.5 < EW Na (Å) < 7.5 and |D MS | < 0.1. We binned the data into 0.5Å-wide bins and computed the median D MS in each. We then fit a straight line through these points, using the reciprocal square root of the number of data points in each bin as the weights. The best-fitting relation between EW Na and D MS , shown in Figure 20 is: We show the resulting photometric metallicity calibration in Figure 21. Our calibration is valid from 2.5 < EW Na (Å) < 7.5, corresponding to −0.7 < [Fe/H] < 0.3 and for 0.2 < H −K S < 0.35. The 1σ uncertainty in EW Na translates to 0.1 dex for EW Na = 7Å and 0.5 dex for EW Na = 3Å. This calibration is particularly useful because it does not require V magnitudes, which are often unreliable, or parallaxes, which are often unavailable. In contrast, accurate JHK S magnitudes are available for the majority of nearby M dwarfs from 2MASS.

RADIAL VELOCITIES FROM NIR SPECTROSCOPY
Absolute wavelength calibration for moderate resolution NIR spectra are typically determined using a lamp spectrum taken at the same pointing as the science spectrum, as done by Burgasser et al. (2007), who measured the radial velocity of an L dwarf binary to 18 km s −1 accuracy using SpeX (R ≈ 2000). An alternative is to take deep sky exposures and use OH emission lines to perform wavelength calibration. This approach was used, for example, by Muirhead et al. (2013), who use the TripleSpec instrument on Palomar (R ≈ 2700) to measure absolute radial velocities for the eclipsing post common envelope binary KOI-256 with typical errors of 4 km s −1 .
We acquired Thorium-Argon spectra regularly throughout the night to track instrumental variations, but it was not possible to obtain them at every telescope position due to the exposure times required. We found that this procedure was not adequate for accurate radial velocity work. We therefore used telluric absorption features to supplement the wavelength calibration by adjusting the velocity zero-points for individual observations, then cross-correlated each spectrum with a standard spectrum to measure its absolute RV ( §8.1). In §8.2, we discuss using precisely measured RVs from Chubak et al. (2012) to investigate random and systematic error. We describe further tests of our method in §8.3.

Radial velocity method
Atmospheric absorption features present in our data provided a natural replacement to arc spectra. By correlating the telluric lines in our spectra with a theoretical atmospheric transmission spectrum (hereafter called simply the "transmission spectrum"), we determined the absolute wavelength calibration. The SpeXtool package includes a transmission spectrum created using ATRANS (Lord 1992). This spectrum was created using environmental parameters typical of Mauna Kea and an airmass of 1.2 and has a resolution five times that of SpeX. We used the wavelength calibration determined by SpeX using ThAr arc spectra as our initial wavelength guess for the nontelluric corrected science spectrum. From this wavelength solution, we created a wavelength vector that was oversampled by a factor of six and linearly spaced in wavelength.
We found that excellent continuum removal was required for the wavelength calibration to be determined through direct cross correlation of the science spectrum and the transmission spectrum. However, the large atmospheric features Stars with EW Na > 7.5Å are plotted in magenta. Grey points are stars of other spectral types other than the range indicated in the top panels. Overplotted are the dwarf (blue) and giant (red) tracks from Bessell & Brett (1988), converted to the 2MASS system using the updated color transformations of Carpenter (2001), which are available online. made this difficult. Instead of attempting to remove the continuum from the M dwarf and subsequently finding the offset between the stellar spectrum and the atmospheric spectrum, we tackled these problems simultaneously. We did this by finding the modifications to the transmission spectra that provided the best match the telluric features observed in the science spectrum. There were three differences between the theoretical transmission spectrum and the telluric features as observed in the science spectrum: the continuum, the strength of the telluric features and the pixel offset between the spec-tra.
The first parameter of our model was a Legendre polynomial as a function of pixel by which the transmission spectrum was multiplied in order to replicate the shape of the spectrum. The curvature of the spectrum was affected by both instrumental effects and the M dwarf spectral energy distribution. We used a 3rd or 4th degree Legendre polynomial and fit for the coefficients. We selected the degree of the polynomial by eye for each order, using the lowest degree polynomial required to model several representative M dwarf spectra.  Figure 19, while symbols indicate NIR spectral type (K7V-M3V as triangles, M4V-M5V as circles, and M6V-M9V as stars). Large filled symbols are our metallicity calibrators. Overplotted are isometallicity contours for our best fit, which relate distance from the main sequence to metallicity via EW Na .
The second parameter was an exponential scaling of the flux, to account for the effects of airmass and atmospheric water vapor on the depths of telluric features. The transmission spectrum represents typical conditions on Mauna Kea, while we observed at air masses from 1.0 to 1.7 with humidity from 85% to less than 15%.
As discussed in Blake et al. (2010), differences in airmass scale the depths of the telluric features (T ) as T = T τ 0 where the optical depth τ scales linearly with airmass. Blake et al. (2010) were able to find a single linear scaling between airmass and τ using a large sample of A0V stars. We attempted to use the same approach, but found substantial scatter and systematic differences in the scaling of different telluric features with airmass. This is likely due to the water absorption features in our spectral region, which are time-variable, and cannot be modeled by a simple function of airmass alone. We therefore chose to take an empirical approach and included the exponential scaling τ as a model parameter.
The third and final parameter was the offset in pixels between the transmission and science spectrum. We modeled the offset as linear in wavelength. To apply the shift, we cre-ated a new wavelength vector that was linearly shifted from the original and interpolated the transmission spectrum onto the new wavelength vector. We constrained the allowable range for the offset because atmospheric features appear at regular spacing and we found that if unconstrained, our fitting program can too often land in a local minimum. We used 0.0008µm as the limit, which is larger than any offset we expected. In our full sample, no shifts beyond 0.0006µm were found, and very few beyond 0.0004µm.
We modeled each order of the non-telluric corrected science spectrum independently, minimizing the difference between our model and the science spectrum using a nonlinear least squares approach, implemented through mpfit (Markwardt 2009). We determined by trial and error the region of each order to use. Regions with high signal to noise and strong telluric features but uncontaminated by strong stellar features were required for optimal performance. Because of these constraints, this method worked better in the J, H and K-bands than it did in Y or Z.
Once we determined the absolute wavelength solutions of science target and an RV standard, we interpolated the telluric-corrected spectra onto a common wavelength vector that was oversampled and uniform in the log of the wavelength (such that a radial velocity introduces a constant offset in pixels). The continuum is different in the telluric-corrected spectrum because telluric correction removed instrumental effects, so we used a spline to remove the continuum. We used xcorl to cross-correlate the two spectra and determine the offset. We used the same standard star (Luyten's star, also known as Gl 273 or LSPM J0727+0513) throughout because it had an accurately measured absolute radial velocity from Chubak et al. (2012) and a NIR spectral type in the middle of our range (M4V).
We took the final RV for each target to be the median of the RVs measured in the J, H and K-bands and applied the heliocentric correction, implemented through the IDL routine baryvel (Stumpff 1980). Our final estimate of the error is the 1σ confidence limit on the RV after 50 trials added in quadrature to 4.4 km s −1 (our internal measurement error, see §8.3). These values are reported in Table A1. This method of measuring radial velocities is applicable to other moderate resolution NIR spectrographs, including TripleSpec, and uses observations of the target star to refine the wavelength calibration. Our method is therefore likely to be useful for instruments where obtaining lamp spectra is expensive. Chubak et al. (2012) presented absolute, barycentriccorrected RVs for 2046 dwarf stars with spectral types from F to M. M dwarf RVs were measured by comparison to an M3.5V RV standard, offset to agree with the measurements from Nidever et al. (2002). No corrections were made for convective or gravitational effects for M dwarfs, and Chubak et al. (2012) report a systematic error of 0.3 km s −1 (random errors are at this level or lower in nearly all cases). Ten of their M dwarfs are in our sample. We chose one of these, LSPM J0727+0513, as our standard star. For the other nine stars, we compare our measurements to those from Chubak et al. (2012) in Figure 22. Considering the RV measured in each order separately, we found that the bluest two bands (Z and Y ) systematically underestimate (Z-band) or overestimate (Y -band) the RV. The wavelength calibration is also subject to failure in those bands. We suggest that this  Chubak et al. (2012), with NIR spectral type on the horizontal axis. Data points show the difference between our adopted RV for each star, which is the median of the RV measured in each of the J, H, and K-bands, and that reported in Chubak et al. (2012). The dashed line shows the mean difference between our measurements and those from Chubak et al. (2012). We also look at how well the RV measured from a single band compares to the values from Chubak et al. (2012); the mean difference for each band is plotted as a colored line. The Y and Z-bands tend to over-and underestimate the RV. A −2.6 km s −1 offset has been applied.

Using precise RVs to investigate errors and systematics
is because in these two orders, the strongest stellar features overlap with the strongest telluric features, compromising the wavelength calibration and therefore the velocity measurement. They were also the orders with the lowest S/N. The RVs reported in this paper are the median of the J, H, and K-band measurements.
We measured RVs for all our targets using each of the ten RV standards from Chubak et al. (2012) in order to determine our internal error and systematic RV offset. The typical standard deviation of RVs measured against an alternative standard relative to that measured against LSPM J0727+0513 was 4.2 km s −1 . We used this value as our internal random error. RVs measured using LSPM J0727+0513 were systematically higher than those measured using other RV standards. Considering M3V-M5V standards, the median offset was 2.6 km s −1 with a standard deviation of 1.5 km s −1 . The values reported in this paper include a −2.6 km s −1 systematic RV correction. Our total internal measurement error is 4.4 km s −1 , which is our internal random error (4.2 km s −1 ) added in quadrature to our internal systematic error (1.5 km s −1 ).
Our choice of a single, mid-M RV standard does not appear to systematically affect the RV measurements or errors of early and late M dwarfs at this level of precision. We investigated the effect of the standard spectral type by comparing the results using LSPM J0727+0513 with using an M2V star, PM I06523-0511 (Gl 250), to measure the RVs of early Mdwarfs, and an M7V star, J1056+0700 (Gl 406), to measure the RVs of late M-dwarfs, finding that these choices did not appear to systematically affect the measured RVs, and that the scatter remained consistent with our estimated uncertainties.

Validating the use of SpeX for radial velocities
To determine the precision of our wavelength calibration method, we used the transmission spectrum to create simulated data in each order, which we then calibrated. We simulated stellar absorption lines of random widths, depths and locations on top of the transmission spectrum and multiplied by a polynomial (drawn from a random distribution) to curve the data. We then offset the spectrum and monitored how well we could recover that offset. The accuracy declined as more stellar absorption lines were added to the spectrum. With 50 added lines, accuracy was better than 5 km s −1 in all orders and better than 1 km s −1 in H-band.
We have multiple observations for 26 stars at different epochs. The time between observations ranges from days to months to years. We compared our RVs for these stars (Figure 23). The mean difference between the observation we elected to keep and the observation we chose to discard with 0.08 km s −1 with a standard deviation of 4 km s −1 , consistent with our calculation of the error.
Finally, we compared the RVs of CPM pairs ( Figure 24). 11 of these stars are separated and were observed independently and 11 were observed together on the slit. These observations were taken close in time, at near-identical conditions and were reduced using the same wavelength calibration and telluric standard. The mean RV difference between the primary and secondary components is 0.2 km s −1 with a standard deviation of 2 km s −1 .

CONCLUSIONS
The MEarth team and collaborators are creating a wellstudied sample of nearby M dwarfs which will be the basis for future studies investigating their fundamental properties, their evolution, and the exoplanets orbiting them. The data set being assembled is diverse, with photometric rotation periods, parallaxes, and optical spectra. In this work, we presented metallicities, NIR spectral types and radial velocities for a fifth of MEarth M dwarfs.
We created a NIR spectral typing routine, determined byeye spectral types and presented spectral standards for M1V-M8/9V dwarfs. We related NIR spectral type to PMSU spectral type, finding the conversion to be metallicity-sensitive. We calibrated a new spectroscopic distance relation using NIR spectral type or H 2 O-K2, which can be used to estimate distances to 14%.
We used M dwarfs in CPM pairs with an F, G or K star of known metallicity to calibrate an empirical metallicity relation. We validated the physical association of these pairs using proper motions, radial velocities and distances (making use of our RV measurements and spectroscopic distance estimates for the secondaries). We explored the NIR for combinations of EWs that effectively trace stellar metallicity, and found that the EW of the Na line at 2.2µm is sufficient. Our metallicity calibration has a standard deviation of 0.12 dex and R ap = 0.78. It is calibrated using 36 M dwarfs with NIR spectral types from M1V to M5V and −0.6 < [Fe/H] < 0.3, and can be extrapolated to [Fe/H] = −1.0 dex. We found no evidence that the calibration breaks down for M dwarfs as late as M7V.
Using our EW Na measurements of 447 M dwarfs and the J − H, H − K S color-color diagram, we calibrated a relationship between an M dwarf's distance from the Bessell & Brett main sequence and its sodium equivalent width. It is valid from 2.5 < EW Na (Å) < 7.5. The standard deviation of our fit is 2Å and has an R 2 ap value of 0.92. Metal-rich M dwarfs can be selected by taking those M dwarfs whose J − K S colors are redder than the Bessell & Brett (1988) M dwarf track in the J − H, H − K S color-color diagram.
We developed a method to wavelength calibrate SpeX M dwarf spectra using telluric features present in the data, and we measured absolute radial velocities for the stars in our sample at a precision of 4.4 km s −1 . We used synthetic spectra, M dwarfs with precise radial velocities from Chubak et al. (2012) and M dwarf-M dwarf binaries to validate our method. Because telluric absorption features are strong in even short exposure data, our method for determining the absolute wavelength calibration requires no information beyond the science spectrum itself. This opens up the possibility of measuring radial velocities for stars with an extant moderate resolution NIR spectrum.
Our measurements, including NIR spectral types, EWs, radial velocities, and spectroscopic distance estimates are presented in Table A1. We also include distances estimated from parallaxes, and radial velocities from PMSU. To facilitate joint use of our datasets, we reproduce spectral measurements for M dwarfs observed by R12 in Table A2, with EWs modified to account for differences between their TripleSpec and our IRTF measurements and [Fe/H] inferred using our calibration; we also include PMSU spectral types and RVs, and the parallaxes reported in R12.
In future work, will continue to explore the use of the NIR as a diagnostic of intrinsic stellar properties, investigating how metallicity relates to rotation period, tracers of magnetic activity, and galactic kinematics.
ERN is supported a National Science Foundation Graduate Research Fellowship. The MEarth team gratefully acknowledges funding from the David and Lucile Packard Fellowship for Science and Engineering (awarded to DC). We thank S. Dhital, A. Dupree, M. Holman, and A. West for helpful conversations. This material is based upon work supported by the National Science Foundation under grant number AST-0807690 and AST-1109468. Based on observations at the Infrared Telescope Facility, which is operated by the University of Hawaii under Cooperative Agreement no. NNX-08AE38A with the National Aeronautics and Space Administration, Science Mission Directorate, Planetary Astronomy Program. This research has made extensive use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center / California Institute of Technology, funded by NASA and the NSF, NASAs Astrophysics Data System (ADS), and the SIMBAD database, operated at CDS, Strasbourg, France. APPENDIX Tables A1 (All M dwarfs from our rotation and nearby samples and the potential calibrators) and A2 (M dwarfs observed by R12) are available online and in the refereed version of this article. These tables contain positions, proper motions, spectral measurements, measured radial velocities and those from the literature, estimated distances and those from the literature, and inferred [Fe/H].