On the Use of Field RR Lyrae as Galactic Probes. V. Optical and Radial Velocity Curve Templates

We collected the largest spectroscopic catalog of RR Lyrae (RRLs) including ≈20,000 high-, medium-, and low-resolution spectra for ≈10,000 RRLs. We provide the analytical forms of radial velocity curve (RVC) templates. These were built using 36 RRLs (31 fundamental—split into three period bins—and five first-overtone pulsators) with well-sampled RVCs based on three groups of metallic lines (Fe, Mg, Na) and four Balmer lines (H α , H β , H γ , H δ ). We tackled the long-standing problem of the reference epoch to anchor light-curve and RVC templates. For the V-band, we found that the residuals of the templates anchored to the phase of the mean magnitude along the rising branch are ∼35% to ∼45% smaller than those anchored to the phase of maximum light. For the RVC, we used two independent reference epochs for metallic and Balmer lines and we verified that the residuals of the RVC templates anchored to the phase of mean RV are from 30% (metallic lines) up to 45% (Balmer lines) smaller than those anchored to the phase of minimum RV. We validated our RVC templates by using both the single-point and the three phase point approaches. We found that barycentric velocities based on our RVC templates are two to three times more accurate than those available in the literature. We applied the current RVC templates to Balmer lines RVs of RRLs in the globular NGC 3201 collected with MUSE at VLT. We found the cluster barycentric RV of V γ = 496.89 ± 8.37(error) ± 3.43 (standard deviation) km s−1, which agrees well with literature estimates.

Pulsating variables are behind numerous breakthroughs in astrophysics. Classical Cepheids (CCs) were used to estimate the distance to M31 and solve the Great Debate concerning the extragalactic nature of the so-called Nebulae (Hubble 1926) and to trace, for the first time, the rotation of the Galactic thin disc (Oort 1927;Joy 1939). The size and the age of the Universe were revolutionized thanks to the discovery of the difference between CCs, RRLs, and type II Cepheids (TIICs, Baade 1956). Indeed, while they were previously thought to represent the same type of variable stars, it became clear that they represented very distinct populations, with the RRLs and TIICs being very old (t≥10 Gyr), and the CCs very young (t≤300 Myr).
Nowadays, CCs are among the most popular calibrators of the extragalactic distance scale (Riess et al. 2019). RRLs, albeit fainter, are excellent standard candles that can provide robust, independent distance measurements even for stellar populations where the young CCs are absent. RRLs obey the well-defined Period-Luminosity-Metallicity (PLZ) relations for wavelengths longer than the R-band (Bono et al. 2003;Catelan 2009). As tracers of purely old stellar populations, they can be used to investigate the early formation and evolution of both the Galactic Halo (Fiorentino et al. 2015;Belokurov et al. 2018;Fabrizio et al. 2019) and Bulge (Pietrukowicz et al. 2015;Braga et al. 2018).
It is noteworthy that we lack general consensus on the Galactic Halo structure, in part because different stellar tracers provide different views concerning its spatial structure and the timescale for its formation. Indeed, Carollo et al. (2007) by using Main Sequence, subgiants and RGs and Kinman et al. (2012) by using RRLs, suggested that the outer halo is more spherical and its density profile is shallower when compared with the inner halo. In contrast, Keller et al. (2008) by using RRLs and Sesar et al. (2010Sesar et al. ( , 2011 by using RRLs plus Main Sequence stars suggested that the outer halo has a steeper density profile when compared with the inner halo. Deason et al. (2011) by using Blue Horizontal Branch stars and Blue Stragglers found no change in the flattening as a function of the Galactocentric distance (Sesar et al. 2011). More recently, Xue et al. (2015), by adopting a global ellipsoidal stellar density model with Einasto profile found that the models with constant flattening provide a good fit to the entire Halo.
The tension between different measurements may be due to the sample selection of each study. On the one hand, the ages of the RRLs cover a narrow range from ∼10 to ∼13 Gyrs. There is evidence that a few RRLs-or stars that mimic RRLs, see Smolec et al. (2013)-are the aftermath of binary evolution, but they only represent a few percent of the populations (Bono et al. 1997;Pietrzyński et al. 2012;Kervella et al. 2019). On the other hand, Red giants (RGs) and main sequence (MS) stars, typically used to investigate the Halo, have only very weak age constraints (Conroy et al. 2021). Indeed, all stellar structures less massive than 2M (older than ∼0.5-1.0 Gyr) experience a RG phase and MS stars also cover a broad range in stellar masses/ages. This means that if the Halo is the result of an in-tense disruption and merging activity (Monachesi et al. 2019) RG and MS stars are far from being optimal tracers of the early formation, because they are a mixed bag concerning the age dis-tribution.
Field RRLs are less numerous when compared with RG and MS stars, but their narrow age distribution makes them uniquely suited for Galactic archaeology. They probe a significant Halo fraction (Galactocentric distance ≤150 kpc) with high accuracy. Their individual distances have uncertainties on average smaller than 3-5% and their accuracy improves when moving from optical to NIR (Longmore et al. 1986;Catelan et al. 2004;Braga et al. 2015). This is a key advantage even in the Gaia era: Gaia EDR3 has an accuracy of 3% for Halo RRLs (G≤15 mag) at 1 kpc and this accuracy will be extended to 2 kpc at the end of the mission (Clementini et al. 2019). RRLs are also valuable targets from the kinematical point of view. In fact, by measuring their velocities, one gets information on the kinematical state of the old population (Halo, Globular Clusters, Bulge). The pioneering work by Layden (1994Layden ( , 1995, based on 302 RRLs, pointed towards a non-steady formation of the Halo, favouring a fragmented accretion scenario (Searle & Zinn 1978). More recently, Zinn et al. (2020) were able to pinpoint the membership of several Halo RRLs to past merger events (Gaia-Enceladus and the Helmi streams, Helmi et al. 1999;Myeong et al. 2018;Helmi et al. 2018). A few Halo RRLs were also associated with the Orphan stream by Prudil et al. (2021), leading to more solid constraints on the origin of the stream itself. Concerning the Bulge, the kinematic properties of RRLs display a duality, with one group of stars associated with the spheroidal component and the other with the Galactic bar (Kunder et al. 2020).
The number of identified RRL is rapidly growing thanks to the enhancements in telescope collecting areas and instrument efficiency. Thanks to long-term optical (Catalina, Drake et al. 2009Drake et al. , 2017ASAS, Pojmanski 1997; ASAS-SN, Jayasinghe et al. 2019;DES, Stringer et al. 2019;Gaia, Clementini et al. 2019;OGLE, Soszyński et al. 2019;Pan-STARRS Sesar et al. 2017) near-infrared (VVV, VVV-X, Minniti et al. 2011) and mid-infrared (neo-WISE, Wright et al. 2010) surveys, more than 200,000 RRLs were identified in the Galactic spheroid. However, RRLs are demanding targets from an observational point of view. Wellsampled time series, meaning at least a dozen, properly sampled, photometric measurements, are required for a solid identification and an accurate characterization. The same limitation applies to the measurement of the RRL barycentric radial velocity (V γ ), because it requires multiple measurements to trace the radial velocity (RV) variation along the pulsation cycle. To overcome this limitation, several authors have used the radial velocity curve (RVC) of X Ari, observed more than half a century ago by Oke (1966), as a pseudo-template. More recently, RVC templates have been developed for fundamental (RRab) RRLs (Liu 1991;Sesar 2012, henceforth, S12). They allow to estimate V γ even with a small number of velocity measurements, provided that the V-band pulsation properties are known. The current RVC templates are affected by several limitations: despite being based on 22 RRab stars with periods between 0.37 and 0.71 days, the template of Liu (1991) was derived from RVCs with-at most-a few tens of points each. These points are velocities obtained from a heterogeneous set of unidentified metallic lines, since they were collected from several different papers. S12 provided templates for both metallic and H α , H β and H γ lines with a few hundreds of RV measurements. However, their Balmer templates do not cover the steep decreasing branch and, even more importantly, the templates were based on only six RRab with periods in a very narrow range (0.56-0.59 days). Finally, no RVC templates are available for first-overtone RRLs (RRc).
This work aims at providing new RVC templates for both RRab and RRc variables by addressing all the limitations described above. We adopted a wide set of specific and well-identified metallic and Balmer lines for both RRab and RRc stars and hundreds of velocity measurement for each template. As the velocity curves of the RRab display some peculiar variations among themselves, we also separated them into three bins according to their specific shape and pulsation period. Thus, we can provide uniquely precise templates that cover a wide range of intrinsic parameters of these variable stars.
The paper is structured as follows. In section 2 we investigate the phasing of the optical light curve and discuss on a quantitative basis the difference between the reference epoch anchored to the luminosity maximum and to the mean magnitude along the rising branch. We present the spectroscopic dataset in Section 3 and provide new RVCs and their properties in Section 4. We put together the RVCs and derive the analytical form of the RVC templates in Section 5, discuss the reference epoch to be used to apply the templates in Section 6 and validate them in Section 7. We provide a practical example of how to use the RVC template on spectroscopic observations of NGC 3201 in Section 8. Finally, in Section 9, we summarize the current results and outline future developments of this project.

OPTICAL LIGHT CURVE TEMPLATES
Light curve templates are powerful tools that model the light curve of a periodic variable star. The templates are parametrized with the properties of the variable stars (pulsation mode, period and amplitude). These come in hand, e.g., to estimate the pulsation properties with a few available data (Stringer et al. 2019), to obtain O − C diagrams that trace the rate of period change (Hajdu et al. 2021), to predict the luminosity of the star at a given phase, and for various other purposes.
The use of both luminosity and RVC templates relies on the use of a reference epoch. This means that the phase zero of the RVC template has to be anchored to a specific feature of the luminosity/RV curve. The most common reference epoch adopted in the field of pulsating variable stars is the time of maximum light in the optical (T opt max ). The RVC templates available in the literature are also anchored to T opt max because it matches, within the uncertainties, the time of minimum in the RVC T RV min . Note that, by "minimum" in the RVC, we mean the numerical minimum, i.e., the epoch of maximum blueshift. This is an approximate choice due to the well-known phase-lag between light and RVC (Castor 1971).
Our group introduced a new reference epoch, namely, the epoch at which the magnitude along the rising branch of the Vband light curve-that is, the section of the light curve where brightness changes from minimum to maximum-becomes equal to the mean V-band magnitude (T opt mean , Inno et al. 2015;Braga et al. 2019). We thoroughly discussed the advantages of using T opt mean versus T opt max in the context of NIR light curve templates for both CCs and RRLs. The reader interested in a detailed discussion is referred to the quoted papers. Here, we summarize the key advantages in adopting T opt mean for RRL variables. i) RRab variables with large amplitudes have RVCs with a "sawtooth" shape, where the maximum can be misidentified by an automatic analytical fit if the phase coverage is not optimal. The rising branch, however, can be more easily fitted. ii) A significant fraction of RRc variables displays a well-defined bump/dip before the maximum in luminosity. A clear separation between the two maxima is not trivial if the phase coverage is not optimal. iii) The estimate of T opt max is more prone to possible systematics, even with well-sampled light curves, because several RRc and long-period RRab variables display flat-topped light curves i.e. light curves in which the maximum is almost flat for a relatively broad fraction of the phase cycle (∼0.10). iv) T opt max is typically estimated either as the top value of the fit of the light cure or the brightest observed point, when the sampling is optimal (e.g., ASAS-SN). This means that T opt max is affected by the intrinsic dispersion of the observations and by the time resolution of photometric data. Meanwhile, T opt mean is estimated by interpolating the analytical fit the mean magnitude (see Appendix C.1), which is a very robust property of the star. Therefore, T opt mean is intrinsically more robust because its precision is less dependent of sampling.
In the following, we address on a more quantitative basis these key issues in the context of optical light curves. For this purpose, we take advantage of a homogeneous and complete sample of V-band light curves for cluster and field RRL variables. In particular, we use visual light curves for RRLs in M4 (Stetson et al. 2014) and in ω Cen (Braga et al. 2016) together with literature observations for RRLs with Baade-Wesselink (BW) analysis, (Braga et al. 2019, and references therein). The RRLs in M4 and in ω Cen have well-sampled light curves, with the number of phase points ranging from hundreds to more than one thousand.

