Self-consistent Stellar Radial Velocities from LAMOST Medium-resolution Survey DR7

The scientific goals of the LAMOST MRS mainly include stellar multiplicity, stellar pulsation, star formation, emission nebulae, Galactic archeology, host stars of exoplanets, and open clusters. The details of the scientific plan and the survey strategy are described in Liu et al. (2020). In this section, we summarize them briefly.

Each scientific goal has a principal investigator who is responsible for its targeting. A planid is assigned to each MRS plate (pointing), which has the form [TD/NT]hhmmss[N/S]ddmmssXnn. The first two characters denote time-domain (TD, will be repeatedly observed during the five-year survey) or non-time-domain (NT, just observed on a single night), hhmmss[N/S]ddmmss represents the equatorial coordinates of the plate center, [N/S] means north/south, the digit X is used to denote the scientific goal (see Table 1), and the last two digits represent a serial number. For example, the planid TD062610N184524B01 means the plate is a time-domain plate dedicated to binarity/multiplicity research. There also exist some irregular planid that are testing fields, such as HIP8426401 and NGC216801. HIP8426401 means the central star of the plate is HIP 84264, and NGC216801 means a plate toward the star cluster NGC 2168.

The MRS uses the local modified Julian minute (LMJM), an 8 bit integer defined as 1440× the local modified Julian date at the beginning of each exposure, as the stamp of each exposure. The typical exposure time is set to 1200 s, while 900 s and 600 s exposures also exist, depending on the brightness of the targets. Each NT plate is observed with three consecutive exposures while each TD plate is observed until it is unobservable (usually five to six 1200 s exposures are allowed). An arc lamp exposure is taken at the beginning of each plate, and at the end of an observing night (every ∼2 hr). The targeting is mostly based on the Gaia DR2 source catalog (Gaia Collaboration et al. 2018). For a 1200 s exposure, the magnitudes corresponding to S/N = 5 in blue- and red-arm spectra are G ∼ 14.5 and 15.2 mag, respectively (see Figure 1 in Liu et al. 2020). Magnitudes brighter than these cover the majority of the objects observed by the LAMOST MRS.

The whole LAMOST MRS DR7 catalog (R ∼ 7500) contains over five million single-exposure spectra of over 800,000 stars obtained from 2017 September to 2019 June. In this work, 3,753,659 spectra with S/N > 5 either in the blue arm or in the red arm for 600,771 stars are selected. Their distributions on the number of exposures and time span are presented in Figure 1. The spectra are oversampled. The typical wavelength steps are 0.11 and 0.14 Å at the blue and red arms, respectively. The sampling rate (λ/Δλ ∼ 45,000) is as high as six times the spectral resolution (R ∼ 7500). In Figure 2, we show a spectrum of a K-type giant star as a demo of the LAMOST MRS. The blue arm and red arm are designed mainly for the Mg i b triplet at around 5175 Å and for Hα at around 6564 Å. The released spectra are not corrected with response curves.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We notice that cosmic rays frequently pollute MRS single-exposure spectra and are neither identified nor removed by the LAMOST pipeline. Therefore, in our method, the first step is to carefully remove cosmic rays in spectra. We smooth spectra with a 21-pixel median filter followed by a 9-pixel Gaussian filter, and remove the original pixels that deviate from the smoothed spectrum by four and eight times the local standard deviation in the upper and lower directions. The parameters of the filters are set empirically so that the absorption lines in the spectra of A-, F-, G-, and K-type stars are not affected. The removed pixels are replaced by linearly interpolated values using neighboring pixels. We also note that the two ends of both the blue and red arms are sometimes tilted and show unreasonable flux values probably due to extrapolation of the sky flux modeling, so we trim 50 pixels at both edges of each arm.

The second step is to normalize spectra to a pseudocontinuum to place the spectral features (usually absorption lines) on the same scale. We iteratively fit the spectrum with a smoothing spline function and clip the pixels away from the median values by three times the standard deviation in each 100 Å window. The number of iterations is set to 3. As shown in Figure 2, for a typical K-type star with a medium S/N, the normalization is quite adequate for the subsequent RV measurements.

