{'fraction_non_alphanumeric': 0.028045863235172813, 'fraction_numerical': 0.019178421182875525, 'mean_word_length': 4.440056839428614, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Retrieval of Precise Radial Velocities from Near-Infrared High Resolution Spectra of Low Mass Stars

Given that low-mass stars have intrinsically low luminosities at optical wavelengths and a propensity for stellar activity, it is advantageous for radial velocity (RV) surveys of these objects to use near-infrared (NIR) wavelengths. In this work we describe and test a novel RV extraction pipeline dedicated to retrieving RVs from low mass stars using NIR spectra taken by the CSHELL spectrograph at the NASA Infrared Telescope Facility, where a methane isotopologue gas cell is used for wavelength calibration. The pipeline minimizes the residuals between the observations and a spectral model composed of templates for the target star, the gas cell, and atmospheric telluric absorption; models of the line spread function, continuum curvature, and sinusoidal fringing; and a parameterization of the wavelength solution. The stellar template is derived iteratively from the science observations themselves without a need for separate observations dedicated to retrieving it. Despite limitations from CSHELL's narrow wavelength range and instrumental systematics, we are able to (1) obtain an RV precision of 35 m/s for the RV standard star GJ 15 A over a time baseline of 817 days, reaching the photon noise limit for our attained SNR, (2) achieve ~3 m/s RV precision for the M giant SV Peg over a baseline of several days and confirm its long-term RV trend due to stellar pulsations, as well as obtain nightly noise floors of ~2 - 6 m/s, and (3) show that our data are consistent with the known masses, periods, and orbital eccentricities of the two most massive planets orbiting GJ 876. Future applications of our pipeline to RV surveys using the next generation of NIR spectrographs, such as iSHELL, will enable the potential detection of Super-Earths and Mini-Neptunes in the habitable zones of M dwarfs.


Introduction
The radial velocity (RV) technique has been extremely successful in the detection and confirmation of exoplanets (Latham et al. 1989;Mayor & Queloz 1995). In keeping with our pursuit of Earth-like worlds, the majority of RV targets have been FGK stars similar to our Sun (Akeson et al. 2013). However, this neglects M dwarfs, which make up 75% of the stars in the Solar neighborhood (Henry et al. 2006). M dwarfs' habitable zones (HZs) are also at smaller semi-major axes due to their lower luminosities; coupled with their lower masses, this results in larger radial velocity perturbations from a companion planet in the HZ than for FGK stars. Surveys of close-by early-and mid-M dwarfs have reached RV precisions down to ∼ 2 m s −1 , sufficient to detect Super-Earths in these stars' HZs (Zechmeister et al. 2009;Bonfils et al. 2013); by comparison, a twenty-fold increase in RV precision is needed to find Earth-mass planets in the HZs of Sun-like stars.
The detection of planets around M dwarfs is hampered by two important factors, however. First, K and M stars are generally more active than G stars, with late-M dwarfs being the most active (West et al. 2004;Basri et al. 2010). Stellar activity, such as starspots, can produce false positive planet signatures by introducing additional power at the frequency of the rotational period of the star and aliases thereof in the Doppler measurements (Robertson et al. 2015;Vanderburg et al. 2015). Second, the low luminosity of M dwarfs at optical wavelengths makes it difficult to obtain high signal to noise (SNR) measurements. A logical solution to these problems is to observe at near infrared (NIR) wavelengths. Not only do M dwarfs emit most of their light in the NIR, the The lower RV precision of NIR RV surveys is partly attributable to the lack of a wavelength calibration method on par with those of optical surveys. Precise wavelength calibration is currently obtained either through extreme environmental stabilization (e.g. HARPS), or the use of a simultaneous common optical path wavelength reference. The latter process involves the comparison of the observations with a set of "standard" well-calibrated spectra. In the previously mentioned NIR RV studies, telluric lines were used for calibration, which were imprinted on the incoming star light from its passage through the Earth's atmosphere; comparisons between the observed telluric lines with models or observations of telluric spectra then provide a wavelength solution. In optical surveys, precise wavelength calibration is realized by placing an Iodine gas cell in the telescope optical path (Marcy & Butler 1992;Butler et al. 1996). As the incoming star light passes through the gas cell, it picks up the spectral signatures of its gas. Comparing these observed spectra with a precisely measured absorption spectrum of the gas cell at higher spectral resolution than the astronomical spectrograph (R ∼ 500,000 obtained with a Fourier Transform Spectrometer) then yields a wavelength solution. Placing the gas cell in the optical path of the telescope also ensures that the spectral lines are obtained simultaneously and under the same physical conditions (i.e. temperature, pressure, humidity) as the instrument and the telescope, which is not the case for telluric lines. Similar gas cells have recently been developed for the NIR, such as the ammonia gas cell used by Bean et al. (2010) on the CRIRES spectrograph at the VLT to obtain an RV precision of ∼ 5 m s −1 for late-M dwarfs. Gas cells filled with hydrocarbons, chlorocarbons, HCl, and isotopologues of HCN, C 2 H 2 , CO, and CH 4 have also been considered (Mahadevan & Ge 2009;Valdivielso et al. 2010;Anglada-Escudé et al. 2012;Plavchan et al. 2013a).
An additional issue with current NIR RV surveys is that they rely on large (8-10 m aperture) telescopes, which generally have fewer nights available for RV monitoring. In contrast, smaller telescopes can provide the high cadence necessary for precise RV observations. For example, Prato et al. (2008); Crockett et al. (2011Crockett et al. ( , 2012; Davison et al. (2015) were able to use the smaller (3 m aperture) NASA Infrared Telescope Facility (IRTF), in conjunction with telluric line and emission lamp wavelength calibration, to obtain RV precision of 60-130 m s −1 for M dwarfs in the K band using the CSHELL spectrograph, which has a spectral resolution R ∼ 46000 (Greene et al. 1993). In addition, many of the instruments currently in development for dedicated NIR RV surveys will be installed on smaller aperture telescopes, e.g. CARMENES/Calar Alto Astronomical Observatory (Quirrenbach et al. 2010), SPIRou/Canada France Hawaii Telescope (Santerne et al. 2013), and iSHELL/NASA Infrared Telescope Facility (Rayner et al. 2012). The two former instruments will rely on extreme stabilization, whereas the latter will use gas cells for stable wavelength calibration. It is essential that we prepare for these next generation NIR RV surveys by focusing on data reduction and analysis methods that can process their results.
In this paper, we describe in detail the NIR RV extraction pipeline used in our RV survey of M dwarfs conducted using the NASA IRTF CSHELL spectrograph (Greene et al. 1993), and presented in Anglada-Escudé et al. (2012); Plavchan et al. (2013a); Gagné et al. (submitted to ApJ). This pipeline uses the "Grand Solution" approach, which derives the stellar template, wavelength solution, and RVs simultaneously (Valenti & Anderson 2010) in an iterative fashion. This differs from other methods used in some previous RV surveys, such as substituting imperfect stellar spectral models as stellar templates. Thus, our pipeline offers a novel way to extract RVs from raw spectra from past, present, and future observation campaigns in the NIR. We test our pipeline on the RV standard GJ 15 A to assess the RV precision obtainable for early-M dwarfs, as well as the M Giant SV Peg to evaluate the pipeline's precision at the high SNR limit. We also test it on the planet host GJ 876 to validate its planet detection capabilities.
We outline our observations of GJ 15 A, SV Peg, and GJ 876 in §2 and describe our data reduction process in §3. In §4 we detail our RV extraction pipeline and in §5 we present our results. We summarize our work and discuss future prospects in §6.

Observations
Three targets were picked from our NIR RV survey to test the RV stability, high SNR behavior, and planet detection capabilities, of the RV pipeline presented here. These targets are summarized in Table 1. GJ 15 A was chosen as a RV standard star due to its brightness, as well as its established RV scatter of a few m s −1 due mostly to an orbiting exoplanet (Endl et al. 2006;Howard et al. 2014), which is below our estimated level of precision (Anglada-Escudé et al. 2012). SV Peg is an M giant star that has already been used as a reliable high SNR target in Anglada-Escudé et al. (2012), and therefore will be used for our high SNR tests as well. GJ 876 has four confirmed planets, with the most massive two causing RV amplitudes ∼ 200 m s −1 , which should be easily detectable with our pipeline (Delfosse et al. 1998;Marcy et al. 1998Marcy et al. , 2001Rivera et al. 2005Rivera et al. , 2010. Data were taken on the 3 m telescope at NASA IRTF from September 2010 to August 2012 for GJ 15 A and GJ 876, and to August 2011 for SV Peg 1 . Table 2 gives the minimum, maximum, median, and total SNR per pixel of the observations of each night, as well as the number of observations (N obs ) obtained that night for the three targets. Single order spectra in the K band (2.309-2.315 µm) were obtained using the CSHELL spectrograph with a 0. 5 slit and a spectral resolution R ∼ 46000. A circular variable filter (CVF) was used for order selection. A methane isotopologue ( 13 CH 4 ) gas cell with 90% continuum throughput was placed into the beam to achieve wavelength calibration (Anglada-Escudé et al. 2012;Plavchan et al. 2013a). To set the wavelength scale at the start of observations, we observed A-type stars with the gas cell in the optical path of the telescope and matched the deepest gas cell line to pixel column 179 on the CSHELL detector. For the brightest A stars a latent image of the target was sometimes left on the detector after switching from imaging mode (used to position the target in the slit) to spectroscopic mode; these spectra were flagged during the data reduction process and discarded. Around 15 flat fields and darks were taken in spectroscopic mode at the end of each night of observations, each with an exposure time of 20 s. Nodding along the slit was used during the September 2010 run, but only the "A" nod position was used here, as systematic errors were introduced into the RVs upon usage of different trace positions. Subsequent runs did not use nodding, as it was deemed unnecessary and detrimental to RV precision. Typical integration times were tuned to avoid nonlinear detector regimes (>2000 ADU per exposure) and to account for seeing variations, ranging from 60-200 s and 180-240 s for GJ 15 A and GJ 876, respectively, allowing for SNR per pixel ∼ 30 for each raw spectrum and ∼ 200 after ∼ 2 hours of combined integration. Seeing usually varied between 0. 5 and 1. 5. Exposure times of 5 s were sufficient to obtain SNR per pixel ∼ 100 for each spectrum of SV Peg under most seeing conditions. By comparison, the limiting magnitude of CSHELL in the K band (∼ 2.2 µm) for SNR per resolution element of 10, resolving power of 21500 (1. slit), a total observing time of 1 hour, and a single exposure integration time of 120 s is 12.8 (Greene et al. 1993). Converting from this standard to our case is difficult without a PSF model, which varies considerably for our observations. The time baseline of the observations analyzed in this work for GJ 15 A and GJ 876 are 817 and 707 days, respectively. For SV Peg, which is variable on a time scale ∼ 145 days with RV amplitude ∼ 1.5 km s −1 (Hinkle et al. 1997;Lebzelter & Hinkle 2002;Winnberg et al. 2008), we focus on intra-night and intra-run RV stability only.

Data Reduction
All spectra were extracted consistently with a custom interactive data language (IDL) pipeline. We found that our ability to achieve RV measurements at precisions < 55 m s −1 strongly relies on (1) a careful construction of the flat fields, (2) efficient detection and correction of bad pixels to prevent them from contaminating the extracted spectra, and (3) a correction of the instrumental fringing in the individual flat fields, which occurred at a ∼ 0.2-0.6% level and are caused by the CSHELL CVF. Instrumental fringing also affects individual science images; however, it is not possible to efficiently correct them directly due to low illumination of most of the detector. For this reason, residual fringing was taken into account in our spectral model when radial velocities are extracted from the reduced spectra (see §4).
Combined spectroscopic flat field images were created for every night by median-combining typically 15 individual flats with exposure times of 20 seconds, from which a dark frame was subtracted. Removal of the instrumental fringing pattern in the flat fields is shown in Figure  1. Fringing subtraction was done by first creating a fringing-free flat field (panel A) by averagecombining a large number of per-night flat fields; fringing gets averaged out since its phase and amplitude vary randomly across observing nights, leaving behind only the permanent detector response. This fringing-free flat field was subsequently divided from each combined flat field (one per night) to make 2D fringing patterns apparent (panel B). We median-combined individual columns of this image to obtain a 1D fringing pattern. A Levenberg-Marquardt least-squares algorithm was then used to fit the 1D fringing pattern with a shifted and renormalized interference function of the form is a normalization factor, A f is the amplitude of the interference pattern, F is the finesse parameter, P o is the pixel grid on the detector, ω p is the spatial frequency (in units of pixel −1 ), φ is the phase, and I 0 and s allow for a linear slope in the flux. The resulting parameters are then used as starting estimates to fit the full 2D fringing pattern using the following model : where y is the column number on the detector, perpendicular to P o , φ y allows for a linear phase shift with respect to rows (resulting in tilted interference fringes), A f,y allows for the amplitude to vary with rows and I 0 , s p , s y and s p,y allow for a 2D plane to be fit to the residual flux. The resulting 2D fringing pattern (panel C), from which the linear plane fit was subtracted, I 2 (P o , y) + 1 − (s p P o + s y y + s p,y P o y) was finally divided from the combined flat field to obtain a corrected flat field to be used in the data reduction process (panel D, multiplied by the fringing-free flat field). Though the fringing is not fully removed at the top and bottom edges of the image, the target spectral trace is always located at the center of the image, and so our technique is sufficient.
Individual science images were read separately and divided by the corrected flat field of the corresponding observing night. We cross-correlated a median spatial profile of the spectrum at each spectral position to determine a linear solution to the trace tilt and corrected it by interpolating in the spatial direction only. This step was repeated three times iteratively.
We re-constructed a median spatial profile of the straightened trace and fitted a Moffat profile, defined by A m /((P o − P o ) 2 /σ 2 m + 1) γ where A m is the amplitude, P o the central pixel position, σ m the characteristic width, and γ an index that controls the size of the wings (Moffat 1969). We used the resulting fitting parameters to set all pixels to zero at positions |P o −P o | > σ m / √ 2γ + 2.3 in the median spatial profile. The median spatial profile was then used in combination with a standard optimal extraction procedure (Horne 1986;Massey & Hanson 2013) to extract the spectrum a first time.
The resulting spectrum and spatial profile were then used to reconstruct a smooth synthetic 2D spectrum. We subtracted this synthetic spectrum from the raw science image, and created a map of small scale structures in the spatial direction by taking the minimum value between the absolute values of the upwards and downwards vertical (i.e., spatial direction) derivatives in the resulting difference image. We normalized each column of this map so that they have a unit median, and flagged all pixels with a resulting flux value larger than 5 as bad pixels. The optimal extraction was then repeated while forcing the masking of these previously identified bad pixels.
The extracted spectra were normalized to a unit continuum and no wavelength calibration was performed at this step. The instrumental blaze function was also not corrected in this step, and is instead treated as a free parameter in our modeling of the science spectra when RVs are computed (see §4). This allows a better determination of the blaze function, especially in the presence of several deep absorption features from the gas cell.
A final bad pixel detection filter was applied in the extracted 1D spectrum by identifying up to 5 pixels with flux values significantly larger than the continuum. We chose this value of 5 pixels as it is the maximal number of bad pixels that we observed in an extracted 1D spectrum. We then counted the number of pixels with flux values f i > t where t is a threshold value that goes from 1 to 1.5. If any range at least as large as ∆t = 0.05 could be identified over which the number of flagged pixels did not change and was lower than 5, then these pixels were flagged as bad and ignored in the RV extraction pipeline. Figure 2 shows typical reduced spectra for each of the three targets.

Radial Velocity Pipeline
The RV extraction pipeline is written in Matlab with a few bookkeeping scripts written in IDL used for determining barycentric corrections (barycentric vel.pro, Wright & Eastman (2014)). The RVs are retrieved by minimizing the difference between a model spectrum and the observed spectrum using the SIMPS Nelder-Mead amoeba simplex algorithm (Bajzer & Penzar 1999). Figure  3 shows a schematic of the RV pipeline and includes all the free parameters that define the model spectrum. We elaborate on these parameters and the calculation of the RVs in the following sections 2 .

Spectral Model
We base the model I obs (λ) on that of Butler et al. (1996), where λ is the wavelength solution of the model. However, additional components are necessary for our NIR observations, such that where I s is the stellar spectrum, which is derived iteratively (see §4.3); LSF is the line spread function/instrumental profile (P SF in Eq. 1 of Butler et al. (1996)); Σ is a sinusoidal function that compensates for the interference fringing left in I obs by the CVF filter; K is a quadratic blaze function that normalizes any curvature in the continuum, with the amplitude of each order (b 0 , b 1 , b 2 ) as free parameters; ∆λ s and ∆λ t are the Doppler shifts of the stellar and the telluric lines, respectively, and both are treated as free parameters; T t is the atmospheric telluric transmission function, with the line depths controlled by an exponent τ t where T to is the NOAO telluric absorption spectrum given by Livingston & Wallace (1991); and T g is the methane isotopologue gas cell transmission function with its line depths controlled in an identical manner: where T go was measured at high resolution using a Fourier Transform Spectrometer (FTS; Plavchan et al. 2013a). τ g is fixed at 0.97 under the assumption that the gas cell is stable over the observational period, with the specific value of 0.97 set through optimization tests. The fact that the fixed value is not 1.0 is likely due to the off-axis angle at which the cell was placed in the FTS, as opposed to in CSHELL, and the resulting slightly different light path lengths through the cell. τ t is allowed to vary as a proxy for the airmass, though this is only valid in our case since most of the telluric lines in our wavelength range are from methane absorption; if water vapor absorption were also present, then their lines would behave differently, and a single τ t would no longer suffice.
The upper panel of Figure 4 shows T go (bottom) and T to (top) in the wavelength range of the observations. This range was chosen to minimize the number of contaminating telluric and OH emission lines while maximizing the number of gas cell lines used for wavelength calibration and the number of stellar CO lines to increase the RV information content of the target spectra. T go and T to have 12 and 5 times higher resolution than that of our observations, respectively.
Even though T go is assumed to be much more accurately measured than T to and I s , we found features in the residuals of our fits of our spectral model to the observations that appeared to correspond to gas cell lines. To identify these features, we observed A type stars with the gas cell and fit to the data our spectral model using a flat line for I s ; the residuals were than averaged to reveal any coherent features, the majority of which should correspond to differences between the gas cell spectrum as observed at NASA IRTF and its laboratory-measured template, since the stellar spectrum has no features and there are only a few telluric lines. This process can be iterated multiple times by adding the averaged residuals to the gas cell and then repeating the fit and residual averaging. We expand on this process in §4.3, where we use it to derive the stellar template. The lower panel of Figure 4 shows the difference between the original T go and the final T go after multiple iterations. The lack of difference at the edges is due to the wavelength extent of the observations. There appears to be some correlation between the high frequency oscillations of the difference and the gas cell lines. The standard deviation of the difference is ∼ 1%, and may be attributable to the differences in the environmental properties between where the gas cell template was measured and where it is being used for observations of our M dwarf targets, though the gas cell temperature is kept steady to within 0.1 K when used for observations (Plavchan et al. 2013a). Other possibilities include the existence of low amplitude, coherent noise in the detector that was not accounted for during the data reduction process, as well as subtle interpolation errors in the RV pipeline. Regardless of the cause, using the new T go reduces our RV scatter compared to using the original T go .
T g is used to set the deviation of the gas cell spectrum, ∆λ g , from the "true" wavelength solution. In other words, we assume that the gas cell wavelength solution is correct, and any deviations of the gas cell lines in the observations from this true correct wavelength solution is due to a shift in the wavelength solution of the observations. We can then define a relative wavelength shift ∆λ = ∆λ s − ∆λ g that is directly related to the barycenter-corrected RV of the star with where c is the speed of light, λ c is the central wavelength of the spectral window, and v b is the barycenter velocity calculated using the IDL program barycentric vel.pro (Wright & Eastman 2014). We perform the stellar wavelength shift in logarithmic space to be consistent with the Doppler equation, as ∆ln λ ∼ ∆λ/λ c .
The LSF is constructed using Hermite functions ψ i (x) (i.e., Arfken et al. 2012), which are derived iteratively using the following recursive relation and zeroth and first degree terms: where x = P f /ω, P f is the pixel grid of the model, and ω is the standard deviation (width) of the Gaussian factor, which is a free parameter in the model. We compute the LSF by summing the first m terms where a i is the amplitude of the ith term relative to the 0th term and are also free parameters. Our nominal model uses m = 4, which we show to be sufficient in §5.1.4. The LSF is normalized to a unit area under curve, and the simplex minimization is constrained such that no part of the LSF becomes negative. We experimented with multi-Gaussian and multi-Lorentzian LSFs, as in previous works (e.g., Valenti et al. 1995;Butler et al. 1996;Johnson et al. 2006;Bean et al. 2007Bean et al. , 2010, but they yielded lower RV precision than our current Hermite function implementation, likely as a result of the relatively low SNR of our observations. The convolution is done on the pixel grid rather than the wavelength grid, though they are equivalent. As described in §3, sinusoidal interference patterns of variable amplitude and phase are present in both the flat fields and the science data due to the CSHELL CVF filter. These fringes become obvious in the spectrum of an A star taken without the methane isotopologue gas cell, as shown for 32 Peg in Figure 5. Fringes of a similar amplitude were also seen in the spectra taken by Brown et al. (2002) using NIRSPEC, and were corrected by parametrization of the sinusoid and linear regression. Meanwhile, Blake et al. (2010) also observed fringes while using NIRSPEC, which they corrected by applying a model sinusoid of fixed amplitude, period, and phase. We use a similar method and parametrize the fringing as a sinusoid function Σ(λ): where the parameters A, B, and C are free to vary.
The spectral model is compared to the data on a common wavelength grid. As with Crockett et al. (2011), we find that the wavelength solution of the data is variable between observations, and that the variability is quadratic in wavelength. We also find that it is dependent on the trace position y of the observed spectrum on the detector. Thus, we define the wavelength solution of the observations, λ o , with respect to λ by where α is a linear function of y and β is constant in y α and β are free parameters, P o is again the pixel grid on the detector, P c o is the central pixel of said pixel grid, and N pix = 256 is the total number of pixels in the data. The coefficients of Eqs. 13 and 14 are derived from fitting a linear trend to the relationship between y and a set of α and β values obtained from fits to several high SNR observations of GJ 537 A (Gagné et al. submitted to ApJ) where no y dependence was assumed (i.e. α = α, β = β) The model I obs (λ) is then interpolated onto λ o . We use 4096 pixels in the model pixel grid P f to reach sufficient resolution for convolution with the LSF, and downsample to 256 pixels when interpolating onto λ d . The downsampling is done by binning every 4096/256 = 16 model "pixels". Any model pixel that is lying partly in two data pixel bins is split between them according to the fraction of the model pixel that is in each of the two data pixel bins. Table 3 lists all of the free parameters that define the nominal spectral model, which are varied by SIMPS to optimize the fit between the model and the observed spectra. The optimization is accomplished by minimizing the average of the RMS and the robust sigma (Hoaglin et al. 1983;Beers et al. 1990) of the residuals after application of the bad pixel mask, which sets the weight of bad pixels to zero (see §3). Including the robust sigma in the minimization process decreases the impact of outliers due to noise, hot pixels, and unmasked bad pixels. However, just minimizing the robust sigma alone without the RMS produces bad fits to the telluric lines due to their similarity to outliers when there are very few of them, as in our case ( Figure 4). All parameters aside from the gas cell optical depth τ g are allowed to vary, while only the relative wavelength shift ∆ln λ, the gas cell wavelength shift ∆λ g , and the phase of the fringing correction C are allowed to vary freely. The other parameters are allowed to vary within bounds determined from optimization tests.

RV Calculation
Due to the low SNR per pixel (Table 2) and small wavelength range (∼ 6 nm) of our individual spectra, the RV precision of a target is calculated from nightly averaged RV values RV i , for the ith night, defined as the weighted mean of the RVs of all the individual spectra taken during that night, weighted by the inverse square of the RMS of residuals of each model fit to those spectra, ignoring the contributions to the RMS from bad pixels, where RV j i , w j i , I j d,i (k), and I j obs,i (k) are the RV, weight, observed flux of the kth pixel, and model flux of the kth pixel of the jth individual spectra on the ith night, respectively, with k bp indicating a bad pixel. The 1σ error bar of the ith nightly averaged RV point, δRV i , is calculated as the weighted standard deviation of the individual RVs, weighted by the same quantity as calculated for the nightly averaged RV values (see Eq. 16), and divided by the square root of the number of spectra N obs taken during that night, The RV precision of a target is then defined as the standard deviation of the nightly averaged RV values, while the reduced chi-square χ 2 red is calculated using where N n is the number of nightly averaged RV points (= number of nights of observations/epochs), and RV is the weighted mean of all the RV i values, weighted by the inverse square of the δRV i values.
Note that our method of calculating the RV precision of a target is different from that of Gagné et al. (submitted to ApJ), our companion survey paper. They define the RV precision as the weighted standard deviation of the nightly averaged RV values, weighted by δRV i , which allows them to better constrain planet mass sensitivity. In contrast, the purpose of the RV precision obtained in our work is to compare to the RV precision obtained by previous works, where it is sufficient to calculate just the standard deviation of the nightly averaged RV values. The method of Gagné et al. (submitted to ApJ) results in lower RV RMS overall, which is indicative of mixing data of different SNRs.

Stellar Template Generation
The retrieval of the original stellar spectrum, I s , for use in Eq. 3 has consistently been a difficult task. Some previous works have used synthetic stellar models that calculated stellar spectra, given effective temperatures and surface gravities (Blake et al. 2010;Crockett et al. 2011;Bailey et al. 2012;Tanner et al. 2012), but this runs the risk of introducing spectral features not present in the science targets into the spectral fitting, which lowers the RV precision. Alternatively, the stellar spectrum can be obtained from deconvolution of high resolution stellar observations using LSFs derived from observations of A and B stars taken with a gas cell (Butler et al. 1996;Bean et al. 2010;Rodler et al. 2012). However, this presumes that the LSF remains stable between the observations of the targets and those of the A and B stars, which is not true for our case due to temperature variations, instrumental flexure, and mechanical disturbances to CSHELL itself (e.g. moving the slit into position) as it moves with the telescope at its Cassegrain mount.
We thus use the target observations themselves to derive the stellar template iteratively using a method similar to that of Sato et al. (2002). However, our method differs from theirs in that the initial guess template is completely flat in order to minimize contamination of features that may not be present in the observed spectra, and that we use all of our target observations to derive the template in order to maximize SNR and sample a large range in barycenter velocities. This ensures that the stellar spectrum is decoupled from "stationary" features, such as the gas cell and telluric lines and the sinusoidal fringing, so that they do not contaminate the stellar template.
The procedure begins with the building of a model using the flat spectrum in place of I s in Eq. 3. The fit produces residuals similar in shape to the actual stellar spectrum, though they are distorted due to the simplex algorithm minimizing the residuals; for example, the stellar CO lines will be made more shallow, or raised above the continuum to reduce their impact on the RMS and robust sigma. We lessen this effect by repeating the first iteration with the CO lines masked out, allowing only the stellar continuum to contribute to the RMS and robust sigma. The CO lines are masked by ignoring all pixels with values one standard deviation below the mean of the normalized template flux derived from the original first iteration; two pixels on each side of each CO line are also masked to further reduce their impact on the fit.
The residuals are then "de-shifted" to the barycenter so that they all have the same Doppler shift, with the assumption that the star has zero RV perturbations. This is valid even for potential planet hosts, as RV perturbations caused by planets are much smaller than a single resolution element. Residual values that are three standard deviations above or six standard deviations below the mean residual value (zero), and those below −1, are set to zero so that they do not contribute to the eventual template (i.e. they are assumed to be bad pixels). We then take the median of these residuals, weighted by the inverse square of the average of the RMS and robust sigma of each fit, and add it to the template from the previous iteration to generate the new template that will be used in the next iteration. Pixels that are flagged by the bad pixel mask are given zero weight. In other words, residuals from better fits are more represented in the stellar template than those from worse fits. We find that variations in the quality of fit is usually due to the quality (SNR) of the data itself, and therefore it is reasonable to rely more on the best data to generate the stellar template.
This process is repeated until the RV precision of the target stabilizes, which usually takes ∼ 10 iterations, though we typically run the pipeline for 20 iterations to confirm that RV stabilization has been reached. We speed up this process by deconvolving the residuals from the first iteration with the best fit LSF corresponding to each residual using the Richardson-Lucy algorithm (Richardson 1972;Lucy 1974), so as to more quickly converge to the "real shape" of the stellar spectrum. Deconvolution is avoided in subsequent iterations, as it would significantly amplify the noise in the data.

Results and Discussion
5.1. RV Stability at Modest Signal-to-Noise: GJ 15 A 5.1.1. Radial Velocities Figure 6 shows example fits of the model spectrum to a high SNR per pixel (top) and a low SNR per pixel (bottom) observed spectrum of GJ 15 A. Note that the high SNR spectrum is not included in our RV calculations, as it was the only spectrum taken that night. It is clear that, at SNR ∼ 100, our spectral model is able to reproduce the observations with high fidelity, while at low SNR there is significant scatter and deviation between model and data, though all of the major features have been captured by the model. As Table 2 shows, the majority of our observations have SNR per pixel closer to ∼ 30-40, thus limiting our RV precision. We will quantify the effect of the SNR on our achievable RV stability in §5.1.3.
The iterative nature of our RV pipeline results in multiple values of RV precision for a single target, one for each iteration, and we accept the lowest RV scatter (RMS) value among all the iterations as the RV precision achievable by our RV pipeline for that target. Figure 7 shows the RV scatter as a function of iteration for GJ 15 A. The first few iterations with high RV scatter result from errors in the stellar template, as it is continually augmented from one iteration to the next. Convergence in both the stellar template and the RVs is reached beyond iteration 9, after which the RV scatter is ∼ 40 m s −1 , with occasional deviations to lower values. Figure 8 shows the RVs from iteration 13, where we are able to achieve a RV precision of 35 m s −1 over a 817 day long baseline. Table 4 gives the nightly RV values and associated 1σ error bars. The obtained RV precision is consistent with the theoretical RV precision given our observational setup and a SNR per pixel of ∼ 100, the effective SNR of each of our nightly averaged RV points (see Table 2 of Anglada-Escudé et al. 2012). If the individual RV points are considered rather than the nightly RVs, then the RV precision is 131 m s −1 , and the SNR of the observations are those between the minimum and maximum SNR values listed in Table 2. Figure 9 shows the progression of the stellar template with increasing iteration. In the initial iterations, the depth of several of the CO lines are much shallower than the central CO line, when they should be similar in depth instead. This is due to some overlap between these stellar lines and the gas cell lines, allowing the gas cell template to partially compensate for them during the spectral fits. However, at higher iterations, improved fits to the observations arising from improved stellar templates lead to similar CO line depths. The higher iterations also show more noise in the templates, which is a drawback of this algorithm, as any coherent noise in the residuals will be added to the template. Therefore, at higher iterations continuous augmentation of the stellar template does not lead to higher RV precision, as shown in Figure 7.

Parameter Correlations
Given our large number of nuisance parameters (Table 3), it is important to investigate any correlations between themselves and between them and the RVs. Correlations are tested by evaluating Pearson's linear correlation coefficient ρ, defined by (Pearson 1895 where cov(X, Y ) is the covariance of the variables X and Y , and σ X and σ Y are their standard deviations. A large degree of correlation corresponds to |ρ| → 1. We calculate ρ by assigning the RVs and different spectral model parameters to X and Y . Figure 10 shows the relationship between every parameter and the individual RVs of GJ 15 A. No obvious trends can be discerned for any of the parameters, while the absolute value of ρ is < 0.4 for all of the RV-parameter pairings. Some parameters have hit their parameter bounds, such as the parameters controlling the shape of the LSF. For example, the amplitude of the first degree Hermite function (a 1 ) hits the bounds in both directions; this is understandable since it shifts the LSF back and forth, which could create a false RV signal, and thus it is necessary to constrain it to a narrow set of values. The FWHM of the LSF (2ω √ 2Ln2) spans the critical sampling resolution at the precision of CSHELL (∼ 2 pixels/resolution element), and therefore we are both under and oversampling our science spectra.
Clustering of points can be seen in some of the panels. For ∆ln λ, the clusters correspond to different epochs with different barycentric corrections (see Eq. 6 and the discussions that follow). For ∆λ g and ∆λ t , the clustering corresponds to different epochs where the central wavelengths were set to slightly different values. Variations in the LSF parameters are due to the variability inherent in observing using a non-stabilized spectrograph. Figure 11 shows the variability in the LSF within a single night (top) and between several nights (bottom). To quantify the variability, we take the average of the standard deviation of the difference between the LSFs and their mean within 4 pixels of the x-grid zero point. The resulting σ LSF values show that the LSF is about 4 times more variable between nights than within a night, which could account for some of the systematic RV differences between nights.
While no parameters show unexplained increasing or decreasing trends with RV, several parameters do have ranges in values where the corresponding RV scatter is greater than that of other values of said parameters. For example, low values of the sinusoidal fringing amplitude (A) appear to correspond to RV points with lower scatter than those that correspond to high values of A. Other examples of this phenomenon include the blaze function polynomial coefficients (b 0 , b 1 , b 2 ) and the wavelength solution parameters (α, β). These trends in RV scatter are shown more clearly in Figure 12, where we bin every 10 RV points in the order of increasing value of the chosen parameters and calculate their RV scatter. We thus see that RV precision appears to increase with increasing continuum level, decreasing (in magnitude) linear and quadratic terms in the blaze function, decreasing sinusoidal fringing amplitude, and decreasing (in magnitude) linear and quadratic terms in the wavelength solution.
The absolute value of ρ is < 0.5 for most of the parameter-parameter pairings. The largest |ρ| (0.9985) is obtained between the gas cell wavelength shift ∆λ g and the telluric wavelength shift ∆λ t , which is expected since they both act as wavelength calibration, even though the gas cell is a more reliable fiducial. Aside from this correlation, the largest |ρ| values ( 0.6) occur between the six parameters shown in Figure 12, and their correlations and ρ values are shown in Figure 13.
For the correlations between the linear and quadratic terms of the blaze function (b 1 and b 2 ) and wavelength solution (α and β), anti-correlations may arise due to each pair being alternating orders in a polynomial. In addition, the change in shape of the continuum slope may play a role as well. Consider a linear slope in the continuum of a spectrum; if the data wavelength solution is expanded in both directions equally, corresponding to a more positive α, then the linear slope of the data will tend towards zero, which, since most of the linear slopes in the sample are negative, results in a more positive model linear slope b 1 . This explains the positive correlation (large positive ρ) between α and b 1 . Similarly, a positive β corresponds to an expansion at longer wavelengths and compression at shorter wavelengths; this results in the generation of positive curvature when the linear slope is negative, corresponding to a positive b 2 , and thus the positive correlation between b 2 and β. We will see in the next section that some of the remaining correlations, such as those between the constant blaze function term b 0 and the fringing amplitude A are related to the SNR of the spectra.

Error Analysis
Typical contributions to the RV scatter of a target include photon noise, wavelength calibration errors, instrumental effects, and astrophysical sources, such as stellar activity (Bailey et al. 2012). GJ 15 A has been shown to be stable down to 3.2 m s −1 (Howard et al. 2014), with the dispersion dominated by the RV signal from a Super-Earth exoplanet (K = 2.94 m s −1 ); this can be regarded as the total error contribution from astrophysical sources, and it is clearly far lower than our calculated RV precision. Therefore, our RV error budget is dominated by the other terms.
We now calculate the error contribution from photon noise, σ RV , due to the RV information content of the stellar spectrum and the SNR of the observations. Following Butler et al. (1996), where dI i /dV i is the slope of the stellar spectrum at the ith pixel, calculated as the change in the normalized flux of the spectrum at the ith pixel, I i , divided by the corresponding change in wavelength, expressed in units of velocity, at the ith pixel; i is the uncertainty in the flux at the ith pixel, given by in the photon-limited case, where N pho is the number of photons and the SNR is that of the observation (Table 2).
Due to our much smaller wavelength range compared to that of Butler et al. (1996), we use Eqs. 20-21 differently. In their work, they split their one observed stellar spectrum of τ Ceti into 704 2Å segments and calculated σ RV for each segment (i.e. the summation in Eq. 20 is over the pixels of a single segment); they then compiled a histogram of σ RV values for the segments, and calculated the RV precision arising from combining all the segments. Our small wavelength range forbids us from splitting our observations into multiple segments, so instead we treat each observation as a single "segment", giving us 103 segment in total for GJ 15 A; we then calculate σ RV for each observation/segment by summing over every pixel in the observation in Eq. 20 and using the SNR per pixel of each observation in Eq. 21.
An additional complication is that we do not have any observations of the star by itself, as even observations taken without the gas cell are contaminated by atmospheric telluric lines. Therefore, we build a synthetic stellar spectrum for each observation using the stellar template for the iteration with the best RV precision. The synthetic spectrum is calculated by (1) downsampling the stellar template to the data wavelength solution derived from the spectral fit of the model to that particular observation, (2) convolving the result with the LSF derived from the same fit, and (3) adding onto the result the residuals from the fit to simulate the appropriate level of noise.
The upper panel of Figure 14 gives a histogram of the calculated σ RV values of the 103 individual observations taken for GJ 15 A. The most common σ RV value is ∼ 120 m s −1 , slightly lower than 131 m s −1 , the actual RV precision for the individual observations (not the nightly coadded RV points) as output by our RV pipeline.
We can repeat steps (1) and (2) above for the gas cell spectrum to find the RV scatter contribution from errors in the wavelength calibration. We omit step (3) as the gas cell does not contribute to the photon noise. The lower panel of Figure 14 shows the result. Given the higher line density of the gas cell spectrum it is not surprising that its σ RV values are much lower than that of the stellar spectrum. The peak σ RV value is ∼ 50 m s −1 . Combining this with the photon-limited errors calculated above in quadrature gives √ 120 2 + 50 2 = 130 m s −1 , essentially equal to the RV precision of the actual individual observations. We can thus conclude that our RV scatter is dominated by photon noise, and to a lesser extent from errors in wavelength calibration.
We can also estimate the theoretical nightly averaged RV scatter by assuming the 103 individual observations are spread out evenly among the 14 nights, such that the RV scatter should decrease by 103/14 ∼ 2.7, so that it equals 131/2.7 ∼ 48 m s −1 , slightly higher than our actual results (Figure 7). This is likely due to the fact that the observations are not spread out evenly among the nights. Figure 14 shows a small local maximum at higher values of σ RV , which correspond to the low SNR observations of the night of July 13th, 2011 ( Table 2). As the RV content and wavelength calibration of each observation are nearly identical, the chief factor causing the higher σ RV values is the low SNR. The upper panel of Figure 15 shows the relationship between RV and SNR per pixel. There is a clear trend showing that RVs of low SNR observations have greater scatter than RVs of high SNR observations. We can thus estimate the impact the SNR has on the RV precision by repeating the calculations of Figure 12, but with the parameter values replaced by SNR values; the result of this calculation is shown in the lower panel of Figure 15. To quantify the relationship between RV scatter and SNR, we fit a power law to the points with the form where σ 0 is the error contribution from all other sources and is a free parameter in the fit, and s 1 and s 2 are constants that are also free parameters. The power law is consistent with a 1/SNR relationship between σ RV and SNR, as expected when the main source of error is due to photon noise.
Given the impact the SNR has on the RV scatter, it is useful to investigate how SNR variations affect the free parameters in the spectral model, which is shown in Figure 16. As with Figure 10, most of the parameters show no clear trends with SNR, except for the blaze function coefficients, fringing sinusoid amplitude, and the linear and quadratic terms in the wavelength solution. It is understandable why these parameters are the most affected. At low SNR, neither the continuum nor the wavelength solution are well constrained due to the noise in the spectra, leading to large deviations from zero (or 1.05, in the case of the constant term in the blaze function). The rise in fringing amplitude with decreasing SNR can be seen as either an actual increase in the relative amplitude of the fringing as the number of photons from the star decreases, or as the model attempting to fit the higher noise level with a sinusoid in low SNR spectra.
These results help to explain some of the correlations between the parameters seen in Figure 13. For example, as the blaze function constant term b 0 can only decrease and the fringing amplitude A can only increase at low SNR, they appear to correlate with each other, whereas the true cause of the correlation is the low SNR. Similar arguments can be made for correlations between these two parameters and the other four parameters of Figure 13.
Sources of error stemming from instrumental effects are likely dominated by contributions from bad pixels, which were abundant in the science images and can be generally defined as due to detector imperfections and cosmic rays. Figure 17 shows the relationship between RV scatter and the number of bad pixels N bp flagged on the 2D trace of GJ 15 A spectra. As with Figure 15, the upper panel plots the individual RV values against N bp , while the lower panel shows the RV scatter for every 10 points binned in increasing N bp . We once again fit a power law to the binned points, though we ignore all points corresponding to N bp > 150 since our knowledge about the RV precision for this range of N bp is limited by the small number of spectra. The fitted curve is the same as in Eq. 22, with N bp in place of the SNR. The initial steep rise in σ RV from N bp = 0 to N bp = 150 shows the significant impact bad pixels have on both the data reduction and retrieval of RVs. Beyond N bp = 150, the RV scatter appears to saturate. Bad pixels not only cause certain parts of the spectrum to be unusable, thereby decreasing the information content of the data, but they also affect the data reduction process by interfering with the initial fitting to the spatial profile of the 2D trace.

Exploring the Multiplicity of LSF and Blaze Function Terms
The LSF and blaze function components of our spectral model are constructed from the summation of multiple Hermite function and polynomial terms, respectively. However, there are no significant a priori constraints on how many terms these components should have. We thus perform sensitivity tests by running our RV pipeline using data from GJ 15 A and various numbers of component terms in the LSF and blaze function and then comparing the resulting RV precision. Figure 18 shows the resulting RV precision when using LSFs constructed with up to 11 Hermite function terms in the spectral model. Since each run of the RV pipeline involved 20 iterations, a spread in RV RMS values is generated. Accordingly, we plot each run as a box plot, where for each box the error bars show the full spread in RV RMS values; the vertical extent of the box covers the upper to lower quartile of RV RMS values; and the red line indicates the median RV RMS value. For these tests the blaze function is set as a quadratic polynomial. It is clear that an LSF generated by summing up the first 5 degrees of Hermite functions correspond to the lowest median and absolute RV RMS. This can be understood as the balancing of two effects: for LSFs made up of fewer Hermite functions, the degree of freedom is too low and the LSF is not able to account for all distortions in the line profile. Conversely, LSFs made up of more Hermite functions run the risk of fitting noise in addition to the line profile, diluting the RV information content of the spectral lines. These results are in contrast to previous works that used single Gaussians for their LSFs (Crockett et al. 2011;Blake et al. 2010;Bailey et al. 2012;Tanner et al. 2012), which corresponds to the 0th degree Hermite function. However, a direct comparison with the latter three works may be inappropriate due to their reliance on a different spectrograph with a reduced spectral resolution. Meanwhile, Crockett et al. (2011) also used CSHELL, and found no improvements to RV precision by using multi-Gaussian LSFs. It is unknown how their results would change if Hermite function LSFs were used instead. An additional complication is the change in LSF shape with wavelength, which was included in the model of Blake et al. (2010). However, this is not an issue for our work due to our short wavelength range (6 nm) compared to theirs (34 nm), such that we can assume the LSF shape is independent of wavelength. Figure 19 shows the resulting RV precision when using blaze functions constructed with polynomials of order 6 or less in the spectral model, in a similar format as for the LSF tests. The test case with a constant blaze function ("order 0") is not shown as its RV RMS value is more than 5 times larger than the plotted cases. For these tests the LSF was set as the sum of the first 5 degrees of Hermite functions. Unlike the LSF tests, the best solution is not immediately clear, though the quadratic blaze function provides the lowest minimum and median RV scatter. Blaze functions with order >4 increase the RV scatter due to the high sensitivity of the blaze function amplitude to the higher order terms. Previous works have used both linear (Valenti et al. 1995;Bailey et al. 2012) and quadratic (Crockett et al. 2011) blaze functions, with the choice likely depending on the degree to which the spectra were already normalized to the continuum during the data reduction phase. For example, the California Planet Search RV pipeline uses three order-6 polynomials of decreasing amplitude to normalize their spectra to the continuum before the RV extraction phase, during which only a constant normalization factor is needed (John A. Johnson, 2013, private communications). This method is valid for optical wavelengths where the continuum of the spectrum is well defined between the spectral lines. However, as the spectral lines dominate over the continuum in our wavelength range (Figure 2), higher order blaze functions could have exaggerated amplitudes due to fitting to deep lines or absorption bands, thereby introducing distortions to the spectral model flux when using them for normalization.

RV Stability at High Signal-to-Noise: SV Peg
The ensemble RVs of SV Peg and their associated 1σ error bars are given in Table 4 and are shown in Figure 20. RV variations with amplitude ∼ 1.5 km s −1 spanning hundreds of days can be seen, which are consistent with the 145 day variability time scale and amplitude from previous observations (e.g., Hinkle et al. 1997) and reaffirms SV Peg's quasi-periodic pulsations as an M giant. Therefore, in order to investigate RV stability using SV Peg observations, we consider its RVs on a shorter time scale. For example, if we assume purely periodic variability, such that the RVs vary by 3 km s −1 within half of the total variability period, 145 days, then within 1 hour the RVs will change by only ∼ 1.7 m s −1 . The actual change is likely to vary by at least a factor of two due to changes in the RV slope occurring on shorter time scales (see, for example, Figure 8 of Hinkle et al. 1997). Figure 21 shows the individual RVs of SV Peg on five nights when > 100 spectra were taken, each with RV RMS ∼ 18-40 m s −1 . As these values are an order of magnitude larger than the predicted RV variations within a few hours, it is not surprising that no long term trends are apparent in the data. The left column of Figure 22 shows the RV precision as a function of increasing SNR per pixel for each observation, calculated from progressively binning an increasing number of individual RV points, up to one-third of all points obtained that night per bin. The resulting RV noise floors lie at ∼2-6 m s −1 . The SNRs needed to reach the noise floor are ∼400-600 for the 5 nights, which reflects differing weather conditions, pointing/guiding/focusing issues, detector artifacts, and contributions from read noise and dark noise. The theoretical integration time to reach the noise floor for SV Peg, calculated as the product of the integration time per individual exposure (∼3-5s) and the square of the ratio of the total SNR required (∼400-600) to the SNR per individual exposure (see Table 2), is ∼200 s. Both the individual and binned RV RMS values are consistent with previous analyses of these data performed by Anglada-Escudé et al. (2012), as well as the level of precision expected given SV Peg's intrinsic RV variability within hours, as calculated in the previous paragraph. In other words, as a result of the high SNR of the observations and the high information content of the deep CO lines of SV Peg stellar spectra, our RV precision for SV Peg may be significantly limited by stellar variability in addition to photon noise or wavelength calibration.
To test this possibility we repeat the calculations shown in the left column of Figure 22 and add an extra step where we subtract a linear trend from the binned RV points at each SNR before calculating σ RV . The linear trend is fitted to the binned RV points weighted by the inverse square of the scatter in each bin and represents the short-term RV variability due to stellar pulsations. The results are shown in the right column of Figure 22, where for all but the night of 11/08/20 there is some reduction in the RV slope. However, this is not definitive evidence that we have uncovered the signature of stellar pulsations, as the subtracted linear trend could be due to some other systematic, such as variations in the gravity vector, focus, and/or temperature. Figure 23 shows the linear RV slope for the subtracted linear trend at each SNR. If the slope is due to pulsation, then they should all have similar magnitudes (a few m s −1 hr −1 ) and be consistent across many values of the SNR. It is clear that this is not the case for most of the nights, where the magnitudes are either too large (10/10/09, 10/10/10, 10/10/13) and/or too variable (10/10/10, 11/08/20). By contrast, for the night of 11/07/10, the slope's magnitude is consistent with stellar pulsations and is stable across many SNR values. In addition, it features the largest decrease in noise floor upon subtraction of the linear trend, and the total duration of the observations (4 hours) may be sufficient for RV variability due to stellar pulsations to be made apparent. However, more precise RV measurements and complementary observations of activity indicators are needed to fully evaluate whether the linear slope is indeed due to stellar variability. Figure 24 shows the correlation between RV and the RV pipeline parameters for SV Peg observations taken on August 20th, 2011. This date was chosen as it had the most RV points. These results show much less structure than those shown in Figure 10, which is due to the much shorter time baseline of the observations. Like for GJ 15 A, no obvious correlations exist, except for the relative wavelength shift, ∆ln λ. This is due to scatter in the retrieved RVs of SV Peg. To illustrate, suppose the retrieved RVs are made up of two components: a scattered component, and an ideal component that, when added to the barycentric correction, results in zero RVs. Therefore, the barycenter-corrected RVs would just be made up of this scattered component, and thus they are linearly correlated with each other. The same correlation exists in the analogous panel of Figure  10, but is isolated to single epochs. Figures 25 and 26 show the nightly RVs for two separate observing runs, one in October 2010 and one in November 2010, respectively. Each run consists of spectra taken on nearly consecutive nights, revealing a short-term linear RV trend resulting from SV Peg's intrinsic variability. Given the ∼ 1.7 m s −1 per hour RV slope calculated above, over several days we expect to observe RV changes of ∼ 100 m s −1 , roughly consistent with the slopes shown in Figures 25 and 26. Once the linear trends are subtracted, the resulting RV scatter is ∼ 3-4 m s −1 , similar to the noise floors of the intra-night results. The detrended nightly RVs of the Oct. 2010 and Nov. 2010 runs and their associated 1σ error bars are given in Table 5. Both our intra-night and intra-run results are consistent with the observed RV stability of other M giant stars on similar time scales (e.g., Seifahrt & Käufl 2008).
The error bars for the SV Peg nightly RVs shown in Table 5 are smaller than those of the same nightly RVs shown in Table 4. This results from the different sets of spectra input into the RV extraction pipeline for these two cases, and the way we derive our stellar templates. In Table  4, all SV Peg spectra were input simultaneously, deriving one single stellar spectrum used to fit all observed spectra despite SV Peg's variability; in Table 5, only the spectra from the individual runs were input, and separately, thus producing two different stellar spectra for the two runs. The smaller error bars for the nightly RVs in Table 5 then suggest that the nightly RV scatter is smaller when using a stellar template derived specifically from the spectra within those runs, rather than one derived using all available SV Peg spectra. In other words, the stellar spectrum of SV Peg is changing over its long-term variability cycle such that no one stellar template can be used to satisfactorily fit all observations. We explore this effect in Figure 27, where five stellar templates derived from the five separately analyzed nights shown in Figures 21 and 22 are compared. Some small differences are expected due to noise in the observations, but larger differences are also present, as marked by the arrows. The wavelengths of these large differences do not correspond to any major gas cell or telluric features and thus are likely due to changes in SV Peg itself. We also do not see this phenomenon in any of the other targets analyzed using our RV extraction pipeline (Gagné et al. submitted to ApJ). Furthermore, the two groups of templates, one derived from 2010 data and one derived from 2011 data, show greater differences between them than between templates within each group, indicating, as expected, that the stellar spectrum changed more over months than over days. Although a full interpretation of these results is beyond the scope of this paper, it would be interesting to ascertain what these changes to the stellar spectrum mean for the chemistry and dynamics of the atmosphere of SV Peg.

Validation of Planet-Detection Capabilities: GJ 876
The top panel of Figure 28 shows our GJ 876 nightly RVs (blue), while Table 4 gives their values and their associated 1σ error bars. The amplitude of the RV variations are consistent with those observed in optical surveys, aside from the outlier at JD-2455455.0 ∼ 300 with the large error bars (Delfosse et al. 1998;Marcy et al. 1998Marcy et al. , 2001Rivera et al. 2005Rivera et al. , 2010. Using the Systemic Console (Meschiari et al. 2009), Gagné et al. (submitted to ApJ) was able to fit a 1-planet solution to the data without any prior constraints and retrieve orbital parameters within 2σ of the published values for planet b, thereby confirming, for the first time, the existence of GJ 876 b using NIR RVs.
To further test the pipeline results, we fit a 2-planet solution to our data using the Systemic Console, corresponding to planets b and c. The best fit RV curve is shown in the top panel of Figure 28 (dashed line). Due to the low cadence of our data and the rapid dynamical evolution of the system, it is difficult to accurately retrieve the orbital parameters of both of these planets. Instead, we fix the planet masses, orbital periods, and orbital eccentricities to the published values and allow the mean anomaly, longitude of periastron, and the velocity zero point to vary. Our best fit solution does not show the alignment of the two planets' longitudes of periastron reported in previous studies (e.g., Laughlin et al. 2005;Rivera et al. 2005), but that is likely due to the low cadence of our observations. The fit results in a residual RMS of 68 m s −1 (bottom panel of Figure  28); removing the outlier at JD-2455455.0 ∼ 300 reduces this value to 35 m s −1 , within the RV precision of our RV standard, GJ 15 A. We do not pursue 3-planet and 4-planet solutions, as the two lower mass planets in the system are beyond the detection limits of our pipeline.
These results demonstrate that our RV pipeline is capable of detecting Jupiter-mass planets around M dwarfs with orbital periods of tens of days using CSHELL. However, it is clear that a confirmed detection will require high cadence, dedicated observations of RV variables in order to better constrain the orbital parameters of any possible planets they may host.

Summary and Outlook
We have constructed a data analysis pipeline that can process into RVs the NIR spectra of M dwarfs observed using CSHELL on NASA IRTF with the aid of a methane isotopologue gas cell by optimizing fits between the observed spectra and a spectral model. The pipeline is able to retrieve the stellar spectrum in an iterative process directly from the science observations without the need for additional observations, and takes into account temporal variations in the instrumental LSF, curvature in the continuum (blaze function), sinusoidal fringing, and wavelength solution. The pipeline is capable of (1) obtaining an RV precision ∼ 35 m s −1 for the RV standard M dwarf GJ 15 A, with the error contributions coming mostly from photon noise, (2) obtaining, for the high SNR target, SV Peg, noise floors of ∼ 2-6 m s −1 on individual nights, and ∼ 3 m s −1 RV precision over several days after subtracting out a linear trend consistent with RV variability due to stellar pulsations, and (3) detecting/confirming at least one Jupiter mass planet around the M dwarf GJ 876 in the NIR. These results are derived from observations with a 22-year old InSb detector within CSHELL, a spectrograph that was not designed with precision NIR RVs in mind, and which is mounted at the Cassegrain focus, thereby introducing significant mechanical, thermal, and pressure variations that all contribute to the RV scatter. Furthermore, the RV information content is restricted to five CO lines across a single order of 256 pixels that covers only 6 nm (Greene et al. 1993). However, despite these limitations, our results compare favorably to previous works that use larger telescopes, newer detectors, and spectrographs with greater mechanical, thermal, and pressure stability. Therefore, it will be useful to evaluate the possible improvements to our RV stability resulting from the reduction or elimination of the aforementioned limitations in newer NIR spectrographs. We will base our discussions on iSHELL, as it is due to replace CSHELL in the near future. iSHELL is a cross-dispersed high resolution echelle spectrograph with R ∼ 70000 at minimum slit width (0. 375) (Rayner et al. 2012). As with CSHELL, it will be fixed to the Cassegrain focus of the 3 m telescope at NASA IRTF. It uses a Teledyne 2048x2048 Hawaii 2RG array as its main detector and enables multiple orders to be observed simultaneously. A single exposure in the K band allows for ∼ 200 nm to be captured, thereby improving our wavelength coverage by a factor of 200/6 = 33; extrapolating from our calculations in §5.1.3 with the assumption that the RV information content stays roughly constant with wavelength, this results in a typical σ RV = 130/ √ 33 ∼ 22.5 m s −1 for the individual RVs due to photon noise and wavelength calibration errors. Given an SNR per pixel of 50 for each observation, and a total SNR of 200 per night, the nightly RV scatter would thus be 22.5/(200/50) ∼ 5.6 m s −1 . This is consistent with the estimates of Reiners et al. (2010) and Bottom et al. (2013) and will facilitate the detection of Super-Earths and Mini-Neptunes in the Habitable Zones of M dwarfs.
On the other hand, while the higher spectral resolution of iSHELL will increase the RV information content of the observed spectra, the SNR per pixel will decrease at the same time due to spreading out the same number of photons across more pixels. The impact of this effect on the RV precision depends on the sharpness of the spectral features of the target, which is a function of its rotational velocity. A target that rotates quickly (vsini ∼ 10 km s −1 ) may be insensitive to a change of R from 46000 to 70000, while a target that rotates slowly (vsini < 1 km s −1 ) may see its RV scatter increase by 50% (Bouchy et al. 2001). In addition, the increase in wavelength coverage will cause severe contamination of the observed spectra by atmospheric water and methane telluric absorption lines and OH emission lines, which are outside of the 6 nm wavelength range of our current CSHELL observations. Bean et al. (2010) corrected for telluric contamination by modeling the absorption and emission lines using line lists from the HITRAN database in conjunction with a time-resolved model of the atmosphere above the observatory and the LBLRTM radiative transfer code. Further improvements to the modeled telluric lines were made by fitting to a telluric standard star in order to determine systematic wavelength shifts and strength variations of the lines compared to their default values in the HITRAN database. A similar method can be applied to iSHELL observations, with additional improvements to the synthetic telluric template made possible by iterative processes akin to what we have done to improve the gas cell template ( §4.1), though it will be more challenging in this case due to the presence of multiple components (methane and water absorption lines and OH emission lines). Alternatively, a synthetic telluric spectrum can be constructed from principle component analysis using a library of observations of telluric standard stars, which was shown by Artigau et al. (2014) to improve RV precision in the R band without the addition of spurious RV signals, though it remains to be seen how well this method works in the H and K bands. Finally, the deepest telluric lines can simply be masked during RV extraction, though this reduces the RV information content of the observations. There exist several factors in addition to the increased wavelength coverage, however, that will help increase the RV precision of targets observed using iSHELL. The upgrades to the filter selection wheel will largely eliminate the sinusoidal fringing and thus its contributions to the RV error. The newer detector will lead to a decrease in the number of bad pixels on the 2D trace, further improving the RV precision by increasing the available information content and reducing the need to interpolate along the spatial direction during data reduction. The increased efficiency of the detector will allow for higher SNR for the same targets currently observable using CSHELL, while also making it possible to observe dimmer targets. Recent tests of iSHELL showed that an SNR per pixel of 100 for targets with K magnitude ∼ 9.5 is achievable within one hour (Peter P. Plavchan, 2015, private communications). Comparisons with the targets of Gagné et al. (submitted to ApJ) show that this will permit future RV surveys with NASA IRTF to observe targets 2-3 K magnitudes dimmer, thereby increasing the sample size of potential RV targets. Additional upgrades to iSHELL, such as the implementation of fiber feeds for improved stabilization of detector illumination and a laser comb for higher-fidelity wavelength calibration could be possible, and have already been demonstrated using CSHELL (Plavchan et al. 2013b;Yi et al. 2015). Given the proliferation of precision NIR RV surveys in the near future, a survey undertaken using iSHELL will nicely complement these efforts.
We thank K. Sung, S. Crawford, B. Drouin, E. Garcia-Berrios, N. S. Lewis, S. Mills, and S. Lin for their effort in the building and setting up of the methane isotopologue gas cell. We thank B. Walp for his help with data collection at NASA IRTF. We thank J. Rayner, L. Bergknut, B. Bus, and Table 3 for the definitions of the parameter symbols.  Table 3 for the definitions of the parameter symbols. Blue dotted lines mark the zero values in the b 1 , b 2 , α, and β plots.    Normalized Flux 10/10/09 10/10/10 10/10/13 11/07/10 11/08/20