Phasing of optical light curves
We selected 57 RRLs (7 RRc, 50 RRab) from our M4, ω Cen and BW RRLs and separated them into four period bins (See Fig. 1). Following Braga et al. (2019) the thresholds are the following: "RRc", "RRab1" (RRab with periods shorter than 0.55 days), "RRab2" (RRab with periods between 0.55 and 0.70 days) and "RRab3" (RRab with periods longer than 0.70 days). See Section 4.4 for a more detailed discussion. We normalized all the light curves by subtracting the mean magnitude and dividing them by their peak-to-peak amplitude A(V), and estimated T opt mean and T opt max for the entire sample and the individual values are listed in columns 5 and 6 of Table 1. T opt mean was estimated as described in Appendix C.1 and T opt max was estimated by converting the phase of maximum light of the model light curve, into an Heliocentric Julian Date.   We visually inspected all the reference epochs derived in this work (see Figure 2). To overcome thorny problems in the phasing of light curves, we manually selected the best value of T opt max as the HJD of the phase point closest to the maximum for the variables where the fit does not follow closely the data around maximum light. In contrast, no manual selection of T opt mean was needed because its estimation is, by its own nature, based on a more robust approach. We anchored the phases to both T opt max and T opt mean and we piled up the light curves into four period bins. We ended up with eight cumulative and normalized light curves: four with τ 0 anchored to T opt max and four anchored to T opt mean . Finally, we adopted the PEGASUS (PEriodic GAuSsian Uniform and Smooth) function (a series of multiple periodic Gaussians, Inno et al. 2015) to fit the cumulative and normalized light curves. The form of the PEGASUS fit is: where A 0 and A i are the zero points and the amplitudes of the Gaussians, while φ i and σ i are the centers and the σ of the Gaussians.  The coefficients of the PEGASUS fits are listed in Table 2. The standard deviations plotted to the right of the light curves (see also the last column in Table 2) bring forward two interesting results. i)-The standard deviations of the light curves phased by using T opt mean are systematically smaller than those phased using T opt max . The difference for the period bins in which the light curves display a cuspy maximum (RRab1, RRab2) is ∼37% smaller, but it becomes ∼45% smaller for the RRc and the RRab3 period bins, because they are characterized by flat-topped light curves. ii)-The cumulative light curves for the RRc and RRab3 period bins phased using T opt max show offsets along the rising branch. This mismatch could lead to systematic offsets of ∼30% in Amp(V) adopted to estimate the mean < V > magnitude. Meanwhile, the cumulative light curves phased using T opt mean overlap better with each other over the entire pulsation cycle. There is one exception: the RRab3 period bin shows a marginal difference across the phases of maximum in luminosity, but the error in the adopted Amp(V) is on average a factor of two smaller (∼15%) than those obtained by using T opt max as anchor. The current circumstantial evidence, based on the same photometric data, indicates that the use of a reference epoch anchored to the phase of mean magnitude along the rising branch allows a more accurate phasing with respect to the phase of the maximum in luminosity.

Phase offset between T opt max and T opt mean
We are aware that large photometric surveys-but also smaller projects focused on variable stars-provide, as reference epoch, T opt max . To overcome this difficulty and to provide a homogeneous empirical framework, we investigated the phase offset between T opt mean and T opt max . In particular, we defined the phase difference where mod is the remainder operator. For this purpose, we could adopt a larger sample of visual light curves of 291 RRLs (54 RRc and 237 RRab) from large photometric surveys (Gaia, ASAS, ASAS-SN and Catalina), from our own photometry of globular clusters (ω Cen, M4), and from the literature (BW sample, see caption of Table 11). We found that the phase difference shows, as expected, a trend with the pulsation period (see Figure 2.2). In particular, the RRab variables show a quite clear linear trend of phase offset with period (∆Φ = 0.043 + 0.099 · P), with an intrinsic dispersion of 0.024. The standard deviation for RRc variables is larger, but there is no clear sign of a period dependency. Therefore, we assume a constant phase difference (∆Φ =0.223±0.036) for RRc variables. We also investigated a possible correlation of phase offset with metallicity by adopting the estimates recently provided by Crestani et al. (2021a), but we found none.
To provide new RVC templates we performed a large spectroscopic campaign aimed at providing RV measurements for both field and cluster RRLs. We reduced and analyzed a large sample of high-, medium-and low-resolution (HR, MR, LR) spectra. This mix of proprietary data and data retrieved from public science archives was supplemented with RVCs of RRLs available in the literature.