We adopt the synthetic grid published by Allende Prieto et al. (2018) based on the ATLAS9 stellar atmosphere model (Kurucz 1979; Mészáros et al. 2012) as our spectral library. We degrade the spectral resolution from 10,000 to 7500 to fit the MRS configuration and convert the air wavelength to vacuum wavelength using the formula proposed by Morton (2000). To limit the computational cost of the subsequent RV measurements to a reasonable amount, we interpolate the synthetic library to generate 100 spectral templates with stellar parameters randomly drawn from a uniform distribution in the ranges of 3500 < T_eff/K < 15,000, $0\lt \mathrm{log}g/\mathrm{dex}\lt 5$, −2 < [Fe/H]/dex < 0.5, and −0.5 < [α/Fe]/dex < 0.7. Extremely metal-poor and extremely hot templates are not considered because the spectral features are not significant. Tests on high-S/N MRS spectra show that the sparsity of our templates induces statistical errors ∼0.10 and 0.20 km s⁻¹ in the blue and red arms, respectively, which are negligible as compared to other sources of uncertainties. Figure 3 shows the distribution of parameters of the 100 spectral templates and a series of PARSEC isochrones (Bressan et al. 2012) with solar metallicity and logarithmic age ${\mathrm{log}}_{10}\tau =7$, 8, 9, and 10. These spectral templates are normalized to a pseudocontinuum in the same way as in Section 3.1.

Zoom In Zoom Out Reset image size

The cross-correlation function (CCF; Tonry & Davis 1979) method is widely used in spectroscopic surveys to measure stellar RVs (e.g., Nidever et al. 2015; Steinmetz et al. 2020b). One important advantage is that the CCF can be accelerated using fast Fourier transformation (FFT) once a spectrum is continuum-subtracted and resampled to a logarithmic wavelength grid. However, the drawback of such a scheme is that the sampling of the resulting CCF is generally very sparse. To evaluate the CCF at a smaller RV step, we do not follow the FFT way. In our implementation, the CCF at an RV of v is evaluated as

$$\begin{eqnarray}&&\mathrm{CCF}(v| {\boldsymbol{F}},{\boldsymbol{G}})=\displaystyle \frac{\mathrm{Cov}\left({\boldsymbol{F}},{\boldsymbol{G}}(v)\right)}{\sqrt{\mathrm{Var}({\boldsymbol{F}})\mathrm{Var}({\boldsymbol{G}}(v))}},\end{eqnarray} \tag{ 1 }$$

where F is the vector of the normalized observed spectrum, G (v) is the vector of the normalized synthetic spectrum shifted by an RV (v) and resampled to the wavelength grid of F , Var represents the variance operator, and Cov represents the covariance operator (see Appendix A for more details).

Deriving the final RV consists of three steps:

• 1.  
  The initial estimates are made from an RV grid of −1500 to 1500 km s⁻¹, with a step of 10 km s⁻¹. The template with the maximum CCF value is selected as the best-match template, and the corresponding RV is adopted as the initial guess of the final RV of the observed star. The parameters of the template (T_eff, $\mathrm{log}g$, [Fe/H], and [α/Fe]) are recorded, which make a good prior for some following analyses, such as for stellar atmospheric parameter determination.
• 2.  
  With the best-match template, we maximize the CCF to determine the final RV (v_obs) using the optimization routine scipy.optimize.minimize with the Nelder–Mead algorithm (Nelder & Mead 1965). The corresponding CCF value is recorded as CCFMAX to assess the likelihood between the best-match template and the observed spectrum. The S/N–CCFMAX relations are shown in Figure 4.
• 3.  
  To obtain the measurement error σ_v,obs, we use a Monte Carlo method—namely, we repeat this process 100 times and each time we add Gaussian random noise to the spectrum according to the flux error. The measurement error is computed using the 16th and 84th percentiles, i.e., σ_v,obs = (v₈₄ – v₁₆)/2. The S/N–σ_v,obs relations are shown in Figure 5.

We avoid Gaussian fitting to the CCF, which is widely used in literature (e.g., Nidever et al. 2015; Wang et al. 2019). Our reasons for doing this include the fact that the CCF peak is obviously non-Gaussian and that at low S/N the fitting process is fragile. The final error of the RV includes the measurement error and a noise floor term that is due to systematics (e.g., background, detector imperfections, temperature changes, or focusing issues) and will be assessed below.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this work, one RV is estimated for the blue arm (v_B) and two for the red arm (v_R and v_Rm). The v_Rm is measured using the Hα-masked red-arm spectrum. As shown in Figures 4 and 5, these RV measurements deteriorate rapidly at S/N < 20. For cool stars (e.g., FGK types), at a given S/N, v_B is more precise than v_R because of rich spectral features in the blue arm (for a detailed discussion of spectral information content in MRS spectra, see Zhang et al. 2020b). However, for hot stars, v_R is more reliable than v_B because of the Hα feature and v_Rm is significantly less precise than v_R due to the absence of Hα. Therefore, we recommend that readers only consider v_Rm when the targets have (possible) Hα emissions.

