We present the results of a precise near-infrared (NIR) radial velocity (RV) survey of 32 low-mass stars with spectral types K2-M4 using CSHELL at the NASA IRTF in the $K$-band with an isotopologue methane gas cell to achieve wavelength calibration and a novel iterative RV extraction method. We surveyed 14 members of young ($\approx$ 25-150 Myr) moving groups, the young field star $\varepsilon$ Eridani as well as 18 nearby ($<$ 25 pc) low-mass stars and achieved typical single-measurement precisions of 8-15 m s$^{-1}$ with a long-term stability of 15-50 m s$^{-1}$. We obtain the best NIR RV constraints to date on 27 targets in our sample, 19 of which were never followed by high-precision RV surveys. Our results indicate that very active stars can display long-term RV variations as low as $\sim$ 25-50 m s$^{-1}$ at $\approx$ 2.3125 $\mu$m, thus constraining the effect of jitter at these wavelengths. We provide the first multi-wavelength confirmation of GJ 876 bc and independently retrieve orbital parameters consistent with previous studies. We recovered RV variability for HD 160934 AB and GJ 725 AB that are consistent with their known binary orbits, and nine other targets are candidate RV variables with a statistical significance of 3-5$\sigma$. Our method combined with the new iSHELL spectrograph will yield long-term RV precisions of $\lesssim$ 5 m s$^{-1}$ in the NIR, which will allow the detection of Super-Earths near the habitable zone of mid-M dwarfs.

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.