Spectroscopic catalog
We collected the largest spectroscopic dataset-both proprietary and public-for RRLs. Preliminary versions of this spectroscopic catalog were already used in studies focused on chemical abundances (Fabrizio et al. 2019;Crestani et al. 2021a,b) and RV (Bono et al. 2020a). In this investigation, we added new spectroscopic data and discuss in detail the spectra used for RV measurements. We ended up with 23,865 spectra for 10,413 RRLs. Figure 3.1 shows that the distribution of the RRLs is well-spread over the Galactic Halo. The key properties of the spectra (spectral resolution and signal-to-noise ratio), the spectrographs and the spectroscopic sample are summarized in Table 3. . Aitoff projection in Galactic coordinates (l, b) of the RRL spectroscopic dataset. Blue and red symbols display RRc and RRab variables. High-, medium-, and low-resolution spectra are marked with large circles, crosses, and small circles, respectively. RRL for which we only have radial velocity measurements from the literature are plotted as diamonds (Baade-Wesselink sample, see Section 4.2). Black stars indicate the RRLs used to build the RVC templates.
The HR sample mainly includes spectra collected with the Las Campanas Observatory du Pont echelle spectrograph (du Pont, 6,208 spectra), plus HR spectra collected from ESO telescopes (277 from UVES@VLT, 320 from HARPS@3.6m, 55 from FEROS@2.2m MPG). We also have 100 HR spectra from SES@STELLA, 81 from HRS@SALT, 10 from HARPS-N@TNG and 34 from HDS@Subaru. We collected MR spectra from both X-Shooter@VLT (121 spectra) and the LAMOST MR survey (1271 spectra). Finally, our spectroscopic dataset includes LR spectra from the LAMOST (9,099 spectra) and from the SDSS-SEGUE (6,289 spectra) surveys. Note-Each row gives either the spectrograph or the spectroscopic dataset (column 1), the total number of spectra (column 2), the number of RRab and RRc variables (column 3 and 4), the typical spectral resolution (column 5) and the typical SNR@3950 • A(column 6).
The main aim of this investigation is to provide RVC templates that can be used to provide V γ for RRLs from a few random RV measurements based on a wide variety of spectra. For this purpose, we selected a broad range of strong and weak spectroscopic diagnostics.

Radial velocity spectroscopic diagnostics
The decision to use multiple spectroscopic diagnostics was made because different lines form at different atmospheric layers. As the RRL are pulsating stars, different lines may trace very different kinematics even when observed at the same phase. The resulting velocity curves for different lines, therefore, may have different shapes and amplitudes. Consequently, combining different hydrogen and/or metallic lines for a single velocity determination would blur the fine detail of the velocity curves and decrease the accuracy of the V γ estimate. With this in mind, we performed RV measurements separately with the following diagnostics: four Balmer lines (H α , H β , H γ and H δ ), the Na doublet (D1 and D2), the Mg I b triplet (Mg b 1 , Mg b 2 and Mg b 3 ) and a set of Fe and Sr lines (three lines of the Fe I multiplet 43 and a resonant Sr II line, Moore 1972).
The laboratory wavelengths of the quoted absorption lines are listed in Table 4. Figure 4 displays the regions of the spectrum of four RRLs where the quoted lines are located. The four RRLs were selected in order to have one RRL for each period bin of the RVC template (see Section 4.4). To measure the RVs for the quoted diagnostics we performed a Lorentian fit to the absorption lines by using an automated procedure written in IDL. The wavelength range adopted by the fitting algorithm is fixed according to the spectral resolution of the different spectrographs. Typically, we selected a range in wavelength that is ten Full Width at Half Maximum (FWHM) to the left and ten to the right. The FWHM was estimated as FWHM=2.355 × λ obs R , where λ is the wavelength of the diagnostic and R is the spectral resolution.
The median uncertainties of the single RV estimates for the adopted spectroscopic diagnostics and the standard deviations of the different datasets are listed in Table 5. Note that the different datasets have median uncertainties, on average, smaller than 1.5 km/s.  Table 5. Typical uncertainties in radial velocity measurements for the adopted diagnostics in the different spectroscopic datasets. To complement our dataset, we collected RVCs of RRLs from the literature (Liu & Janes 1989;Skillen et al. 1993b,a;Jones et al. 1988a,b;Cacciari et al. 1987;Jones et al. 1987a;Clementini et al. 1990;Fernley et al. 1990;Jones et al. 1987b). During the 80s and 90s of the previous century, several bright RRLs were observed both photometrically (optical and NIR) and spectroscopically (velocities from metallic lines) to apply the Baade-Wesselink method (BW, Baade 1926;Wesselink 1946) in order to obtain accurate distance determinations. Therefore, we label the set of RVCs from these works as the BW sample. Unfortunately, Figure 5. From top to bottom, HR spectra collected with du Pont for RRLs adopted in the four bins of the RVC template. The flux units are arbitrary. The Gaia ID (alternative name in parentheses), iron abundance ([Fe/H]) and the pulsation period are labeled. The four portions of the spectrum display, from left to right, the metallic lines, Mg I b triplet, Na doublet and the Balmer lines. Each line is marked in blue and labeled. All the spectra presented in this figure were taken at phases near one third of the rising branch of the RV curve and are only minimally affected by nonlinear phenomena. The first row shows the spectrum of the RRc variable CS Eri. Second row: same as the top, but for the RRab HH Pup. Third row: same as the top but, for the RRab RV Oct. Fourth row: same as the top but, for the RRab AT Ser.

Radial velocity curves from the literature
it was not possible to collect the spectra, therefore we adopted the RV estimates as provided in the quoted papers. Overall, the BW sample includes 2,725 RV measurements for 36 RRab and 3 RRc.
Although this dataset is inhomogeneous and based on a mix of weak metallic lines, it is extremely useful to complement our own measurements. Some of the works mentioned above included optical light curves, which we used to validate the robustness of the reference epoch used in the phasing of the RVC template.

Estimate of barycentric velocities, radial velocity amplitudes and reference epochs
To derive the analytic form of the RVC templates, it is necessary to know the pulsation period (P), the reference epoch, V γ and the RV amplitude (Amp(RV)) of the RRLs with a well-sampled RV curve. The former two are needed to convert epochs into phases and the latter two are used for the normalization of the RV curve. The normalization is a crucial step because the RVC templates have to be provided as normalized curves, with zero mean and unit amplitude.
Our data set includes RV measurements for more than 10,000 RRLs, but only 74 of them have a well-sampled pulsation cycle. A good phase coverage is necessary for the determination of the pulsation properties required for the creation of the RVC template. Reference epochs and Amp(RV) are particularly sensitive to the quality of the pulsation cycle sampling. We neglected all the RRLs displaying a clear Blazhko effect (a modulation of the pulsation amplitude, both in light and in RV) that would introduce a large intrinsic spread in the RVC templates. Because of this exacting quality control, we derived the analytic form of the RVC template using only a subset of three dozen RRL in the spectroscopic template (31 RRab, 5 RRc). They cover a broad range in pulsation periods (0.27-0.84 days) and iron abundances (-2.6 ≤ [Fe/H]≤-0.2). We label these stars with the name a The iron abundances are all taken from homogeneous metallicity estimates in Crestani et al. (2021a). The only exception is Cl* NGC 6341 SAW V1, for which we have adopted the abundance published in Kraft & Ivans (2003) "Template Sample" (TS) and their properties are listed in Table 6. The individual RV measurements for the TS variables are given in Table 7.
The RRLs in the TS have well-covered RVCs for all the adopted spectroscopic diagnostics. The only exception is a cluster star (Cl* NGC 6341 SAW V1) that has good RVCs only for Fe and Mg lines. The number of calibrating RRLs adopted in this work is six times larger than the RRL sample adopted by S12. Moreover, S12 only included RRab variables covering a limited range in pulsation periods (0.56-0.59 days). To estimate V γ , Amp(RV), the epoch of mean velocity on the decreasing branch and the epoch of minimum velocity (T RV mean and T RV min , both for Fe and H β RVCs), we fitted the RVCs with the PLOESS algorithm, as described in Bono et al. (2020a). Then, we derived V γ as the average of the fit and Amp(RV) as the difference between the maximum and the minimum of the fit. The estimates of V γ , Amp(RV) and their uncertainties are provided in Tables 8 and 9. Note that we provide these estimates for both the Balmer lines and for the averaged RVCs of Fe, Na and Mg. By using V γ and Amp(RV), we normalized all the RVCs of the TS RRLs and derived the Normalized RVCs (NRVCs).

Period bins for the RVC templates
The shape of both light curves and RVCs of RRLs depends not only on the pulsation mode but also on the pulsation period. To improve the accuracy of the RVC templates available in the literature, we provide independent RVC templates for RRc and RRab variables. Moreover we divide, for the first time, the typical period range covered by RRab variables into three different period bins. This improvement is strongly required by the substantial variation in pulsation amplitudes (roughly a factor of five) when moving from the blue to the red edge of the fundamental instability strip (Bono & Cassisi 1999) and in the morphology of both light curves and RVCs (Bono et al. 2020a;Braga et al. 2020).
With this in mind, we adopted the same period bins that were used for the NIR light curve templates in (Braga et al. 2019) and for the optical light curve templates in Section 2. The reasons for the selection of these specific thresholds were already discussed in Braga et al. (2019): they are basically to maximize the number of points per bin without disregarding the change of the curve shape with period, and to separate RRLs with/without Blazhko modulations. The same arguments hold for the current investigation, with the additional advantage that, by adopting the same bins, the whole set of RVC templates and NIR light curves  templates are rooted on homologous sub-samples of RRL variables. The Bailey diagram and its velocity amplitude counterpart in Fig. 4.4 show their Amp(V), Amp(RV), pulsation periods, and adopted period bins. We provide 28 RVC templates in total, by considering the combination of seven different spectroscopic diagnostics (H α , H β , H γ , H δ , Mg, Na, and Fe+Sr, see Table 4) and four period.