LAMOST has 16 spectrographs of which each has 250 fibers (4000 fibers in total). Excluding a few tens of sky fibers and a few problematic fibers, each spectrograph typically produces ≲200 spectra in an exposure, depending on the targeting, the condition of the instrument, the data quality, and the reduction algorithm. Let i denote the exposure epoch or LMJM, j the spectrograph ID, and k the fiber ID; ideally, we seek the solution of the RVZP for each fiber, for each spectrograph, and exposure by exposure, namely, Δv_i,j,k . This scheme is infeasible because the true/reference RVs v_i,j,k,abs of the targets are not always known.

In this work, assuming that the fibers in a spectrograph in one exposure (hereafter, we refer to it as a spectrograph exposure unit, or an SEU) share similar RVZPs, the systematic RVZP Δv_i,j can be determined as long as a homogeneous reference set of RVs can be found for a fraction of fibers in that SEU. The assumption is quite reasonable given the fact that the wavelength calibration of a multifiber spectrograph is done by fitting a 2D grating equation. And as we will see in Section 4.2, the Gaia DR2 RVs (Katz et al. 2019) meet our needs for the reference set.

As a contrast, both the LAMOST pipeline and Wang et al. (2019) calculate Δv_j assuming the RVZPs for a specific spectrograph do not vary with time, and, therefore, get around the temporal variation of RVZPs. However, we notice that there exist some weird absolute RVs (as we will see in the results in Section 5).

Thanks to the European Space Agency's Gaia mission (Gaia Collaboration et al. 2016), we now have access to the largest RV data set, which matches the LAMOST MRS in terms of velocity precision and magnitude limit. The spectral resolution of Gaia-RVS (R ∼ 11,500; Cropper et al. 2018) is slightly higher than that of the MRS (R ∼ 7500). The magnitude limit of the Gaia DR2 RV catalog (Katz et al. 2019) is at G_RVS = 12 mag or G ∼ 14 mag depending on spectral type and line-of-sight interstellar extinction, which is slightly brighter than the LAMOST MRS magnitude limit. Gaia DR2 contains qualified median RVs for 7,224,631 stars derived from the Gaia-RVS spectra, with T_eff in the range [3550, 6900] K, excluding large RV-variant stars (see Katz et al. 2019 for details). At the faint end, G_RVS = 11.75 mag, the precisions for T_eff = 5000 and 6500 K are 1.4 and 3.7 km s⁻¹, respectively.

Aiming at studying time-domain astrophysical phenomena, e.g., spectroscopic binaries, we proceed to carrying out the second-best scheme—Δv_i,j . Cross-matching the LAMOST MRS DR7 catalog with Gaia DR2, we find 1,582,948 out of the 3,753,659 single-exposure spectra (42.1%) have Gaia RVs, and the common objects usually have good S/N in the MRS because they are relatively bright in that survey. The number of objects in Gaia DR2 is ∼1000 times larger than that in catalogs of RV standard stars, such as that of Huang et al. (2018). The challenge arises because not all of the seven million objects in the Gaia-RVS catalog are RV-invariant, i.e., quite a number of them are pulsating stars or binary/multiple systems that have periodic/nonperiodic RV variations. Below, we demonstrate a robust method that can determine the RVZP self-consistently for each spectrograph in each exposure (Δv_i,j ) by comparing the observed RVs to the Gaia DR2 RVs without identifying RV standard stars.

Assuming that the RV variables vary with random periods at random phases or nonperiodically, and are not the majority of the observed stars, we can regard them as outliers and use a robust method, e.g., least absolute residual (LAR) regression (or least absolute deviation regression; see Press et al. 2007), to estimate Δv_i,j (the common RV bias shared by the objects in an SEU). From a Bayesian perspective, LAR regression originates from an exponential likelihood while least squares (LSQ) regression comes from a Gaussian likelihood. Utilizing the LAR technique, extreme values have less influence on the fit compared to those in LSQ regression. Besides, since we aim at time-domain analysis, as long as our RVZPs are temporally self-consistent, the absolute scales are not very important.

Using the indices proposed in Section 4.1, for each group of pointings, we construct a global cost function f as follows:

$$\begin{eqnarray}\begin{array}{rcl}f\left({\boldsymbol{\Delta }}{\boldsymbol{v}}\right) & = & {{\rm{\Lambda }}}_{1}\displaystyle \sum _{i}\displaystyle \sum _{j}\displaystyle \sum _{k}\displaystyle \frac{| {v}_{i,j,k,\mathrm{obs}}+{\rm{\Delta }}{v}_{i,j}-{v}_{i,j,k,\mathrm{GAIA}}| }{\sqrt{{\sigma }_{\min }^{2}+{\sigma }_{i,j,k,\mathrm{obs}}^{2}+{\sigma }_{i,j,k,\mathrm{GAIA}}^{2}}}\\ & & +{{\rm{\Lambda }}}_{2}\displaystyle \sum _{i}\displaystyle \sum _{j}\displaystyle \sum _{k}\displaystyle \frac{| {v}_{i,j,k,\mathrm{obs}}+{\rm{\Delta }}{v}_{i,j}-\overline{{v}_{\cdot ,\mathrm{obs}}}| }{\sqrt{{\sigma }_{\min }^{2}+{\sigma }_{i,j,k,\mathrm{obs}}^{2}+{\sigma }_{\cdot ,\mathrm{obs}}^{2}}},\end{array}\end{eqnarray} \tag{ 3 }$$

where Δ v is the vector of {Δv_i,j } for all relevant SEUs in the group of pointings; v_i,j,k,obs and σ_i,j,k,obs are the RV and associated measurement error of the kth star in the SEU {i, j}; v_i,j,k,GAIA and σ_i,j,k,GAIA are the Gaia RV and associated uncertainty of the kth star in SEU {i, j}; σ_min is the noise floor of the measured RV, which indicates the stability of the wavelength calibration, i.e., the dispersion of Δv_i,j,k in SEU {i, j}; Δv_i,j is the RVZP correction value of SEU {i, j}, which is a free variable to be solved; $\overline{{v}_{\cdot ,\mathrm{obs}}}$ and σ_.,obs are the median and scatter of the (RVZP-corrected) measured RVs of the star {i, j, k} in other SEUs; and Λ₁ and Λ₂ are the regularization parameters of the two terms. In this scheme, the first term guarantees that the absolute scale of our RVZP-corrected RVs is close to that of Gaia DR2, while the second term makes use of multiple exposures and guarantees that the relative RVZPs are self-consistent. We set Λ₁ = Λ₂ = 1 so that the final correction of each SEU is determined as an average of the two effects. Then the vector Δ v is determined by minimizing the cost function f, i.e., ${\boldsymbol{\Delta }}{\boldsymbol{v}}=\arg \,\min f$. This algorithm can be implemented by minimizing

$$\begin{eqnarray}\begin{array}{rcl}{f}_{i,j}\left({\rm{\Delta }}{v}_{i,j}\right) & = & {{\rm{\Lambda }}}_{1}\displaystyle \sum _{k}\displaystyle \frac{| {v}_{i,j,k,\mathrm{obs}}+{\rm{\Delta }}{v}_{i,j}-{v}_{i,j,k,\mathrm{GAIA}}| }{\sqrt{{\sigma }_{\min }^{2}+{\sigma }_{i,j,k,\mathrm{obs}}^{2}+{\sigma }_{i,j,k,\mathrm{GAIA}}^{2}}}\\ & & +{{\rm{\Lambda }}}_{2}\displaystyle \sum _{k}\displaystyle \frac{| {v}_{i,j,k,\mathrm{obs}}+{\rm{\Delta }}{v}_{i,j}-\overline{{v}_{\cdot ,\mathrm{obs}}}| }{\sqrt{{\sigma }_{\min }^{2}+{\sigma }_{i,j,k,\mathrm{obs}}^{2}+{\sigma }_{\cdot ,\mathrm{obs}}^{2}}}\end{array}\end{eqnarray} \tag{ 4 }$$

for each SEU {i, j} sequentially and iteratively, where Δv_i,j is the RVZP correction value for SEU {i, j}. Therefore, the problem is to solve the vector Δ v that has N_SEU elements, where N_SEU is the number of related SEUs. We claim that Δ v is determined when the L_∞ norm of the difference between the solutions in the lth and (l + 1)th iterations is less than a specified value, i.e., $\max (| {\boldsymbol{\Delta }}{{\boldsymbol{v}}}_{l}-{\boldsymbol{\Delta }}{{\boldsymbol{v}}}_{l+1}| )\lt \epsilon$, where is the tolerance and is set to 0.075 km s⁻¹.