The reference epochs of the RVC templates
The photometric data available in the literature for the TS RRLs were not collected close in time with our spectroscopic data. Therefore, we cannot anchor the RVC templates to the photometric reference epochs (e.g., T opt mean ). Small period variations and/or random phase shifts might significantly increase the dispersion of the points in the cumulative RVCs. The phase coverage of the TS RRLs is good enough to provide independent estimates of both the pulsation period and of the reference epoch. Moreover, T opt mean matches T RV(Fe) mean within 5% of the pulsation cycle (see Section 6 for more details), therefore we can adopt the latter to compute the RVC templates of the metallic lines (Fe, Mg, Na). This choice allows anyone to adopt T opt mean to phase the RV measurements and then use our templates (see Appendix C for detailed instructions). Note that, to compute the RVC templates of the Balmer lines, we use T RV(Hβ) mean because there is a well defined difference in phase between T RV(Hβ) mean and T RV(Fe) mean . To provide a solid proof of our assumptions, we performed the same test discussed in Section 2.1. Table 9. Barycentric radial velocities and RV amplitudes based on metallic lines. We derived T RV mean and T RV min from the average of the RVCs for both the Fe group lines and the H β line, taken as representative of the Balmer lines. As expected, T RV min matches in first approximation T opt max . Indeed, the latter was adopted by Liu (1991) and by Sesar (2012) to anchor their RVC templates. The working hypothesis behind this assumption is that the minimum in the RVC of the metallic lines takes place at the same phases at which the RVCs based on the Balmer lines attain their minimum, i.e. that T RV(Fe) min matches T RV(Hβ) min . However, we checked this assumption and found that the mean difference in phase between the two epochs is 0.036±0.051, with the minimum in the radial velocity curves of Fe group lines leading the H β minimum. Although the difference is consistent with being zero, its standard deviation is not negligible: a systematic phase drift of one twentieth around the minimum of the RVC might lead to offsets in the estimate of V γ of the order of 10 km/s. We point out that TS RRLs have multiple estimates V γ and of Amp(RV)-three from the metallic lines plus four from individual Balmer lines-but they only have two reference epochs: one for the Fe group lines, representative of the metallic RVCs and one for H β . The individual estimates of the reference epochs for the two groups of lines are listed in Table 6.
The basic idea is to have the different RVC templates phased at reference epochs originating from similar physical conditions. Weak metallic lines and strong Balmer lines display a well-defined RV gradient and their RVCs are also affected by a phase shift since the former form at high optical depths, and the latter at low optical depths (Liu & Janes 1990;Carney et al. 1992; Bono et al. 1994). We verified that the reference epochs for the weak metallic lines (Fe, Mg, Na) are the same within ∼3% of the pulsation cycle, while those for the Balmer lines, anchored to the H β RVC, are the same within ∼5% of the pulsation cycle. Figures 7 and 8 display the cumulative and normalized RVCs (CNRVCs) based on the Fe group lines and on the H β line, respectively. Data plotted in these figures display two interesting features worth being discussed in detail. i)-The residuals between observations and analytical fits for the CNRVCs based on the Fe group lines and phased using the reference epoch anchored to T RV(Fe) mean are systematically smaller than the residuals of the same CNRVCs anchored to T RV(Fe) min . The difference ranges from ∼30% for TS RRLs in the period bin RRab1 to ∼40% for TS RRLs in the period bin RRab3. ii)-The impact of the two different reference epochs is even more evident for the Balmer lines (Fig. 8). Indeed, the difference in the standard deviation ranges from ∼15% in the period bin RRc to ∼45% for the period bin RRab2. Moreover, the CNRVCs for the RRc and the RRab3 period bin show quite clearly that the reference epoch (vertical dotted line) anchored to τ 0 =T RV(Fe) mean takes place at phases that are slightly earlier than the actual minimum (see right panels). This difficulty is associated with the shape of both light curves and RVCs and it causes larger and asymmetrical residuals when compared with the CNRVCs of the same period bin anchored to the mean magnitude/systemic velocity (see the histograms plotted in the panels to the right of the CNRVCs).
The current circumstantial evidence indicates that RVC templates based on Fe group RV measurements and anchored to τ 0 =T RV mean provide V γ that are on average ∼30% more accurate than the same RVCs anchored to τ 0 =T RV min . The improvement in using T RV mean compared with T RV min becomes even more relevant in dealing with the RVCs based on Balmer lines. Indeed, uncertainties are smaller by up to a factor of three (see Section 7). This further supports the use of T RV mean as the optimal reference epoch to construct RVC templates.

RADIAL VELOCITY CURVE TEMPLATES
Before deriving the RVC templates, two pending issues need to be addressed: are the RVCs for the different lines in the Fe group, in the Mg Ib triplet and in the Na doublet, within the errors, the same? Are the mean RVCs of these three groups of lines affected by possible systematics?

Mean RVC templates for the three different groups of metallic lines
Our dataset is large enough to derive RVC templates for each single absorption line listed in Table 4. However, our goal is to provide RVC templates that can be adopted as widely as possible. Therefore, we aim to provide one RVC template for each of the Balmer lines and one for each of the metallic groups (Fe, Na, Mg), making a total of seven different sets of RVC templates. This is feasible only if the RVCs of the lines belonging to the same group have, within the errors, the same intrinsic features, i.e. the same shape, amplitude and phasing.
In principle, the RVC derived with a specific absorption line is different with respect to the one derived from any another line, because different lines may form in different physical conditions of the moving atmosphere. This is quite obvious for the Balmer series, with Amp(RV) progressively increasing by ∼60-70% from H δ to H α (S12, Bono et al. 2020a). This is the reason why independent RVC templates have to be provided for each Balmer line. However, for the Fe, Mg and Na groups, it is not a priori obvious whether different lines of the same group (e.g., Fe1 and Fe2) display, within the uncertainties, similar Amp(RV) and RVC shapes. Therefore, we verified whether the RVCs within the Fe, Mg and Na groups agree within uncertainties.
To investigate on a quantitative basis the difference, we inspected the residuals of the RVCs of each line with respect to the average RVC of the group. The left panels in Fig. 9 display the residuals of the RV measurements based on Fe1, Fe2, Fe3 and Sr lines with respect to the average of the four lines. The middle and right panels are the same, but for the two Mg and the two Na lines. Note that we discarded all the Mg b1 RV measurements because this line is blended with an Fe I line (5167.50 • A). This iron line can be as strong as the Mg b1 line itself or even stronger depending on the Mg abundance and on the effective temperature. Therefore, even using high-resolution spectra the velocity measurements with the Mg b1 line refer to an absorption feature with a center that changes across the pulsation cycle. Figure 9 displays, both quantitatively and qualitatively, that, for Mg and Na, there is no clear trend within the dispersion of the residuals. The maximum absolute offset is vanishingly small, being always smaller than 0.50 km/s, which is also smaller than our RV uncertainty. This means that these lines trace the dynamics of the same atmospheric layer and they can be averaged in order to derive a single set of RVC templates for both Mg and Na groups.
The same outcome does not apply to Fe and Sr lines. Indeed, although the average offsets are all smaller than the dispersions, they seem to follow a trend. More specifically, the average offset of the Fe1 curve is always smaller than the other, while the average offset of the Sr curve is typically larger. This is mostly due to the interplay of a small difference in Amp(RV) (generally increasing from Fe1 to Sr) and a mild trend in the average velocity (generally increasing from Fe1 to Sr). However, all these offsets are smaller than the dispersions (≤ 2.5 km/s) and similar to the intrinsic dispersion of the RVC templates (see Section 5.2). We also note that the standard deviation of the points around the offsets is, within the uncertainties, constant along the pulsation cycle. Indeed, a minimal increase around phase zero is only present for RV measurements based on iron in the period bin RRab1.
In the light of these results, we opted to derive three independent mean RVCs for the Fe, the Mg and the Na groups of lines. Selected RVCs for these three metallic diagnostics and for the four individual Balmer lines, anchored to τ 0 =T RV(Fe) min and τ 0 =T RV(Hβ) min , are shown in Figures 10 and 11.