SEUs with RVZP correction values ∣Δv_i,j ∣ > 50 km s⁻¹ or associated uncertainties ${\sigma }_{{\rm{\Delta }}{v}_{i,j}}\gt 10$ km s⁻¹ (see Section 4.5) are excluded from the iterations. These results are generally due to (1) the number of Gaia DR2 objects being ≤10, (2) the spectral S/N being too low, or (3) bad spectra due to saturation or instrumental problems. We do not think that for these SEUs our scheme and assumptions are valid, so we only keep their initial guesses of Δv_i,j (see Section 4.4) and their uncertainties (see Section 4.5) in our catalog (see Section 6).

Several tricks are used to accelerate the algorithm. The first trick is to get a good initial estimation of Δ v . A good approximation can be made by ignoring the second term in Equation (4), so that with Gaia DR2 RVs we can roughly estimate Δv_i,j by minimizing

$$\begin{eqnarray}&&{f}_{i,j,\mathrm{init}}\left({\rm{\Delta }}{v}_{i,j}\right)=\sum _{k}\displaystyle \frac{| {v}_{i,j,k,\mathrm{obs}}+{\rm{\Delta }}{v}_{i,j}-{v}_{i,j,k,\mathrm{GAIA}}| }{\sqrt{{\sigma }_{\min }^{2}+{\sigma }_{i,j,k,\mathrm{obs}}^{2}+{\sigma }_{i,j,k,\mathrm{GAIA}}^{2}}}.\end{eqnarray} \tag{ 5 }$$

We note that if an SEU has only a few objects in common with the Gaia catalog, the estimation is risky. Therefore, we require that an SEU at least have 10 objects in common with Gaia DR2 to proceed; otherwise we calculate the initial guess Δv_i,j,init but exclude this SEU in the iteration process.

The second trick is to cut down N_SEU by separating physically detached SEUs, which hastens the index evaluation in each iteration. In Equations (3) and (4), the second terms contain cross-terms, meaning that when solving the N_SEU elements, an iteration process is needed to guarantee that the solution of Δ v is stable. However, the evaluation of indices is computationally expensive when N_SEU grows. Therefore, before the optimization process, physically detached sky areas can be separated, so that the many index evaluation processes can be accelerated by cutting down the array sizes. We group the pointings of LAMOST MRS DR7 using a friends-of-friends algorithm with a 5° linking length, which is double the radius of the field of view of LAMOST in case of any possible common stars between them. Eventually, 137 groups are obtained, as shown in Figure 6. Initial guesses of Δv for each SEU are made by optimizing the cost function of Equation (3) for each group of plates.

Zoom In Zoom Out Reset image size

In practice, we find that if is too small, there is the possibility that Δ v jumps back and forth between two solutions and does not converge. Further analysis shows that this is an optimization-method-related problem (we use the Nelder–Mead solver; changing it to the Powell solver does not solve the problem but the two solutions are different). We guess that this might be due to the numerical problem of the optimization routine. To avoid such a situation, we then add random processes into the algorithm, i.e., if in the lth iteration the solution is Δ v _l and after looping over all related SEUs the solution is Δ v _l,opt, we evaluate the (l + 1)th solution by Δ v _l+1 = η(Δ v _l,opt − Δ v _l ) + Δ v _l , where η is the learning rate randomly generated between η₀ and η₁. We set η₀ = 0.5, η₁ = 1.0, and = 0.075; considering the expectation for η is 0.75, the effective tolerance of our solution of Δ v is 0.10 km s⁻¹, which is acceptable when compared to the typical precision of measured RVs (e.g., ∼1.5 km s⁻¹ as reported by Wang et al. 2019).

Finally, the RVZP corrections for all the 137 groups of pointings converge after several tens of iterations with a Dell Precision R740 workstation with two Intel Xeon Platinum 8260 CPUs (2.40 GHz), among which the longest solution takes ∼10 hr. Compared to the computation of RV measurements using the CCF for the 3.8 million spectra, including the blue and red arms (∼1 week with the same machine), computing the RVZPs is quite fast.

Rigorous uncertainties are very difficult to obtain for our RVZP corrections. The uncertainties of the RVZP correction values consist of two parts, namely the tolerance in the iteration process and the formal error. The tolerance is = 0.1 km s⁻¹ as mentioned above. For the latter part, based on the discussion presented in Appendix B, we use the 16th and 84th percentiles to construct a fiducial error of our RVZP correction values Δv_i,j divided by an empirical correction ξ to construct the formal error. Hence, the total uncertainties of the RVZP correction values are evaluated via

$$\begin{eqnarray}&&{\sigma }_{{\boldsymbol{\Delta }}{v}_{i,j}}^{2}={\left(\displaystyle \frac{q{84}_{i,j}-q{16}_{i,j}}{2\xi \sqrt{{N}_{i,j}}}\right)}^{2}+{\epsilon }^{2},\end{eqnarray} \tag{ 6 }$$