Analytical fits of the RVC templates
To construct the RVC templates, we normalized the RVCs by subtracting V γ and dividing by Amp(V r ). Once normalized, we stacked the RVCs within the same period bin of RVC template, thus obtaining the CNRVCs (see Appendix B).
To provide robust RVC templates we decided to fit the CNRVCs with an analytical function. We discarded the Fourier series as a fitting curve because the number of phase points is not large enough (at least in the RRab3 period bin) to avoid non-physical bumps and spurious secondary features in the fits. We adopted the PEGASUS fit as for the light curve templates in Section 2.   Second row: Same as the top, but for the RRab1 variable HH Pup. Third row: Same as the top, but for the RRab2 variable RV Oct. Fourth row: Same as the top, but for the RRab3 variable AT Ser.
Note-We will provide this table only after the paper will be officially published on ApJ. Table 10 provides the coefficients of the PEGASUS functions obtained from the fitting procedure. When possible, we favored the lowest possible order, especially for the RRc template, as first overtone pulsators have more sinusoidal RVCs. These are also the coefficients for the analytical form of the RVC templates. Figures 16 and 17 display the analytical fits of the RVC templates together with the observed RV measurements.
The largest standard deviations are those for the H γ and H δ RVC templates for the RRc period bin (∼0.08 and 0.13, respectively). To convert the σ into the uncertainty on V γ , one has simply to factor in Amp(RV). Since the typical Amp(RV) for RRc variables in H γ and H δ range between 10 and 30 km/s (Bono et al. 2020b), the largest possible uncertainty introduced by the RVC template is of ∼4 km/s. However, the largest absolute uncertainties are associated with the H α and H β RVC templates for RRab1 period bin (σ ∼0.05). Since Amp(RV) is much larger for these diagnostics, the absolute uncertainties on V γ based on these RVC templates are of ∼7 km/s. For metallic lines, all these uncertainties are on average smaller than ∼3 km/s. These results concerning both metallic and Balmer lines indicate that the current RVC templates can provide V γ for typical Halo RRLs with an accuracy better than 1-3%.

REFERENCE EPOCH TO APPLY THE RADIAL VELOCITY CURVE TEMPLATES
A crucial aspect of templates is that they are used especially when the number of RV measurements is small. A first consequence is that, in a realistic scenario, the RV data is insufficient to accurately estimate T RV mean for either metal or Balmer RVCs. However, optical photometry usually is conducted before spectroscopic observations, and a good knowledge of the pulsation period and of Amp(V) are required to apply the RVC template. This means that we are typically dealing with a fairly well sampled V-band light curve, and in turn, with an accurate estimate of T opt mean . With this in mind, it is necessary to check whether T RV mean and T opt mean in RRLs take place, within the errors, at the same phase along the pulsation cycle. If this is the case, we could safely use of T opt mean to phase the spectroscopic data and to apply the RVC template.
We phased a subset of RRLs with a good sampling of the pulsation cycle both in the V-band and in metallic RVCs. Fortunately, among the RRLs of the TS, there are several objects that were used for the Baade-Wesselink (BW) analysis and for which there are available V-band light curves and RVCs collected at relatively close epochs (at most within ∼3 years). This is an important advantage for this consistency test because possible changes in phases (phase drifts) and the effect of period derivatives are small. We derived both T RV(Fe) mean and T opt mean and they are listed in Table 11.
Data listed in this table clearly show that T RV(Fe) mean and T opt mean trace, within the errors, the same phase along the pulsation cycle. Indeed the average phase difference is 0.007±0.019 and always smaller than 0.05 (see column 6 in Table 11) 1 This means the two reference epochs provide the same phasing. As a consequence, the photometric T opt mean can be safely adopted to anchor the RVC templates.
We already mentioned that there is a difference in phase between T RV(Fe) mean and T RV(Hβ) mean . This means that, when adopting T opt mean to use the template on Balmer lines, it is necessary to first shift the phases by an offset ∆τ Fe Hβ =Φ(T RV(Fe) mean -T RV(Hβ) mean ). For this reason, we adopted the data listed in Table 6 and found a linear trend of ∆τ Fe Hβ as a function of the pulsation period (see Fig. 12). The plausibility of the phase difference between metallic and Balmer lines is further supported by the empirical evidence that the standard deviation of the relation is vanishing (0.008). Indeed, it is almost one order of magnitude smaller than the standard deviation of the phase offset between T RV(Fe) min and T RV(Hβ) min (see Section 4.5). Large photometric surveys, however, often provide T opt max but not T opt mean . To overcome this limitation and to facilitate the use of the RVC templates, we provide relations for ∆Φ in Section 2.2 that allow the template user to easily convert the phases anchored on τ 0 =T opt max can into phases anchored on τ 0 =T opt mean .

VALIDATION OF THE RADIAL VELOCITY CURVE TEMPLATES
A solid validation of the RVC templates requires that photometric and radial velocity measurements be as close as possible in time. This methodological approach provides the unique opportunity to derive accurate photometric (pulsation period, V-band mean magnitude, Amp(V), reference epoch) and spectroscopic (V γ ) properties of the RRLs adopted for the validation. We have already mentioned in Section 4.5 that the temporal proximity of both photometric and spectroscopic data is only available for a small number of field RRLs. To overcome this limitation we decided to select from the calibrating sample one RRL per period Note-In column 7, the references for the RVC and light curves adopted to derive ∆Φ are listed in the following order: bin: YZ Cap for the RRc, V Ind for the RRab1, W Crt for the RRab2 and AT Ser for the RRab3. We label these four RRLs as the Template Validation Sample (TVS) and they are the benchmark for the RVC template validation.
Ideally, the validation should be done with an RRL sample independent from the one adopted to construct the RVC templates. However, we have verified that, by removing the four TVS RRLs from the calibrating RRLs, the coefficients of the analytical fits are only minimally affected. Note that the selection of the validating variable might bias the uncertainties of the result because there is a small degree of internal variation in the shape of the RVCs of the different period bins. To investigate on a quantitative basis the dependence on the validating variable, we performed several tests by using as validating variable one after the other all the RRLs included in each period bin. Interestingly enough, we found that the medians (see Tables 12 and 13) agree within 1σ, while the variation of the dispersions is on average smaller than 20%. This means that the selection of the validating variable has a minimal impact on the accuracy of the validation. We adopted the V γ of the TVS RRLs obtained by directly fitting their RVCs (see Tables 8 and 9) and assumed these as the best estimates for the systemic velocity (V γ(best) ) to be used in the validation process. However, to use the RVC template we need to convert Amp(V) into Amp(RV) and then use the latter to rescale the normalized analytical function. The ratio between Amp(V) and Amp(RV) together with the equations for the conversion are thoroughly discussed in Appendix A.
The analytical form of the RVC templates can be used both as curves to be anchored to a single phase point and as functions to fit a small number of phase points (three or more). Therefore, we followed two different paths for the validation process, based either on a single phase point approach (Section 7.1) or on multiple phase points (Section 7.2). The key idea is to estimate the accuracy of the light-curve templates from the difference ∆V γ between V γ(best) and the systemic velocities estimated by adopting the RVC template (V γ(templ) ). Furthermore, the accuracy of the current RVC templates is quantitatively compared, using TVS RRLs, with similar RVC templates available in the literature (S12).

Single phase point
As a first step, we generated a grid of 100 phase points for the seven RVCs (Fe, Mg, Na, H α , H β , H γ , H δ , that we label, respectively, with a j index running from 1 to 7). We interpolated the analytical fits of the RVCs for the TVS RRLs on an evenlyspaced grid of phases (φ i =[0.00, 0.01, ... 0.99], where i runs from 1 to 100). For each φ i and each RVC( j), we generated a RV measurement with the sum RV i j = RV( f it(φ i )) j + rσ. The two addenda are: i) RV( f it(φ i )) j , which is the value of the fit of the RVC( j) interpolated at the phase φ i ; ii) rσ, which simulates random noise: σ is the standard deviation of the phase points around the fit and r is a random number extracted from a standard normal distribution. Figure 13. Left: Black crosses show the Fe RV measurements that we generated from the grid for the RRc variable YZ Cap. Gray dashed lines display the RVC templates from this work (labeled as TW) associated with each phase point. The ID of the RRL is labeled. Note that the comparison with the RVC template provided by S12 is missing because YZ Cap is an RRc variable. Middle: Same as the left, but for the RRab variable V Ind belonging to the RRab1 period bin. Top: Radial velocity curves based on the RVC templates provided by S12. Bottom: same as the top, but based on our RVC templates. Right: Same as the middle, but for the RRab variable SX For belonging to the RRab2 period bin.
We applied the RVC templates on each of the simulated points, thus deriving 100 estimates of the systemic velocity (V templ γ(i, j) ) for each RVC( j). To provide a quantitative comparison, these estimates were performed using both our own and S12 RVC templates. The S12 RVC templates were not applied to YZ Cap, since they are valid only for RRab variables. In discussing the difference between our own RVC templates and those provided by S12, there are three key points worth being mentioned: i)-the comparison for Fe, Mg and Na RVCs was performed by using the generic metallic RVC template from S12, because they do not  Note-Medians (mdn) and standard deviations (σ) of for the ∆V γ based on a single phase point validation, for both our and S12 RVC templates.
provide specific-element RVC templates; ii)-the RVC templates by S12 were anchored to T opt max , while our own were anchored to T opt mean ; iii)-the RVC templates by S12 were rescaled, for internal consistency, by using their conversion equations from Amp(V) to Amp(RV) and not our own equation for the amplitude ratio.
Note that, due to the paucity of RRab3 variables, it was not possible to find one with both a reliable estimate of the reference epoch and with RV measurements close in epoch to the optical light curve. Therefore, the RVC templates for the RRab3 period bin were not validated with the single phase point approach. Nonetheless, we anticipate that we successfully validated the RRab3 RVC templates with the three-point approach (see section 7.2). Figure 13, displays the simulated RV(Fe) phase points and the RVC templates applied to them. Finally, we derived the offsets ∆V γ(i, j) = V templ γ(i, j) − V best γ( j) for each point and each RVC template. Table 12 gives the median and standard deviations of ∆V γ for each RVC template.
We note that that the median of the ∆V is always smaller, in absolute value, than 1.0 km/s and 1.5 km/s for metallic and Balmer RVC templates. In all cases, the standard deviations are larger than the residuals, meaning that the latter can be considered vanishing within the dispersion. The largest standard deviations are found for the H α and H β RVC templates and progressively decrease when moving to H γ , H δ and metallic lines.
The comparison between the current and the RVC templates provided by S12 brings forward that the standard deviations of the former ones are systematically smaller by a factor ranging from ∼1.5 to ∼3. The higher accuracy of the current RVC templates is due to an interplay of a more robust estimate of T mean with respect to T max and a more accurate optical-to-RV amplitude conversion (note e.g., in the upper-right panel of Fig. 13, Amp(RV) is clearly overestimated for the S12 RVC template).

Multiple phase points
To apply the RVC templates to single phase points, it is necessary to know four parameters with great accuracy: i) the pulsation period, ii) the pulsation mode, iii) Amp(V) and iv) the reference epoch of the anchor point (T opt mean ). The last one is particularly delicate, not only because a good sampling of the light curve is needed but also because, when the spectroscopic data were collected several years before/after the photometric data, even small phase shifts or period changes can affect the phasing of the RV measurements. Note that for RRLs affected by large Blazhko modulations, these parameters-especially Amp(V) and reference epoch-cannot be accurately estimated. Therefore we suggest to avoid the application of the templates to RRLs with evident Blazhko modulations.
To overcome this limitation, it is possible to use the RVC templates as fitting functions if at least three RV measurements are available. In this empirical framework, only three parameters are required in order to apply the RVC template, namely the pulsation period, the pulsation mode and Amp(V). We performed a number of tests by assuming that three independent RV measurements were available. In this approach, the Amp(V) has to be converted into the Amp(RV) and then perform a leastsquares fit of the RV measurements by adopting the RVC template as the fitting function. The minimization of the χ 2 is based on two parameters: ∆φ (a horizontal shift) and ∆V γ (a vertical shift). The function to be minimized is: Table 13. Results of the validation for the multiple phase point approach. Note-Medians (mdn) and standard deviations (σ) for the ∆V γ obtained from the three-points validation, for both our and S12 RVC templates.
To simulate three RV measurements, we extracted three random phases over the pulsation cycle and generated three RV measurements following the same approach adopted in Section 7.1. To rely on a wide statistical sample, the process was repeated 5000 times to generate 5000 different triplets of RV measurements for each RVC template. Henceforth, we label each of the triplets with a subscript k.
We estimated the systemic velocity V templ γ(k, j) by fitting each of the triplets with both S12 and our own analytical form of the RVC templates. Analogously to Section 7.1, we derived the offsets ∆V γ(k, j) and their median and standard deviations and they are listed in Table 13.
We note that the uncertainties are larger for the three-point approach with respect to the single point one discussed in the former section. This happens because, when the three points are randomly extracted, it may happen that two (or all) of them are too close in phase, and the fitting procedure provides less solid estimates. This is a realistic case in which one might not have control over the sampling of the spectroscopic data. Also, the very fact that we are fitting the shift in phase, and not anchoring the RV data to the true reference epoch of the variable, means that the phasing of the template is not fixed, and this naturally leads to larger uncertainties. We have verified that, if the reference epochs were available, we would obtain standard deviations that are ∼10-30% smaller and medians that are up to ∼50% smaller than in the case of the single-epoch approach.
Note that we did not put any restriction on which RV measurements were chosen to form a triplet. More specifically, we did not remove triplets in which two phase points were very close in phase, thus mimicking a realistic situation where either only two RV independent measurements were obtained or where the points are really close in phase due to the sampling of the spectroscopic data. Although this scenario usually leads to flawed estimates of V γ , they are a minority, with less than ∼10% of the triplets having points closer than 0.05 in phase. Moreover, it is not only the difference in phase that matters but also the distribution along the pulsation cycle. Indeed, phase points close in phase located along the decreasing branch of RRab variables produce larger uncertainties when compared with other phase intervals. The decision to keep even these troublesome triplets in our tests allowed the computed errors to take into account all the possible drawbacks of real observations, including the scenario of spectroscopic surveys that scan the sky by taking multiple consecutive measurements.
Data listed in Table 13 show that, with the only exception of the RRab3, for which the standard deviations are similar, our RVC templates provide smaller standard deviations of the ∆V compared to S12 templates. This is true even for the three phase points validation, where the median offsets are smaller than the standard deviations. This means that the RVC template provides V γ estimates that are more accurate than the simple average of the three RV measurements.  Piersimoni et al. 2002;Layden & Sarajedini 2003;Arellano Ferro et al. 2014) This cluster has a very high RV: 494.5±0.4 km/s based on giant and subgiant stars and 496.47±0.11 km/s based on turn-off, subgiant and red giant stars (Ferraro et al. 2018;Wan et al. 2021, respectively). Note that Magurno et al. (2018) adopted the S12 RVC templates and applied them to eleven RRLs in NGC 3201 obtaining a cluster RV of 494±2±8 km/s which is within 1σ of literature estimates. Interestingly, by taking advantage of accurate proper motion measurements based on Gaia DR2 Massari et al. (2019) suggested that NGC 3201 is likely an object accreted during either the Sequoia or the Gaia-Enceladus merger.

RR LYRAE IN NGC 3201 AS A TEST CASE
To investigate on a more quantitative basis the cluster properties, we used the MUSE@VLT spectra collected by Giesers et al. (2019) for seven cluster RRab variables. From these spectra, we measured H α and H β RVs (see data listed in Table 15) and the RVCs are displayed in Figure 14. Unfortunately, the optical light curves that we adopted from Arellano Ferro et al. (2014) do not fully cover the pulsation for all the cluster RRLs. More specifically, we could not derive a reliable estimate of the reference epochs for V5 and V100. Due to this, we are not able to apply the templates to these two stars. In principle, we could apply the three-point method, but we wanted to keep our V γ estimates as homogeneous as possible.
After deriving-with the PLOESS algorithm-the best estimate for the systemic velocity (V γ(best) , displayed in Table 14) by fitting the RVCs, we applied the RVC template by using the single phase point approach. Note that, in this case, we do not generate a grid of synthetic points, but we simply extract the points one by one from the RVC. Table 16 shows the average V γ(templ) estimates for all the lines, together with the standard deviation and the median of the offsets. The results are similar to those found in Section 7.1, meaning that all the medians are smaller than the standard deviations and, on average, the standard deviations obtained by using our own RVC templates are smaller than those coming from the use of the S12 RVC templates.
Our final estimate for the velocity of the whole system, by using both H α and H β is 496.89±8.37 (error) ±3.43 (standard deviation) km/s by using the RVC templates and 495.21±3.23 (error) ±4.32 (standard deviation) km/s by simply fitting the RVCs, which agree quite well, within 1σ, with similar estimates available in the literature.