where i and j index the SEUs, q16 and q84 denote the 16th and 84th percentiles of the residuals of the Gaia DR2 RVs and the RVZP-corrected LAMOST MRS RVs, N_i,j is the number of Gaia DR2 objects with RVs, and ξ is the empirical correction factor for small-number statistics.

In Figure 7, we present the Δv_i,j for each SEU for Sc and ThAr arc lamps versus the date. The Sc lamp was in use until 2018 October, after which it was replaced by the ThAr lamp. The mean uncertainties of Δv_B, Δv_R, and Δv_Rm are all ∼0.38 km s⁻¹, which is quite good. The median uncertainties are even 20% smaller. The Δv_B and Δv_R have different patterns while Δv_R and Δv_Rm are very similar. In Figure 8 we show the distribution of Δv_B, Δv_R, and Δv_Rm of the SEUs solved. The μ here is estimated using the median, and σ is estimated using (q84 − q16)/2. The μ are 0.49 and 6.47 km s⁻¹ for the ThAr and Sc lamps in the blue arm, while in the red arm they are 0.26 and 4.91 km s⁻¹, respectively. The σ are 1.07 and 1.06 km s⁻¹ for the ThAr and Sc lamps in the blue arm, and in the red arm they are 0.85 and 0.68 km s⁻¹. Generally, the different systematics for the Sc- and ThAr-lamp-calibrated data are consistent with Wang et al. (2019). Despite the large systematics, we find that the precision of the Sc-lamp-calibrated data is no worse than that of the data calibrated by the ThAr lamp. The Rm results are very similar to those of R, except their μ have a 0.3 km s⁻¹ difference. We note that this is reasonable considering that the systematics vary with wavelength as shown in Ren et al. (2021). In addition, ∼4000 single-exposure spectra calibrated using an Ne lamp are also found in DR7 v1.1. We confirm that the Ne lamp was used to calibrate the LRS spectra and these mistakenly calibrated data will be removed from the internationally available version of DR7, so we exclude these data in the following analysis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is clear that the RVZPs are reasonably stable except after around 2019 May 1, which seems to be correlated with the arc lamp exposure flux. ¹⁴ We plot the Δv for each spectrograph for the time interval with available arc lamp intensity in Figures 9 and 10. Since our mean uncertainty of Δv_i,j is 0.38 km s⁻¹, we regard Δv_i,j larger than 1 km s⁻¹ as significant values, including the 7th, 12th, 13th, and 15th spectrographs of the blue arm and the 7th, 9th, and 15th spectrographs of the red arm. It is currently not known whether these shifts are due to the fact that the new ThAr lamps were brought into use at around 2019 May 1 or to some other issues induced in the maintenance. Probably in DR8, with the RVZP data in a longer time baseline, we can address the problem.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also validate our RVs with stars common to other data sets, namely, Gaia DR2 (Katz et al. 2019), APOGEE DR16 (Column VHELIO_AVG; Jönsson et al. 2020), the RV standard stars from Huang et al. (2018), GALAH DR3 (Column rv_galah; Buder et al. 2021), and RAVE DR6 (Column hrv_sparv; Steinmetz et al. 2020b). The mean μ and scatter σ derived from Gaussian fitting to the residuals are tabulated in Table 2 and shown as functions of S/N in Figure 11. Also shown are the results of the LAMOST pipeline and Wang et al. (2019). Note that, in DR7 v1.1, the LAMOST pipeline only provides two RVs measured from the blue arm using ELODIE empirical templates (Moultaka et al. 2004) and ATLAS9 synthetic templates (Castelli & Kurucz 2003), respectively. In future MRS data releases, the RVs using ELODIE templates will be removed, and the RVs based on ATLAS9 will be provided for both blue and red arms for better performance. Therefore, we use calibrated RVs of blue arms based on ATLAS9 templates in our comparison (Column rv_ku1). The Wang et al. (2019) catalog is a subset of ours (including data taken during the first year and a half with an S/N cut at 10). Spectra without our measurements (i.e., either S/N < 5 or Δv is invalid) are excluded from this comparison, so that it is fair to the other two RV sources. At the high-S/N end (50–100), we find the standard deviations derived from Gaussian fitting for our results can reach 1.00/1.10, 0.84/0.80, 0.69/0.74, 0.72/0.78, and 1.77/1.88 km s⁻¹ with respect to Gaia DR2, APOGEE DR16, Huang et al. (2018), GALAH DR3, and RAVE DR6 in the blue/red arm. The common stars among LAMOST MRS DR7, Gaia DR2, and the other four reference sets are used to calculate the fiducial accuracy of Gaia DR2. At the high-S/N end (50 < S/N < 100), the LAMOST MRS gets close to the performance of Gaia DR2 (see Figure 11). Note that, in our algorithm, the Gaia DR2 RVs are used as a reference set, so the comparison with Gaia DR2 is not an independent validation. These results of the comparisons are quite fascinating. At the high-S/N end, we outperform Wang et al. (2019) and the LAMOST pipeline by ∼20% and 30%, respectively, according to the blue-arm results compared to APOGEE DR16. At the low-S/N end (10–20), the advantages increase to 58% and 47%, indicating that our algorithm of RV measurements and RVZP determinations is quite efficient. Since the Huang et al. (2018) sample is from APOGEE DR14, the comparison to Huang et al. (2018) has a 0.4 km s⁻¹ systematic bias, which does not exist in our comparison to APOGEE DR16. This may be due to the update of the APOGEE data release. GALAH DR3 has a 0.23 km s⁻¹ systematic bias as compared to Gaia DR2. The comparison with RAVE DR6, whose spectral resolution is the same as that of the LAMOST MRS but is lower than that of Gaia-RVS, APOGEE, and GALAH, shows a large scatter but it is still reasonable.