SUMMARY AND FINAL REMARKS
We provide accurate and homogeneous RVC templates for both RRab and RRc variables using for the first time spectroscopic diagnostics based on well-identified metallic (Fe, Mg, Na) lines and on Balmer (H α , H β , H γ , H δ ) lines. In the following we discuss the approach adopted to construct the RVC templates and summarize the most relevant results.
V-band light curve templates -To demonstrate on a quantitative basis the difference between the reference epoch anchored to the maximum brightness and to the mean magnitude along the rising branch, we collected homogeneous V-band photometry for cluster (ω Cen , M4) and field RRLs. We have grouped them into four period bins (the same that we have used for the RVC templates) and found that the anchoring to the epoch of the mean magnitude on the rising branch (T opt mean ) provides smaller standard deviations on the light curve templates than the anchoring to the maximum brightness (T opt max ). The decrease is of the order of 35% for the period bins with cuspy light curves (RRab1, RRab2) and of 45% for those with flat-toopped light curves. These finding strongly supports the results obtained by Inno et al. (2015) and by Braga et al. (2019) in using T opt mean to phase the NIR light curves of both classical Cepheids and RRLs.
Spectroscopic catalog -In this work, we present the largest collection of RRL spectra ever compiled in the literature. Altogether, we collected 23,865 spectra for 7,070 RRab and 3,343 RRc variables. These measurements were secured using eleven different spectrographs, ranging from low (2,000) to very high spectral resolution (115,000). To build the RVC templates, the  Table 14. For V5 and V100, we have assumed an arbitrary reference epoch, since we could not derive an accurate estimate from the optical light curves. most important dataset is the du Pont (see Table 3). Spectroscopic observations at this telescope were specifically planned to cover the entire pulsation cycle of several bright RRLs (< V >∼10-11 mag for the majority of them). We also collected RV measurements and V-band time series from the literature (the Baade-Wesselink [BW] sample), which were crucial to investigate the reference epoch and the ratio between RV and optical amplitudes.
Amplitude ratio -To apply the RVC templates, it is necessary to have prior knowledge on the optical amplitude of the variable, to correctly rescale the RVC template itself and to optimally match the true RVC. We provide new equations for the ratio of Amp(RV) over Amp(V). Those available in the literature (S12), for RVCs based on metallic and Balmer (H α , H β , H γ ) lines, were constructed using six RRab variables covering a very narrow range in pulsation period (0.56-0.59 days). In this investigation we provide different RVC templates for both RRc and RRab variables based on metallic and Balmer (H α , H β , H γ , H δ ) lines. Even more importantly, our relations are based on three dozen variables covering a wide range in pulsation periods (0.27-0.83 days) and metal content (-2.6 [Fe/H] 0.0).
Reference epoch -When applying the RVC template to single RV measurements, it is necessary to anchor the RVC template at the same epoch of the observations. The RVC templates are applied to RRLs with well sampled optical light curves and a few spectroscopic measurements. Therefore, the only pragmatical possibility to phase the spectroscopic data is to derive a reference epoch from the light curve. By using the light curves and the RV curves of our BW sample, we demonstrated that T opt mean is, within 5% of the pulsation cycle, identical to the time of mean velocity on the decreasing branch of the Fe RVC (T RV(Fe) mean ). This means that RVC templates based on Fe, Mg and Na lines can be safely anchored to T RV(Fe) mean , and this approach does not introduce any systematic error when T opt mean is adopted to apply the RVC template. We also found that RVCs based on Balmer lines display a well defined phase lag across the phases of mean RV. Therefore, we decided to anchor the current H α , H β , H γ and H δ RVC templates to the time of mean velocity on the decreasing branch of the H β RVC (T RV(Hβ) mean ), taken as representative of the Balmer lines. Additionally, we found a clear trend in the phase difference between T RV(Fe) mean and T RV(Hβ) mean as a function of the pulsation period. The new analytical relation gives the phase difference required to use the Balmer RVC templates, because T opt mean does not match T RV(Hβ) mean . To discuss the concerning pros and cons of the reference epochs anchored to T RV mean and to T RV min into a more quantitative framework, we performed a series of numerical experiments. Interestingly enough, we found that RVC templates based on metallic RV measurements and anchored to τ 0 =T RV mean provide V γ that are on average ∼30% more accurate than the same RVCs anchored to τ 0 =T RV min . Even more importantly, we found that the improvement in using T RV mean compared with T RV min becomes even more relevant in dealing with RVCs based on Balmer lines. Indeed, the uncertainties decrease up to a factor of three (see Section 7). The current circumstantial evidence further supports not only the use of two independent reference epochs for metallic and Balmer lines, but also the use of T RV mean as the optimal reference epoch to construct RVC templates. Finally, we investigated the correlation in phase between T opt mean (Φ mean ) and T opt max (Φ max ). This search was motivated by the fact that the current photometric surveys only provide reference epochs anchored to the time of maximum light (T opt max ). We found that the difference between the two (∆Φ max−min = Φ max − Φ min ) is constant-within the dispersion-for RRc variables and it follows a linear trend for RRab variables. We provide these relations for those interested in using the current RVC templates to RRLs for which only T opt max reference epochs are available. Radial velocity curve templates -We provide a total of 28 RVC templates of RRLs: these are divided into four different period bins (one for the RRc and three for the RRab variables) and seven diagnostics (Fe, Mg, Na, H α , H β , H γ and H δ ). The analytical form of the templates is provided in the form of a PEGASUS series (fifth to ninth order). Analytical fits based on PEGASUS series, when compared with Fourier series, do not display spurious ripples when there are neither gaps in the phase coverage nor in the case of uneven sampling. The RVC templates have intrinsic dispersions that lead to errors smaller than 10 km/s in the worst case (H α and H β for high amplitude RRLs) and one order of magnitude smaller for the RVC templates with the smaller intrinsic dispersion (metallic lines RVC templates). To maximize the reach of the results of this work, we provide, in Appendix C, clear instructions on how to apply the RVC templates in different real-life observational scenarios.
Template validation -To validate the current RVC templates we performed a detailed comparison with the RVC templates provided by S12 for RRab variables and based on metallic and Balmer (H α , H β , H γ ) RVCs. We performed these tests on a sub-sample of variables that were used to build the RVC template (YZ Cap, V Ind, W Tuc, AT Ser). The validation process was performed using both a single phase point approach, where the knowledge of the reference epoch is mandatory, and with a three phase points approach, where the RVC template is used as a fitting function. We found that the median offset of the V γ estimates from the RVC templates is always smaller than 1.5 km/s (one point approach) and 6 km/s (three point approach). The medians are smaller than the standard deviations, meaning that systematic errors are negligible with respect to the statistical errors. We also found that our RVC templates provide V γ estimates that have a dispersion smaller by a factor of 1.5-3 than those based on the RVC templates provided by S12.
RRLs in NGC 3201 -We reduced the MUSE spectra already presented in Giesers et al. (2019) and obtained H α and H β RVCs. We derived V γ both by fitting the RVC and by extracting the measurements one by one, and by adopting the RVC templates. Their results based on these RV measurements are very similar to those of the validation process, with offsets smaller than 6 km/s and standard deviations that are smaller than those on the S12 RVC templates. Our estimate of the V γ of the cluster is 496.9±8.4 (error) ±3.4 (standard deviation) km/s from the RVC templates and 495.2±3.2 (error) ±4.3 (standard deviation) km/s from the fit with the RVC templates. They both agree, within 1σ, with literature estimates.
In the next coming years the ongoing (OGLE, ASAS-SN, ZTF, sky-mapper, PanSTARRS) and near future (VRO) ground-based and space-based (Gaia, WISE, WFIRST) observing facilities will provide a complete census of evolved variable stars associated with old stellar tracers in the Milky Way and in Local Group galaxies. This will open new paths in the analysis of the early formation and evolution of the Galactic spheroid. However, firm constraints on the formation mechanism, namely, the dissipative collapse (Eggen 1962), dissipation-less mechanism (Searle & Zinn 1978) and Cold Dark Matter cosmological models (Monachesi et al. 2019) require detailed information concerning the kinematics and the metallicity distribution of the adopted stellar tracers. This is the approach already adopted to fully characterize the stellar streams and the merging history of the Galactic Halo (Helmi et al. 2018;Vasiliev 2019;Prudil et al. 2021).
Upcoming and ongoing low-(LAMOST-LR, Su et al. 1998 Su et al. 1998;RAVE, Steinmetz et al. 2006;WEAVE, Dalton et al. 2012) and high-resolution (APOGEE, Majewski et al. 2017;MOONS, Taylor et al. 2018;PFS, Tamura et al. 2018) spectroscopic surveys will provide a wealth of new data for large samples of dwarf and giant field stars. In this context, old variable stars play a crucial role because their individual distances can be estimated with a precision of the order of 3-5% within the Local Group. Recent spectroscopic investigations based on high resolution spectroscopy Sneden et al. 2017;Crestani et al. 2021a) indicate that detailed abundance analysis can be performed with spectra collected at random phases. Unfortunately, the typical limiting magnitudes for spectroscopic investigations are roughly five magnitudes brighter than photometric ones, with current spectrographs available at the 8-10m class telescope allowing us to reach limiting magnitudes of the order of V∼20-21 mag. However, even with the usage of large telescopes, the estimate of V γ is time consuming, because it requires a spectroscopic time series of at least a dozen points. The RVC templates developed in this investigation provide new solid diagnostics to provide accurate V γ determinations by using a small number (three or less) of RV measurements, based on low-resolution spectra. Highly accurate estimates of line-of-sight velocities of stream stars are imperative for constraining the dark matter distribution (Bonaca & Hogg 2018); these RVC templates will provide just that for the numerous RRL harbored by MW streams. The current diagnostics are focused on spectroscopic features located either in the blue or in the visual wavelength regime. A further extension into the red and the near-infrared regime is mandatory to fully exploit the most advanced space (Gaia, Gaia Collaboration et al. 2016;JWST, Gardner et al. 2006;Roman Spergel et al. 2015) and ground-based (GMT, ELT, TMT) observing facilities. mounting empirical evidence that amplitude ratios using NIR and MIR bands follow a mild trend with metallicity (Mullen et al. 2021). Before deriving the relation between Amp(RV) and Amp(V), we checked whether we could assume the Fe, Mg and Na amplitudes as equivalent. Figure 15 shows the residuals of Na and Mg Amp(RV) versus the Fe Amp(RV). It is clear that there is no trend with period and also the offset is well within the Amp(RV) uncertainty. This means that we can provide only two Amp(RV) vs Amp(V) relations (one for RRc and one for RRab variables) that hold for both Fe, Mg and Na. We also collected RV curves of RRLs from literature, all derived from metallic lines (BW sample, see caption of Table 11). We found that the BW Amp(RV) displays a very small difference with our Fe Amp(RV) (0.23±3.73 km/s). This allowed us to merge the two sets of Amp(RV) and derive a more solid relation for the metallic Amp(RV), based on a larger number of objects (12 RRc and 60 RRab).
The right panels of Figure 15 display the trend of Fe and of the Balmer Amp(RV) with Amp(V). A steady increase of the slopes for the Balmer Amp(RV), from δ to α is clear enough, both visually in the figure and quantitatively from the coefficients listed in Table 17. This is also the first time that different trends have been found for RRc and RRab. By inspecting the CNRVCs it is clear that, within the bins RRab1 and RRab2, the morphology of the RVCs is, on a first approximation, dichotomic. More specifically, within the RRab1 bin, SW Aqr, ST Leo, VX Her, V Ind and V0440 Sgr display a local maximum around the phase 0.70-0.75, instead of the more or less steady rising behavior of the other RRab1 RRLs. Moreover, TY Gru, CD Vel and SX For (RRab2 bin), do not display a local minimum around phase 0.7-0.8, as the other RRab2 variables do. This happens for all the diagnostics, although it is more evident for the Balmer lines. We remark that these features are also present in the optical light curves of these stars.
We checked whether these features can be associated with either pulsation or physical properties of the stars. While it is true that the RRab1 variables with a more prominent local maximum are located in the High-Amplitude Short-Period (HASP, Fiorentino et al. 2015) region, they are not the only HASP variables in our sample. Their iron abundances range between -1.9 and -1.4, which is around the peak of the distribution of field RRLs (Crestani et al. 2021a) and it was not even possible to constrain a morphological class of RVCs based on the Fourier parameters R2, R31, φ 21 and φ 31 of their light curve.
To sum up, there is no quantitative way to predict, either from the light curve or from the physical properties, which is the RVC morphology of RRab1 and RRab2. This has several consequences: the first and most obvious is that we cannot split these bins and provide more RVC templates because we cannot provide criteria for using one or the other. This may seem to be a disadvantage, but luckily, this dichotomy introduces a 5 km/s offset in the estimate of V γ with a probability of ∼10-15% (i.e., the fraction of pulsation cycle where the RVCs have a different behavior).