Zoom In Zoom Out Reset image size

Besides the precision, a check on self-consistency is necessary and important before using RVs in time-domain research. As a demonstration, in Figures 12 and 13 we validate the temporal RV variations of the standard stars (whose RVs are assumed invariant) from Huang et al. (2018) in two pointings, namely, planid = HIP8426401 and HIP4312101, which have 10 and 17 exposures and contain 27 and 41 RV standard stars, respectively. The RVs from the LAMOST pipeline show large fluctuations, which are as expected from the comparison in Section 5.2. It turns out that for many stars, the RVs of multiple exposures from the Wang et al. (2019) catalog have exactly the same values. This is due to the failure of their Gaussian fitting process and then they fall back on the best estimation from the 1 km s⁻¹ RV grid, which is a defect in the algorithm. By comparing the measured RVs and the RVZP-corrected RVs, we find that the RVZP-corrected RVs are much cleaner, indicating that the RVs do benefit from our algorithm for Δv_i,j determination. The LAMOST pipeline and Wang et al. (2019) assume that the RVZP for a spectrograph is static so that the RVZP corrections are basically a shift of the RVs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To quantify the performance in self-consistency, we select RV standard stars from Huang et al. (2018) with at least eight exposures having valid RVZP-corrected RVs in our catalog and calculate their standard deviation, empirically corrected for the small-number statistic effect, which is discussed in Appendix B. The cumulative distribution function (CDF) of the standard deviations is shown in Figure 14. The absolute RVs from the blue arm (B) show the best self-consistency, followed by those from the red arm (R) and the Hα-masked red arm (Rm). By calculating the 50th, 90th, and 95th percentiles of the CDFs for the blue-arm results, we find our absolute RVs have significant advantages over those of the LAMOST pipeline, namely 16.5%, 35.5%, and 37.5% better (see Table 3), while Wang et al. (2019) are at nearly the same level as the LAMOST pipeline. This reveals the excellent self-consistency of our absolute RVs, which will be used in the time-domain analysis. Note that these estimations of advantages are very conservative due to the fact that the Wang et al. (2019) data set cuts the S/N at 10 but LAMOST and our sample cut at 5, for example. Besides, we do not exclude any of the "exactly the same" RVs from Wang et al. (2019) as readers may find that the CDFs for LAMOST and Wang et al. (2019) jump at around the 0 km s⁻¹ position.

Zoom In Zoom Out Reset image size

With this excellent self-consistency, we select 10,320 candidate RV standard stars by requiring that in the blue arm (B) and red arm (R)

• 1.  
  their numbers of exposures are at least 8,
• 2.  
  their absolute RVs have standard deviations (corrected for small-number statistics) of less than 1.45 and 1.85 km s⁻¹ (corresponding to the 95% level of the CDF, or to 95% completeness with respect to Huang et al. 2018), and
• 3.  
  their time baselines are longer than at least 180 days.

These stars can be useful for the RV calibration of low-resolution surveys, such as the LAMOST LRS (R ∼ 1800).