C. HOW TO USE THE NEW RVC TEMPLATES
In this section of the Appendices, we provide precise instructions on how to use the new RVC templates with the aim of estimating V γ in different realistic cases: i)-when only one RV measurement and a well-sampled light curve that allows the estimate of the reference epoch are available; ii)-when a few (we assume three) RV measurement and a well-sampled light curve that allows the estimate of the reference epoch are available; iii)-when a few (we assume three) RV measurement and a light curve that does not allow an accurate estimate of the reference epoch are available.  C.1. Estimate of T opt mean Before describing how to use the templates, we want to instruct the reader on how to derive T opt mean , which is not as common as T opt max , therefore it might not be straightforward to estimate. First of all, the light curve must be phased to an arbitrary reference epoch, e.g., HJD=0, as in Figure 18. Then, the phased light curve must be fit with a model. It is crucial that the model fits well the rising branch. After that, < V > is derived by converting each point of the model into arbitrary flux (F = 10 −0.4 * mag ), integrating the flux model and finding the mean flux < F > and converting back to magnitude (< V >= −2.5 · log 10 (< F >)).
If the model is analytical (e.g., Fourier or PEGASUS), one can easily find the phase on the rising branch at which the model intersects < V >, namely, φ mean . If the model is not analytical (e.g., spline or PLOESS), it is necessary to interpolate < V > with the model sampled on an even grid of phases. A convenient choice for the step of the grid is between 0.001 and 0.01.
Once φ mean is known, the next step is to select any phase point of the light curve. This will be characterized by an epoch of observation (t V(i) , where i indicates the i-th point of the light curve) and a phase (φ V(i) ). Finally, T opt mean = t V(i) − (φ V(i) − φ mean ) · P can be derived.

C.2. Single RV measurement
This is the most classical situation when using any type of template. In this case, only one RV measurement is available. Note that, if the spectral range of the instrument is large enough, it is possible to have one RV measurement per diagnostic (e.g., Fe, Mg and H γ ) but still no more than one epoch is available. This means that any RVC template can be applied to only one point.
The quoted RV measurement consists of an epoch of observation (t), a velocity (RV) and its uncertainty (eRV). In this case, a decent sampling of the light curve is needed, in order to estimate its period (P), amplitude (Amp(V)), mean magnitude (< V >) and the reference epoch T opt max or T opt mean . The first step is to anchor the RV measurement to the same reference epoch as the template. As demonstrated in the body of the paper, it is possible to assume T opt mean = T RV(Fe) mean and derive the phase as φ = t − T opt mean P mod 1. If only T opt max estimates are available, as is often the case of data releases of large surveys, it is necessary to derive the phase anchored to T opt max and then apply the offsets provided in Section 2.2: φ max = t − T opt max P mod 1 and then φ = φ max + 0.223 for RRc variables or φ = φ max + 0.043 + 0.099 · P for RRab variables. In this case, an uncertainty must be associated with φ, namely, the σ of the quoted relations: 0.036 for RRc and 0.024 for RRab. By using φ, the RV measurement is now anchored to the same reference of metallic (Fe, Mg, Na) templates. If, however, the RV measurement is from a Balmer line, it is necessary to convert φ by using the relation provided in Figure 12, that is φ Hβ = φ + 0.023 − 0.096 * P.
The second step is to rescale of the normalized template. For this end, it suffices to convert Amp(V) into Amp(RV) by using the relations provided in Table 17.
The third step is to derive the analytical form of the template both rescaled by Amp(RV) and shifted in zero-point to pass through the RV measurement. For this step, the right coefficients for the RVC template must be selected from Table 10. Afterwards, these coefficients are substituted into Equation 1 to calculate the value of the template (T(φ)) at the phase φ. Finally, V γ is obtained as V γ =RV − Amp(RV) · T (φ).
Of course, if RV measurements are available for more than one diagnostic, the quoted steps can be separately applied to all the diagnostics, providing multiple V γ estimates that can be averaged.

C.3. Multiple RV measurement with reference epoch
In this case, more than one RV measurement per diagnostic and a light curve enough well-sampled are available. The procedure is qualitatively identical to that described in C.2, but having more than one RV measurement per diagnostic allows the averaging of the V γ estimates for the same diagnostic. Also in this case, if RV measurements are available for more than one diagnostic, these can be averaged to constrain V γ on a broader statistical basis.
In principle, this method could be applied to any number of RV measurements. However, when these are ten-twelve (or more) and they are more or less evenly sampled, with a good knowledge of the period, it is possible to just fit the points and directly derive a V γ estimate as accurate as the template itself, or even more if the variable has experienced some phase drift or period change during the time elapsed from the collection of light curve and RV data.

C.4. Multiple RV measurement without reference epoch
In this case, more than two RV measurements per diagnostic are available, but the light curve is only modestly sampled. This means that period and Amp(V) can be estimated from photometric data, but the reference epoch cannot. In this situation, the templates can be used as fitting curves.
This approach is qualitatively different from the other two because one does not have to anchor the template but to fit it to the data. The steps are the following First of all, phases must be derived by adopting an arbitrary reference epoch T 0 (e.g., T 0 =0 or T 0 =2,400,000) φ = t − T 0 P mod 1. Still, if one wants to use the Balmer templates, the conversion φ Hβ = φ + 0.023 − 0.096 * P must be applied.
Secondly, Amp(V) has to be rescaled into Amp(RV) by using the relations provided in Table 17. This step is analogous to that described in C.2.
The third step, is the selection of the RVC template coefficients from Table 10 and perform a χ 2 minimization when fitting the RV data with Equation 2. The minimization must be performed on the two free paramaters ∆φ (a horizontal shift) and ∆V γ , while all the others remain fixed.
Finally, V γ is simply derived by integrating the final analytical form of Equation 2. Also in this case, if RV measurements are available for more than one diagnostic, these can be averaged to constrain V γ on a broader statistical basis. a Calculated at the given phase by using Equation 1 and the coefficients provided in Table 10