We can minimize an L₁-norm or L₂-norm cost function, namely, ${\sum }_{i}| {x}_{i}-\hat{\mu }|$ or ${\sum }_{i}{\left({x}_{i}-\hat{\mu }\right)}^{2}/2$, respectively, to get an estimate of the mean $\hat{\mu }$. The true deviation is by definition

$$\begin{eqnarray}&&{\delta }_{\mathrm{true}}=| \hat{\mu }-\mu | .\end{eqnarray} \tag{ B1 }$$

However, in practice we do not know μ when we tackle such a problem. A fiducial deviation associated with $\hat{\mu }$ can be constructed using the 84th and 16th percentiles (or the interquantiles; see Lupton 1993; Ivezić et al. 2014), i.e.,

$$\begin{eqnarray}&&{\delta }_{\mathrm{est}}=\displaystyle \frac{{\{{x}_{i}\}}_{q84}-{\{{x}_{i}\}}_{q16}}{2\sqrt{N}},\end{eqnarray} \tag{ B2 }$$

where N is the sample size. To obtain an empirical relation between δ_est and δ_true, we assume the following form:

$$\begin{eqnarray}&&{\delta }_{\mathrm{true}}={\delta }_{\mathrm{est}}/\xi (N),\end{eqnarray} \tag{ B3 }$$

where ξ(N) is the empirical correction factor and is a function of N. Then, we draw mock data from a standard Gaussian distribution with the numpy.random module. In each experiment, we draw N samples and calculate δ_est and δ_true, and derive ξ. We repeat this experiment 3000 times for each N, which ranges from 2 to 250 (which is enough for our purposes), and show the 16th, 50th, and 84th percentiles of the results for each N in the left panel of Figure B1. This fiducial deviation is an underestimation of δ_true. The relation between the medians of logξ and logN is fitted with a fifth-order polynomial function, whose coefficients are tabulated in Table B1. With this relation, we can scale the fiducial deviation to a standard that is less affected by the sample size N.

Zoom In Zoom Out Reset image size

Table B1. The Best-fit Coefficients of the Fifth-order Polynomials for the Empirical Relationship between log₁₀ ξ and log₁₀ N

Cost Function	β5	β4	β3	β2	β1	β0
${\sum }_{i}| {x}_{i}-\hat{\mu }|$	0.07349721	−0.60647022	1.97806105	−3.22994084	2.72007585	−0.92812989
${\sum }_{i}{\left({x}_{i}-\hat{\mu }\right)}^{2}/2$	0.08434813	−0.69694429	2.28047526	−3.73821453	3.15612997	−0.99158461

Note. Our definition of a polynomial is poly(x∣β) = ∑_i β_i xⁱ .

Download table as:  ASCIITypeset image

For a Gaussian distribution P(x∣μ, σ), we can use a sample-based standard deviation and the 16th and 84th percentiles to estimate the true standard error σ, i.e.,

$$\begin{eqnarray}&&{\sigma }_{\mathrm{est}}=\sqrt{\displaystyle \frac{\sum _{i}{\left({x}_{i}-\hat{\mu }\right)}^{2}}{N-1}},\end{eqnarray} \tag{ B4 }$$

and

$$\begin{eqnarray}&&{\sigma }_{\mathrm{est}}=\displaystyle \frac{{\{{x}_{i}\}}_{q84}-{\{{x}_{i}\}}_{q16}}{2},\end{eqnarray} \tag{ B5 }$$

respectively. We define ζ by

$$\begin{eqnarray}&&{\sigma }_{\mathrm{true}}={\sigma }_{\mathrm{est}}/\zeta (N).\end{eqnarray} \tag{ B6 }$$

The results clearly show that both methods underestimate the standard error. Similar to the procedures in the previous test, we fit the relation with a fifth-order polynomial to the medians of ζ and ${\mathrm{log}}_{10}N$; the best-fit polynomials are shown in the right panel of Figure B1 and the coefficients are tabulated in Table B2. Compared to Equation (B4), Equation (B5) is more robust to outliers but suffers from more significant underestimation when N is small.

Table B2. The Best-fit Coefficients of the Fifth-order Polynomials for the Empirical Relationship Between ζ and log₁₀ N

Estimator	β5	β4	β3	β2	β1	β0
(q84−q16)/2	0.05712779	−0.47050592	1.56888875	−2.75649024	2.71900493	−0.28637106
$\sqrt{\tfrac{1}{N-1}\sum {\left({x}_{i}-\hat{\mu }\right)}^{2}}$	0.08344418	−0.68630504	2.21084202	−3.50529186	2.77854093	0.08576048

Download table as:  ASCIITypeset image