© The Authors 2024

1 Introduction

Pulsating variables of the RR Lyrae class are horizontal branch giants associated with old stellar populations (Catelan 2009; Savino et al. 2020) periodically changing their brightness with periods ranging from approximately 5 to 24 hours a day (Catelan & Smith 2015). They are one of the cornerstones for studies of the structure and dynamics of the Local Group (e.g., Layden et al. 1996; Clementini et al. 2001; Contreras Ramos et al. 2013; Jacyszyn-Dobrzeniecka et al. 2020; Prudil et al. 2020; Savino et al. 2022). They serve as distance indicators toward smaller stellar systems, such as globular clusters and dwarf galaxies (e.g., Martínez-Vázquez et al. 2016, 2019; Bhardwaj et al. 2021). The RR Lyrae stars are also used to trace tidally disrupted stellar systems often found in the Galactic halo (e.g., Hendel et al. 2018; Mateu et al. 2018; Koposov et al. 2019; Price-Whelan et al. 2019; Prudil et al. 2021) and can help investigate the shape and mass profile of the Milky Way halo (e.g., Medina et al. 2018, 2023; Erkal et al. 2019; Prudil et al. 2022). Their pulsation properties permit estimating distances with a precision better than 5 percent (e.g., Catelan et al. 2004; Braga et al. 2015; Neeley et al. 2017; Muraveva et al. 2018) and deriving photometric metallicity estimates for a population of RR Lyrae stars (e.g., Jurcsik & Kovacs 1996; Smolec 2005; Dékány et al. 2021).

The obvious advantages of utilizing the pulsation properties of RR Lyrae stars come at the cost of the intricacies of deriving information from stellar spectra compared to nonvariable stars. A given RR Lyrae star (depending on its pulsation amplitude) undergoes changes in its temperature (~1000 K), and surface gravity (~1 dex) (see, e.g., Cacciari et al. 1992; Skillen et al. 1993; For et al. 2011; Pancino et al. 2015). In addition, its observed line-of-sight velocity varies as well (~60 km s⁻¹, depending on the spectral lines used for estimating the line-of-sight velocity) within its single pulsation cycle (e.g., Clementini et al. 1990; Jeffery et al. 2007; Sesar 2012; Braga et al. 2021). This complicates the co-addition of nonconsecutively observed spectra, and the rapid changes in the line-of-sight velocity put constraints on the exposure time (longer exposures would result in asymmetric lines). Particularly, the variation in line-of-sight velocity where the pulsation of the atmosphere is convolved with the center-of-mass motion (from here on referred to as the “systemic velocity”) needs to be taken into account before one uses measured line-of-sight velocities of RR Lyrae for kinematical studies.

There are several approaches to estimating the systemic velocity of a given RR Lyrae star. First, one can sample the entire line-of-sight velocity variation and calculate the mean velocity. This method is costly with regard to the observation time and may not be possible for some RR Lyrae variables with long pulsation periods or periods of 0.5 days, especially with ground-based observations. Second, assuming that the pulsation ephemerides such as pulsation period and reference time are known, one can time the observing window for a phase around 0.4–0.5 where the line-of-sight velocity is close to the actual systemic velocity (Liu 1991; Sesar 2012). Although less costly in observation time, this method requires preparation before observation and the capability to modify the observation plan, which is often not possible for extensive spectroscopic surveys (de Jong et al. 2014; Dalton et al. 2014; Zasowski et al. 2017; Blanton et al. 2017; Kollmeier et al. 2017). The third option is to observe the star at any point in its pulsation cycle and then later shift the observation along a velocity template to obtain its systemic velocity. This method, similar to the second one, requires prior knowledge of the pulsation ephemerides and the pulsation amplitude as well (from photometric observations, often from the V passband; e.g., Liu 1991; Sesar 2012; Braga et al. 2021). Unlike the previous two methods, it does not constrain the spectroscopic observations. One can determine a systemic velocity from a single line-of-sight measurement at any point during the pulsation cycle.

The third approach seems to be the most used in the past decade (e.g., Sesar et al. 2013; Kunder et al. 2016, 2020; Liu et al. 2020; Hanke et al. 2020; Prudil et al. 2021; Medina et al. 2023). In addition to the pulsation ephemerides, the last method relies on the scaling relation between amplitudes of photometric and line-of-sight velocity curves and a line-of-sight velocity template that models the behavior of line-of-sight velocities during the pulsation cycle. From the previous studies, scaling relations and templates are available for the most prominent spectral hydrogen lines (Balmer lines Hα, Hβ, Hγ, Hδ) and metallic lines (e.g., Fe, Mg, etc.). In this work, we focus on expanding the scaling relations and line-of-sight velocity templates toward redder wavelengths, particularly toward the calcium triplet and H-band covered by the surveys described below. We selected these two regions since we aim to use newly derived scaling relations and templates to study the Galactic bulge, as most of the spectroscopic data for RR Lyrae stars cover this approximate redder wavelength range.

One of the two spectroscopic surveys used in this work is the Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et al. 2017; Wilson et al. 2019; Beaton et al. 2021). APOGEE is a large-scale infrared spectroscopic survey of the Milky Way focusing mainly on the highly reddened regions, such as the Galactic disk and Galactic bulge. The observations were conducted at the Apache Point Observatory (APO) in New Mexico, USA, and at the Las Campanas Observatory (LCO) in the Atacama Desert in Chile. The APOGEE survey covers the H passband (1.51–1.70 μm wavelength range) with a resolution of ~22 500 (SDSS; Eisenstein et al. 2011; Blanton et al. 2017). The main targets of the APOGEE survey are giant stars, for which APOGEE derived stellar atmospheric parameters, line-of-sight velocities, and chemical abundances for numerous elements (Nidever et al. 2015; Abdurro’uf et al. 2022). Among their targets toward the inner Galaxy, over 6000 RR Lyrae stars were observed, some of which have several epochs of spectroscopic measurements. Therefore, obtaining their systemic velocities could aid in studies of the Galactic bulge formation history.

The second survey we used is the Gaia space telescope mission focused on measuring objects’ precise positions and motions (Gaia Collaboration 2016, 2023). One of the instruments on board Gaia is the Radial Velocity Spectrometer (RVS), which aims to provide multi-epoch line-of-sight velocities. The RVS resolving power is R≈11 500, with a wavelength coverage ranging from 8450 Å to 8720 Å, thus covering the calcium triplet (Ca T) region (Cropper et al. 2018; Katz et al. 2023). The latest data release of the Gaia catalog provides individual line-of-sight velocities for 1096 RR Lyrae stars, with 1086 having more than seven observations (Clementini et al. 2023). The uncertainties on individual line-of-sight velocities range from ≈0.5 km s⁻¹ up to ≈85 km s⁻¹ (Clementini et al. 2023). This catalog provides exquisite data for calibrating scaling relations between photometric and line-of-sight velocity amplitudes and creating line-of-sight velocity templates. The tools mentioned above can be used for estimating systemic velocities for RR Lyrae stars observed by the Bulge RR Lyrae Radial Velocity Assay (BRAVA-RR, Kunder et al. 2016, 2020).

Unlike previous studies that provided linear (with both a slope and intercept) scaling relations between photometric and line-of-sight velocity amplitudes (e.g., Jones et al. 1988; Liu 1991; Sesar 2012; Braga et al. 2021), in our work the intercept of the scaling relations was set to zero. This conserves the physical representation of the relations and avoids deriving transformation equations that would predict line-of-sight velocity changes with zero photometric amplitude variation. Section 2 describes the derivation of the line-of-sight velocity templates and the scaling relations from the APOGEE spectra. Section 3 focuses on a similar task but for the Gaia RVS spectra that cover the Ca T region. In Sec. 4, we discuss nonlinearity in scaling relations for fundamental-mode RR Lyrae stars, and we compare our observational data with stellar pulsation models. In Sec. 6, we outline an example of the determination of a systemic velocity using our derived templates and scaling relations together with their testing. Last, Sec. 7 provides a summary of our results.

2 Line-of-sight velocity templates for APOGEE spectra

This section describes our approach to deriving the amplitude scaling relations and line-of-sight velocity templates for RR Lyrae stars. We focus on single-mode RR Lyrae stars – those pulsating in the fundamental mode, RRab, as well as those pulsating in the first overtone, RRc – by using APOGEE near-infrared spectra for RR Lyrae stars in the Solar neighborhood. In deriving the aforementioned tools, we have drawn inspiration from the studies by Sesar (2012) and Braga et al. (2021) that undertook a similar process but in the visual part of the spectra.

APOGEE provides various spectroscopic products. In particular, we were interested in the individual visits (apVisit), which are a combination of approximately eight individual observations (observed back-to-back), each with 500 s exposures. These individual observations comprise the foundation for creating the line-of-sight velocity templates and constructing the amplitude scaling relations for RR Lyrae pulsators observed by the APOGEE survey.

2.1 Obtaining the line-of-sight velocities

The APOGEE survey intends to observe approximately 10000 RR Lyrae stars both on the southern and northern sky (Bowen & Vaughan 1973; Gunn et al. 2006; Holtzman et al. 2010; Zasowski et al. 2017). In the northern hemisphere, a small sample of 30 nearby (H < 10 mag) RR Lyrae variables were observed by APOGEE to obtain, at least for some of them, sufficient coverage of the entire pulsation cycle that would then serve in creating the line-of-sight velocity template (see Table A.1 for their full list). These RR Lyrae stars have been marked with the rrlyr label in the APOGEE data products. The APOGEE survey serendipitously also observed other nearby RR Lyrae stars (listed as RR Lyrae stars in the International Variable Star Index, VSX, database¹, Watson et al. 2006). We decided to utilize some of these RR Lyrae stars (we used only those with at least five visits and S/N > 50) to create the line-of-sight velocity templates for both RR Lyrae subclasses and the relation between photometric amplitudes and line-of-sight velocity amplitudes (a list of these RR Lyrae stars used in the creation of the line-of-sight-velocity templates is given in Table A.2).

As an initial step, we acquired the individual exposures from the SDSS ftp data access². The APOGEE spectra are provided in vacuum wavelengths (Majewski et al. 2017); therefore, we transformed them into the air wavelengths using a relation from the SDSS website³. For each exposure, we determined the barycentric observation time T_bjd. To determine reliable barycentric line-of-sight velocities (υ_los), we used the iSpec package Blanco-Cuaresma et al. (2014); Blanco-Cuaresma (2019) to synthesize a grid of synthetic spectra that cover common physical properties of RR Lyrae stars (based on stellar parameters found in the literature, For et al. 2011; Sneden et al. 2017; Preston et al. 2019). We used the following set of physical properties to create synthetic spectra: [Fe/H] = (-2.5,-2.0,-1.5,-1.0,-0.5,0.0) dex; T_eff = (6000., 6500., 7000., 7500.) K; log g = (2.0, 2.5, 3.0) dex; and a microturbulence velocity of ξ_turb = 3.5 km s⁻¹.

The spectra were synthesized using the radiative transfer code MOOG (February 2017 version, Sneden 1973) implemented in iSpec, the ATLAS9 model atmosphere (Castelli & Kurucz 2003), the solar reference scale from Asplund et al. (2009), and the line list from the VALD database⁴. The spectral range was tailored to APOGEE spectra (from 15000 Å to 17 000 Å). Thus, we covered the entire APOGEE spectral range and did not focus on small spectral features, such as specific spectral lines (e.g., studies by Sesar 2012; Braga et al. 2021). This decision was motivated by the intensities of the hydrogen lines, members of the Paschen and Brackett series that dominate this region of spectra in RR Lyrae stars. These hydrogen lines are the only ones seen in low signal-to-noise APOGEE spectra of RR Lyrae variables toward the Galactic bulge.

To determine the line-of-sight velocity (υ_los), we cross-correlated each spectrum with synthesized templates, selecting the one that minimized the error in υ_los. This template was then used in a Monte Carlo simulation (200 iterations for each measurement) to address the uncertainties in the observed spectra. We adjusted flux values according to their normally distributed errors. Finally, the bootstrap method was applied to this data to calculate the average and standard deviation of υ_los and its variation, σ_υlos. To remove erroneous measurements, we applied a condition on the υ_los and on its uncertainty: $$\left|v_{\mathrm{l}\mathrm{o}\mathrm{s}}/{\sigma }_{v_{\mathrm{l}\mathrm{o}\mathrm{s}}}\right|>5\text{ or }{\sigma }_{v_{\mathrm{l}\mathrm{o}\mathrm{s}}}<10\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}\text{. }$$(1)

2.2 Estimating the line-of-sight velocity amplitudes

Once we acquired υ_los and T_hjd for the nearby RR Lyrae stars, we focused on obtaining photometric information to constrain the basic pulsation properties of our RR Lyrae template dataset. We used photometric data from the Gaia RR Lyrae catalog (Clementini et al. 2023), which provides photometry in G, G_rp, and G_bp as well as pulsation properties (pulsation period P, times of brightness maxima E, and amplitudes of the light and velocity changes Amp_G, Amp_GRP, Amp_GBP, and Amp_RV).

Once we obtained the Gaia ephemerides (P and E) for the stars in our dataset, we phased the observed υ_los using the element-wise remainder of division of the product, and we assessed the quality of the line-of-sight velocity curves for our RR Lyrae template dataset. We selected only those stars for which both the light descending side and light ascending side of the line-of-sight velocity curve were covered based on visual inspection of the data. In the end, we had 13 RR Lyrae stars (eight RRab and five RRc stars) that were suitable for estimating line-of-sight velocity amplitudes (Amp_los) and subsequently creating the line-of-sight velocity templates. We refer to these stars as our APOGEE template dataset, and they are marked with an asterisk in Tables A.1 and A.2. Our dataset of 13 template stars is somewhat similar in number to that of Sesar (2012), who used six RRab stars to construct radial velocity templates, and Braga et al. (2021), who used 31 RRab but five RRc stars to obtain radial velocity templates.

In the next step, we estimated the line-of-sight velocity amplitudes. Despite using individual exposures, our line-of-sight velocity curves do not cover the entire phase curve without small gaps, and the curves can have a low number of points (e.g., 27 reliable measurements for the prototype of the RR Lyrae variable class). Therefore, a direct Fourier fitting (as used in Braga et al. 2021) could lead to problems in obtaining a reliable Amp_los. Thus, similar to Dékány et al. (2021), we turned to the Gaussian process regression (GPR; Rasmussen & Williams 2005) as a means to describe phased curves of RR Lyrae stars (but in our case, we work with the line-of-sight velocity curves). We used phase-folded line-of-sight velocities in conjunction with the Gaussian process (GP) regressor implemented in the george module (Ambikasaran et al. 2015)⁵. For the GPR, we used the following kernels: $$\begin{array}{rl}\text{ kernel }=& \text{ ExpSquaredKerne1(hyperparameters) }\\ & \times \text{ Matern32Kernel(hyperparameters). }\end{array}$$

The hyperparameters for each kernel were optimized through the marginalized log-likelihood function⁶.

In Fig. 1, we depict a phased line-of-sight velocity curve based on the determined υ_los from individual exposures, line-of-sight velocities from the co-added visits, and a fit of the GP regressor for estimating the amplitude of the line-of-sight velocity changes. The GP model for each individual star does introduce some curves along the rising branch, likely due to the gaps in the observations along the full pulsation cycle, but much of this is smoothed out when the full sample of stars is combined in the generation of the final RRab and RRc template. In Table 1, we list the basic pulsation properties and the line-of-sight velocity amplitudes⁷ for stars used in the creation of the amplitude scaling relation and the υ_los curve template.

Fig. 1 Example of a single exposure spectrum (top panel) and line-of-sight velocities (bottom panel) used in this work. The top panel depicts a spectrum of TU-UMa taken at HJD = 2 457 484.880667 day. The bottom panel shows the line-of-sight velocities derived in this work (black dots) and the velocities extracted from APOGEE data products (red squares). The gray dots represent estimated velocities that did not fulfill the criteria in Eq. (1).

2.3 Scaling relation between Amp_los and Amp_G based on APOGEE spectra

Using the photometric amplitudes obtained from the Gaia catalog, Amp_G, and the acquired Amp_los from the APOGEE, we derived the scaling relation between both parameters. Previous studies, such as Sesar (2012) and Braga et al. (2021), used a linear scaling relation between the photometric and velocity products. We derived separate amplitude relations for both RR Lyrae subclasses. In the case of RRc-type stars, we also used a linear relation, but for RRab, we decided to use a second-degree polynomial relation. This is supported by the observed trend for RRab stars in APOGEE data and subsequently confirmed in Sect. 3 for Gaia RVS spectra.

We implemented a fitting procedure aimed at optimizing the following relations: $${\mathrm{A}\mathrm{m}\mathrm{p}}_{\mathrm{l}\mathrm{o}\mathrm{s}}^{\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}=x_0\cdot {\mathrm{A}\mathrm{m}\mathrm{p}}_G^2+x_1\cdot {\mathrm{A}\mathrm{m}\mathrm{p}}_G$$(2) $${\mathrm{A}\mathrm{m}\mathrm{p}}_{\mathrm{l}\mathrm{o}\mathrm{s}}^{\mathrm{R}\mathrm{R}\mathrm{c}}=x_0\cdot {\mathrm{A}\mathrm{m}\mathrm{p}}_G,$$(3)

where x₀ and x₁ are parameters of the fit (we omitted the intercept to preserve the physical meaning of the relation). In the search for optimal parameters and their covariances, we utilized the emcee module (Foreman-Mackey et al. 2013). For each RR Lyrae subclass, we ran the Markov chain Monte Carlo simulation with 200 walkers for 30000 samples. For deriving the best-fit values, we thinned the sample by τ = 20 and marked the initial 10000 samples as burn-in.

Figure 2 shows the scaling relations for both single-mode RR Lyrae stars. The scaling relations for the RRab and RRc pulsators, together with their uncertainties, can be expressed with the following equations: $$\text{ For RRab stars: }{\mathrm{A}\mathrm{m}\mathrm{p}}_{\mathrm{l}\mathrm{o}\mathrm{s}}=-64\cdot {\mathrm{A}\mathrm{m}\mathrm{p}}_G^2+127\cdot {\mathrm{A}\mathrm{m}\mathrm{p}}_G$$(4) $$\text{ For RRc stars: }{\mathrm{A}\mathrm{m}\mathrm{p}}_{\mathrm{l}\mathrm{o}\mathrm{s}}=62(2)\cdot {\mathrm{A}\mathrm{m}\mathrm{p}}_G\cdot$$(5)

For RRab variables, we also include a covariance matrix for the uncertainty estimation: $${\mathrm{Cov}}^{\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}=\left[\begin{array}{cc}120& -98\\ -98& 81\end{array}\right].$$(6)

2.4 Line-of-sight velocity templates for APOGEE spectra

When creating the line-of-sight velocity templates for RRab and RRc stars, one needs first to normalize the line-of-sight velocity curves of the stars. We utilized the phased υ_los curves and their associated GPR models. We used GPR models for the individual calibration stars to normalize the phase-folded υ_los curves (between −0.5 and 0.5). The individual normalized phased curves were then co-added into two normalized curves, one for each RR Lyrae subclass. Unlike previous studies, we determined the systemic velocity directly from the shift of the line-of-sight velocity template and not from a specific phase point (e.g., 0.38 or 0.5; see Liu 1991; Sesar 2012).

These template curves, as we refer to them, were then modeled using Fourier decomposition with the following relation: $$f(\phi )=a_0+{\displaystyle \sum _{k=1}^n}{a}_k\cdot \mathrm{sin\hspace{0.17em}}[2\pi k\varphi ]+a_{k+1}\cdot \mathrm{cos\hspace{0.17em}}[2\pi (k+1)\varphi ],$$(7)

where ϕ represents the time element (in this case, phase) and a_k represents the amplitudes of the individual Fourier components. For RRab and RRc stars, we selected the fifth (n = 5) and third degree (n = 3) of the Fourier series, respectively. In Fig. 3, we depict the phase-folded normalized line-of-sight velocity curves with associated Fourier templates. In Table 2, we list the Fourier coefficients for both models of RR Lyrae subclasses. The a₀ coefficients are zero due to the prior normalization of the line-of-sight velocity curves.

In the bottom panels of Fig. 3, we observed a noticeable scatter in the template curves, particularly between phases of rapid contraction (0.8 to 1.0). We included this scatter in υ_sys determination (see example in Sect. 6.1) as a root-mean-square scatter in the template, which is similar to the approach described in Sesar (2012). This can be represented as a spline function describing the scatter along the pulsation phase (see Table A.3).

Fig. 2 Scaling relations between the photometric and line-of-sight velocity amplitudes for the fundamental-mode RR Lyrae stars (left-hand panel) and first-overtone RR Lyrae variables (right-hand plot). Both RR Lyrae subclasses are depicted using blue points. The best-fit parameters are shown with a solid red line, while gray lines represent their variations.

Fig. 3 Line-of-sight velocity templates for the RRab (left-hand panel) and RRc (right-hand plots) pulsators. The black points in the top panels represent individual measurements of υlos, and the solid red lines depict the Fourier model for each normalized line-of-sight velocity curve. The bottom panels show the difference, Δ, between the model (Fourier fit) and the phased υlos data for both RR Lyrae subclasses.

3 Gaia calcium triplet spectra

The Gaia RR Lyrae catalog provides all the necessary data to develop scaling relations and line-of-sight velocity templates for the Ca T region. The abundance of data allowed for the creation of a tailored method for Ca T spectra observed by the BRAVA-RR survey for RR Lyrae pulsators toward the Galactic bulge (Kunder et al. 2016, 2020).

Using photometric amplitudes from the Gaia survey for calibration of the scaling relations comes with a trade-off. Gaia provides exceptional photometry for hundreds of thousands of RR Lyrae stars taken simultaneously with the Gaia RVS spectra. However, approximately 20 percent of the RR Lyrae stars in the Galactic bulge (Clementini et al. 2023), which are covered by the APOGEE and BRAVA-RR surveys, do not have photometric amplitudes in Gaia passbands. The main source of photometric amplitudes for these stars comes from the fourth data release of the Optical Gravitational Lensing Experiment (OGLE-IV; Soszyński et al. 2014, 2019; Udalski et al. 2015). OGLE-IV provides a detailed sample of RR Lyrae stars in the Galactic bulge, including abundant photometry in the V and I passbands.

Thus, we decided to create conversion relations between the V, I, and G passbands to enable transitions to and from the Gaia G-band in cases where amplitudes in other passbands are available. To this end, we crossmatched the Gaia RR Lyrae catalog with the OGLE-IV catalog using equatorial coordinates (within a 1 arcsecond radius). This crossmatch provided us with photometry for Galactic bulge RR Lyrae stars in the V and I passbands. We utilized our fitting routine (similar to the one used in Eq. (7)) to decompose the OGLE-IV light curves using the Fourier series. The following Fourier light curve decomposition was optimized: $$m(t)=m_V+{\displaystyle \sum _{k=1}^n}{A}_k\cdot \cos (2\pi k\vartheta +{\phi }_k).$$(8)

In Eq. (8), m_v represents the mean apparent magnitude, and A_k and φ_k stand for amplitudes and phases, respectively. The term n denotes the degree of the fit that we adapted for each light curve using the same approach as in Prudil et al. (2019). Here, ϑ represents the phase function, which is defined as: $$\vartheta =(\mathrm{H}\mathrm{J}\mathrm{D}-E)/P,$$(9)

where HJD is the Heliocentric Julian Date of the observation and E stands for the time of brightness maximum. The time of brightness maxima was used as E to set a common light curve reference point.

To create conversion relations, we followed the same procedure as in Sect. 2.3 for the first-overtone pulsators, selecting only a linear fit without an intercept. To exclude potentially blended stars, we applied the condition ipd_frac_multi_peak < 5. We fitted amplitude comparisons for RRab- and RRc-type stars separately, and the derived relations are listed below: $$\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}:{\mathrm{A}\mathrm{m}\mathrm{p}}_I=0.80(0.01)\cdot {\mathrm{A}\mathrm{m}\mathrm{p}}_G$$(10) $$\text{ RRab: }{\mathrm{A}\mathrm{m}\mathrm{p}}_V=1.32(0.01)\cdot {\mathrm{A}\mathrm{m}\mathrm{p}}_G$$(11) $$\text{ RRc: }{\mathrm{A}\mathrm{m}\mathrm{p}}_I=0.74(0.02)\cdot {\mathrm{A}\mathrm{m}\mathrm{p}}_G$$(12) $$\text{ RRc: }{\mathrm{A}\mathrm{m}\mathrm{p}}_V=1.33(0.04)\cdot {\mathrm{A}\mathrm{m}\mathrm{p}}_G.$$(13)

The derived conversion relations were also used in the cases where we needed to compare scaling relations optimized on V passband.

3.1 Scaling relation between Amp_G and Amp_los derived from Gaia spectra

Our method for obtaining scaling relations and line-of-sight velocity templates closely mirrors the process we utilized for the APOGEE dataset. We extracted pulsation properties from the Gaia catalog (pulsation period and epoch of the maximum of the light curve in the G-band) for all single-mode RR Lyrae stars with publicly available individual υ_los and their associated uncertainties (in total, 1086 variables). We implemented the same criteria for the line-of-sight velocities as for APOGEE (see Eq. (1)). Moreover, we selected RR Lyrae stars that had at least 20 Gaia υ_los measurements (after applying criteria on the line-of-sight velocities), and we only selected RR Lyrae stars for which the measured Amp_los fulfilled the following criteria: $$\left|{\mathrm{A}\mathrm{m}\mathrm{p}}_{\text{los }}/{\sigma }_{{\mathrm{A}\mathrm{m}\mathrm{p}}_{\text{los }}}\right|>4\text{ or }{\sigma }_{{\mathrm{A}\mathrm{m}\mathrm{p}}_{\text{los }}}<10\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}\text{. }$$(14)

Our last criterion was to only use stars above Gaia flag ipd_frac_multi_peak < 5. The last criterion was followed in order to reduce the possible blending of stars, as it would decrease pulsation and velocity amplitudes. In total, for RRab-and RRc-type stars, we had 245 and 88 Amp_los to derive scaling relations, respectively.

As in the APOGEE case, we used a GP regressor to calculate Amp_los. While modeling the amplitude scaling relation, we noticed a non-negligible scatter in the data, and we included intrinsic scatter in the fit (for both RRab and RRc variables, ε_RRab and ε_rrc). To further explore the nature of the scatter, we crossmatched our fundamental mode dataset with a catalog of known RRab stars exhibiting the Blazhko effect (Skarka 2013). We found that approximately 9 percent of our RRab sample shows signs of modulation, and thus, some of the scatter in the RRab scaling relation may be due to the Blazhko effect. Including or removing the identified Blazhko stars did not have an effect on the scatter in the scaling relation to RRab stars. This was partially caused by the fact that we do not have complete information on the identification of Blazhko stars in the Solar neighborhood.

In Fig. 4, we depict the scaling relation for both RR Lyrae subclasses. It is typically assumed that there is a linear relation between Amp_G and Amp_los (for both RR Lyrae subtypes), but it is apparent that the amplitudes do not follow a linear relation (similar to the case we saw in Sect. 2.3). Instead, there is a break at Amp_G ≈ 0.4 mag, where stars with an Amp_G smaller than 0.4 have a steeper scaling relation than those with Amp_G > 0.4 mag. For comparison, we overplotted Gaia data and our derived relation for RRab stars with linear relations found in the literature (Liu 1991; Sesar 2012; Braga et al. 2021, Fe lines;). For interested readers, in the appendix, we enclosed a Table A.4 with 34 RRab variables that are located at the nonlinear part of the scaling relation with some of their basic properties. We observed that for an RR Lyrae variable with Amp_G = 1.0 mag, the scaling relations from the literature (based on Fe lines and converted from Amp_V to Amp_G; Liu 1991; Sesar 2012; Braga et al. 2021) together with our relation predict Amp_los ≈ 68 km s⁻¹. On the opposite side of the amplitude distribution, for low amplitudes (e.g., Amp_G = 0.25 mag), our relation predicts amplitudes lower than Amp_los ≈ 10 km s⁻¹ (for Liu 1991, Sesar 2012 and Braga et al. 2021). One consequence of a steeper relation between Amp_G and Amp_los at a smaller Amp_G is that the intercept of the scaling relation naturally approaches zero, which is what is physically expected – a photometric amplitude of zero should also correspond to no Amp_los. On the other hand, for the first-overtone pulsators, we did not see a rapid decline of Amp_los toward smaller photometric amplitudes. The range of both G-amplitudes and Amp_los is considerably smaller than for the RRab pulsators, which may be a factor contributing to more constant scaling relations. Also, the scaling relation of RRc pulsators naturally passes through an intercept of zero without the need for a change in slope. More RRc-type variables at the low-amplitude end would help solidify the use of a linear relation instead of a higher degree polynomial.

The scaling relations between Amp_los and Amp_G, their associated covariance matrix, and uncertainty are described in the following equations: $$\begin{array}{r}\text{ For }\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}\star :{\mathrm{A}\mathrm{m}\mathrm{p}}_{\mathrm{l}\mathrm{o}\mathrm{s}}=39\cdot {\mathrm{A}\mathrm{m}\mathrm{p}}_G^3-133\cdot {\mathrm{A}\mathrm{m}\mathrm{p}}_G^2\\ +159\cdot {\mathrm{A}\mathrm{m}\mathrm{p}}_G.\end{array}$$(15) $$\text{ For RRc }\star :{\mathrm{A}\mathrm{m}\mathrm{p}}_{\mathrm{l}\mathrm{o}\mathrm{s}}=62(2)\cdot {\mathrm{A}\mathrm{m}\mathrm{p}}_G,{\epsilon }_{\mathrm{R}\mathrm{R}\mathrm{c}}=4\text{. }$$(16) $${\mathrm{Cov}}^{\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}=\left[\begin{array}{ccc}172& -271& 97\\ -271& 438& -163\\ 97& -163& 64\end{array}\right]{\epsilon }_{\text{RRab }}=5.$$(17)

Fig. 4 Same as Fig. 2 but for line-of-sight velocity amplitudes determined based on Gaia calcium triplet spectra. In addition, the green solid, dashed, and dotted lines represent linear relations from Braga et al. (2021), Sesar (2012), and Liu (1991, using p-factor equal to 1.35), respectively. We note that the literature relations were scaled using Eqs. (11) and (13).

3.2 Line-of-sight velocity templates for Gaia spectra

To generate line-of-sight velocity templates from measured υ_los, we proceeded in the same way as for the APOGEE spectra (see Sect. 2.4). In addition to the procedure for the APOGEE spectra, thanks to the large number of measurements, we also included a criterion on the maximum acceptable root-mean-square error (RMSE; measured using available data and the GP regressor for each velocity curve) for a given line-of-sight velocity curve in the following form: RMSE < 4. We again normalized the phased line-of-sight velocity curves and decomposed them with Eq. (7) using the fifth and third degree for the fundamental and first-overtone variables, respectively. Following Braga et al. (2021) and taking advantage of the large dataset, we divided our RRab sample into three groups based on their pulsation periods: RRab-1 (pulsation periods below 0.55 day), RRab-2 (pulsation periods above 0.55 day and below 0.7 day), and RRab-3 (pulsation periods above 0.7 day). This division of RRab variables in period bins is based on the arguments presented in Braga et al. (2019, see their Sect. 3.1 for details). In principle, the photometric amplitudes are not linearly dependent on pulsation periods, which are tightly connected to stellar parameters (e.g., mass, luminosity, metallicity).

In addition to the criteria mentioned in the previous subsection, we selected only the highest quality line-of-sight velocity curves. These include those on the bright end of the mean magnitude distribution (V < 12.0 mag) and those with at least 25 measurements of υ_los. We made an exception in the case of the RRab-3 category, relaxing the criterion for the number of measurements to at least 20 measurements of υ_los. Final templates for all four categories can be found in Fig. 5, and their Fourier coefficients are listed in Table 3. As in the case of APOGEE templates, we estimated the RMSE along the templates and list it in Table A.5. The RR Lyrae variables used for template creation are listed in Table A.6.

4 Nonlinearity of the scaling relations

One unanticipated finding illustrated in Fig. 4 is that the scaling between photometric and line-of-sight velocity amplitudes for fundamental mode pulsators is not linear. Instead, a break at ≈0.4 mag (≈45 km s⁻¹) is apparent, which will primarily affect the shorter period RRab stars. In this section, we attempt to understand the underlying processes governing this behavior through the use of nonlinear pulsation models.

To compare the nonlinear trend in the RRab scaling relation, we used photometric data from the All-Sky Automated Survey for Supernovae (ASAS-SN; Shappee et al. 2014; Jayasinghe et al. 2018). The ASAS-SN survey uses 24 small telescopes (14 cm in diameter) located all around the globe, continuously scanning the night sky. Photometry in the V-band and g-band are provided for many RR Lyrae stars, with typical uncertainties of 0.02 mag for stars brighter than 12.0 mag in V-mag. We used the ASAS-SN V-band photometry to establish the basic pulsation properties (pulsation periods, times of brightness maxima, and amplitudes of the light changes) for our RR Lyrae template dataset. We proceeded in a manner similar to our analysis in Sect. 3 by optimizing the Fourier light curve decomposition through Eq. (8). For this analysis, we used only RRab stars that were utilized in Sect. 3.1. In addition to the criteria for scaling relations, we used only stars above b > 10deg. This criterion was selected to reduce the possible blending of stars, which would decrease pulsation amplitudes, particularly since we were using the ASAS-SN photometry, which has large pixel sizes (8 arcsecs).

The pulsation models were computed with the Radial Stellar Pulsation (RSP) tool implemented as part of the Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023). We used the MESA r21.12.1 public release. The RSP tool computes Lagrangian nonlinear convective pulsation models with time-dependent convection-pulsation coupling given by the Kuhfuss (1986) model. In general, numerical implementation follows the work of Smolec & Moskalik (2008a). First, RSP builds a static equilibrium model. Then, linear stability analysis is conducted, which yields linear pulsation periods and growth rates of the radial modes. Finally, after applying initial perturbation following the line-of-sight velocity eigenvector of a given mode, the model is integrated in time until it approaches a finite amplitude, single-periodic limit cycle pulsations.

We built a chemically homogeneous envelope with 180 cells for all model computations. Sixty surface cells have an equal mass down to the anchor zone where the temperature is fixed to 11 × 10³ K. Then, the mass of the cells increases geometrically inward down to a temperature of 2 × 10⁶ K. The envelope is nonrotating, with constant luminosity at the radiative, rigid bottom boundary. The OPAL opacities were used (Iglesias & Rogers 1996), and at lower temperatures, they were supplemented with the Ferguson et al. (2005) opacity data. The distribution of metals followed the solar abundance distribution as given by Asplund et al. (2009). Our basic model sequences adopted convective parameters of set A, as given in Table 4 of Paxton et al. (2019).

We were interested in computing fundamental-mode models covering a full spectrum of pulsation amplitudes, starting with those exhibiting the smallest amplitudes, which are expected to be close to the edges of the fundamental-mode instability strip (IS). Therefore, we considered model sequences with a constant mass and luminosity and an increasing effective temperature that start close to the red edge of the fundamental-mode IS. The red edge for the fundamental mode is well beyond the IS for the first overtone, which is located at hotter temperatures. Consequently, the models developed full-amplitude fundamental-mode pulsation, which is the regime of interest for this case. If we constructed our sequence of models near the blue edge of the fundamental-mode IS, the models would also be inside the first-overtone IS. Consequently, near the blue boundary for the fundamental mode IS, the model would switch to pulsations in the first overtone instead of developing pulsations in the fundamental mode. Pulsations in the fundamental mode would be possible for lower temperatures, further from the blue edge of the fundamental mode IS, but then the amplitudes would already be significant. (For more information on the mode selection, see the works of, e.g., Szabó et al. 2004; Smolec & Moskalik 2008b; Paxton et al. 2019).

Altogether, we considered eight model sequences: four with [Fe/H] = −1.0 dex (Z = 0.0014, X = 0.75) and four with [Fe/H] = −2.0 dex (Z = 0.00014, X = 0.75). For each metallicity, the adopted mass is either 0.65 M_&ocir;, or 0.75 M_&ocir;, and the adopted luminosity is either 45 L_&ocir;, or 50 L_&ocir;. We used Fourier amplitudes (A₁) instead of peak-to-peak amplitudes to compare with observational data. While pulsation models reproduce the low-order Fourier parameters of the light curves well, they fail when it comes to sharp features in short wavelengths (see, e.g., Marconi et al. 2015; Das et al. 2018; Paxton et al. 2019). In fundamental mode, RR Lyrae models quite often have sharp spikes that develop close to the maximum light and affect the peak-to-peak amplitude, but these sharp spikes are not observed in the light curves. In Fig. 6, we compare the model amplitudes with observations. For the line-of-sight velocity amplitude, we considered photosphere pulsation velocities converted to observed line-of-sight velocities using a fixed projector factor of 1.35 (Kovács 2003). The overall agreement between models and observations is very good; in particular, the correlation between the velocity and light curve amplitudes is well reproduced. Theoretically, at low amplitudes, the relation between the photometric and line-of-sight velocity is linear. It then breaks at ≈0.2 mag but soon becomes linear again, with a shallower slope. The spread of the line-of-sight velocity amplitudes at a given light amplitude may be explained as being due to different physical parameters of the stars (e.g., different masses, luminosities, metallicities). While the overall agreement is satisfactory, we note that the model relation seems to be shifted toward higher photometric amplitudes in a low-amplitude regime; conversely, the predicted line-of-sight velocity amplitude is too small at a given light amplitude. Different projection factors cannot rectify this shift, especially at the low-amplitude regime.

Physical parameters (mass, luminosity, metallicity) are not the only ones to affect light and line-of-sight velocity amplitudes, or more broadly their shape. Parameters of the convection model are also well known to affect the light and line-of-sight velocity curves. To check their effect, for one set of physical parameters (0.65 M_&ocir;, 45 L_&ocir;, [Fe/H] = −2.0 dex), we considered three more model sequences with convective parameters corresponding to those of sets B, C, and D of Paxton et al. (2019, their Table 4) in addition to set A, for which previously described model sequences were computed. Sets B and D include the effects of radiative losses, while in sets C and D, the effects of turbulent pressure and kinetic turbulent energy flux are included. The results are illustrated in Fig. 7. We observed that in the low-amplitude regime, the relation between amplitudes is similar for sets A and B and qualitatively different than for sets C and D. For the latter two, the agreement with observations is significantly better, that is, for a given light amplitude, we observed higher line-of-sight velocity amplitudes. The observations in the low-amplitude range are well matched with sets C and D. These two sets include the effects of turbulent pressure and turbulent kinetic energy flux. It seems crucial, however, that the models of low amplitude in sets C and D are significantly cooler than the models of sets A and B. This is because the red edge of the IS, which is where we started the model sequences, shifts toward lower effective temperatures by about 150 K for sets C and D. For a higher amplitude range, we observed that the models including radiative losses (sets B and D) lead to lower line-of-sight velocity amplitudes at a given light amplitude than sets that neglect radiative losses. We note that there are no strong theoretical constraints for the parameters of the convective model. While some sets are commonly used in the literature, a comprehensive calibration employing numerous observational constraints (e.g., full shapes of the line-of-sight velocity and multi-band light curves) has been performed only sporadically (e.g., Di Fabrizio et al. 2002). Our results illustrate that simultaneous matching of the light and line-of-sight velocity curves is a possible avenue for calibrating convection models.

Fig. 5 Same as Fig. 3 but for the line-of-sight velocity curves determined based on Gaia calcium triplet spectra.

Fig. 6 Relation between the photometric and line-of-sight velocity Fourier amplitudes for RRab stars (gray dots) along with nonlinear pulsation models of [Fe/H] = −1.0 dex (top panel) and [Fe/H] = −2.0 dex (bottom panel) and different masses and luminosities, as indicated in the legend.

Fig. 7 Relation between the photometric and line-of-sight velocity Fourier amplitudes for RRab stars (gray dots) along with four sequences of nonlinear pulsation models. We used the following fixed physical parameters (M = 0.65 M⊚, L = 45 L⊚, [Fe/H] = −2.0 dex, with varying effective temperature) and assumed different sets of convective parameters.

5 Cross-survey analysis: APOGEE-Gaia velocity templates

In this section, we explore the feasibility of exchanging derived line-of-sight velocity templates and scaling relations between the APOGEE (15 000–-17 000 Å) and Gaia (8450–8720 Å) surveys. This approach is motivated by the limited number of calibration RR Lyrae stars in the APOGEE survey, specifically both RRab and RRc pulsators. This limitation results in increased scatter across the template. We also evaluate how previously derived templates, particularly from Braga et al. (2021) and Liu (1991), align with our line-of-sight velocity templates, with a focus on the Braga et al. (2021) template for RRc pulsators. The interchangeability of the line-of-sight velocity curves has been explored in Clementini et al. (2023), and the authors found that the line-of-sight velocity measurements from the Gaia RVS spectra are in good agreement with the line-of-sight velocities in the literature despite the different spectral ranges.

In Fig. 8, we present a comparison between the newly derived line-of-sight velocity templates based on APOGEE and Gaia spectra. In addition, we show a comparison with the previously derived templates in Liu (1991) and Braga et al. (2021) for Fe lines. The top panel shows that the line-of-sight velocity templates from APOGEE and Gaia are largely interchangeable, particularly when considering the scatter associated with each template. This indicates that templates developed using Gaia spectra can effectively be employed for determining systemic velocities from APOGEE data. Furthermore, the bottom panel of Fig. 8 illustrates that the line-of-sight velocity templates derived from iron lines in Braga et al. (2021) are also applicable to APOGEE and Gaia data. This compatibility is particularly notable for the RRab templates, where all three template sources (APOGEE, Gaia, and Braga et al. 2021) exhibit remarkably similar patterns. The only notable deviation occurs with the RRab template from Liu (1991), which shows an offset that is also in phase; adjusting it by −0.05 would align it closely with the APOGEE template. This comparative analysis of the templates implies that various sources for the line-of-sight velocity templates are effectively interchangeable within the errors.

We focus next on the substitutability of the scaling relations between surveys and the spectral lines used in the determination of line-of-sight velocities. Previously, in Fig. 4, we demonstrated that the scaling relations for the RRab variables from Liu (1991) and Braga et al. (2021) accurately represent the higher amplitude regime. For a specific RRab star with a given photometric amplitude of Amp_G = 0.8 mag (Amp_v = 1.06 mag, using Eq. (11)), the line-of-sight velocity amplitudes calculated using the relations from Liu (1991), Sesar (2012), Braga et al. (2021), and Eq. (15) fall within the range of 61 to 62 km s⁻¹. Using the scaling relation for RRab stars with APOGEE data (see Eq. (4)) yields a higher Amp_los = 61 km s⁻¹.

A similar conclusion can be drawn for RRc stars. We compared the amplitude scaling relations derived for the RRc stars based on Gaia, APOGEE, and the Braga et al. (2021) relation (based on iron lines of first-overtone pulsators). This comparison revealed that a typical RRc with an Amp_G = 0.35 mag (Amp_v = 0.47 mag, using Eq. (13)) has a predicted Amp_los equal to 22, 22, and 28 km s⁻¹ for relations from Braga et al. (2021), Eqs. (4) and (16), respectively. Our predictions suggest a slightly lower Amp_los in comparison with the prediction by Braga et al. (2021). There could be several reasons for this, especially the differences in covered amplitude space. Our Gaia dataset covers both a higher and lower Amp_G range than Braga et al. (2021). In the APOGEE dataset, we do not have as many data points as were used in Braga et al. (2021) and in our Gaia dataset.

The comparisons outlined above demonstrate the possibility of substituting templates and scaling relations across surveys. For example, one can use line-of-sight velocity templates and scaling relations based on Gaia spectra to derive systemic velocities for line-of-sight velocities estimated from APOGEE spectra. Despite differences in wavelength coverage and spectral lines between the surveys, the estimated systemic velocities using either method are comparable, as long as we are not in a low-amplitude regime for RRab stars. This yields an advantage, especially for RRab stars that fall into the low-amplitude space, where we can use the well-described scaling relation based on the Gaia dataset across different line-of-sight velocity sources. One can even use the Gaia line-of-sight velocity templates for velocities calculated based on iron lines. In principle, all the variations should lead to equivalent predicted Amp_los and systemic velocities within the uncertainties.

6 An example of determining v_sys and testing the newly derived templates

In what follows, we describe a method for estimating systemic velocities for a selected RR Lyrae variable using the derived templates and scaling relations presented above. We assess the resulting systemic velocity and apply our method to the large portion of the Gaia RR Lyrae catalog with υ_los⁸.

Fig. 8 Comparison between line-of-sight velocity templates derived in this work (top panels, based on APOGEE and Gaia data) and those derived in previous studies (bottom panels, together with APOGEE line-of-sight velocity template). For each line-of-sight velocity template derived in this work, we also included a shaded region (in top and bottom panels) depicting the measured scatter along the template (see Tables A.3 and A.5).

6.1 An example of determining υ_sys using APOGEE data

Most of the RR Lyrae stars observed by the APOGEE survey have well-known pulsation properties (pulsation period, amplitude, and time of brightness maxima) but only a few spectroscopic observations (1–3), so the coverage of the pulsation cycle is incomplete. In cases similar to this, the line-of-sight velocity template and scaling relation need to be used to derive the υ_sys and σ_υsys. For this exercise, we chose BH Aur, which belongs to the RRab subclass, with P = 0.4560898 day, Amp_G = 0.81 ± 0.03 mag, and E = 2 456 907.7783865286 day (based on Gaia photometry). It has only four APOGEE subvisit observations and a well-determined υ_sys = 51.4 ± 1.5 km s⁻¹ in the Gaia survey⁹ due to a large number of spectroscopic observations (28 in total, see Gaia RR Lyrae catalog Clementini et al. 2023).

In the determination of υ_sys and σ_υsys, we first converted the photometric amplitude in the G-band to the Amp_los using the relation in Eq. (4) and estimated the uncertainty on the Amp_los, ${\sigma }_{{\mathrm{A}\mathrm{m}\mathrm{p}}_{\mathrm{los}}}$. This allowed us to assemble a dataset D for BH Aur that consisted of its pulsation properties, measured υ_los (including their associated uncertainties), and the time of observation represented by the Barycentric Julian Date (based on APOGEE subvisit spectra): $$\mathbf{D}=\{P,E,{\mathrm{A}\mathrm{m}\mathrm{p}}_{\mathrm{l}\mathrm{o}\mathrm{s}},{\sigma }_{{\mathrm{A}\mathrm{m}\mathrm{p}}_{\mathrm{l}\mathrm{o}\mathrm{s}}},v_{\mathrm{l}\mathrm{o}\mathrm{s}},{\sigma }_{v_{\mathrm{l}\mathrm{o}\mathrm{s}}},T_{\mathrm{B}\mathrm{J}\mathrm{D}}^{\mathrm{A}\mathrm{P}\mathrm{O}}\}.$$(18)

To find υ_sys and its uncertainty given the dataset D, we optimized a subsequent relation for likelihood p(D | υ_sys) defined as $$p(\mathbf{D}\mid v_{\text{sys }})=\mathcal{N}(v_{\text{los }}^{\text{model }}\mid v_{\text{los }},{\sigma }_{\text{data }}^{\text{model }}),$$(19)

where 𝒩 represents the normal distribution. The $v_{\mathrm{l}\mathrm{o}\mathrm{s}}^{\text{model }}$ and ${\sigma }_{\text{data }}^{\text{model }}$ are described by two following equations: $$v_{\text{los }}^{\text{model }}={\mathrm{A}\mathrm{m}\mathrm{p}}_{\text{los }}\cdot {\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}^{\mathrm{R}\mathrm{R}}(\mathrm{\Phi })+v_{\text{sys }}$$(20) $${\sigma }_{\text{data }}^{\text{model }}=\sqrt{{\sigma }_{v_{\text{los }}}^2+{({\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}_{\mathrm{e}\mathrm{r}\mathrm{r}}^{\mathrm{R}\mathrm{R}}(\mathrm{\Phi })\cdot {\mathrm{A}\mathrm{m}\mathrm{p}}_{\text{los }})}^2+{({\mathrm{T}\mathrm{e}\mathrm{m}\mathrm{p}}^{\mathrm{R}\mathrm{R}}(\mathrm{\Phi })\cdot {\sigma }_{\text{data }}^{\text{model }})}^2}.$$(21)

In these equations, Φ^k stands for a pulsation phase of the observation calculated using the BH Aur ephemerides and $T_{\text{HJD }}^{\text{APO }}$ in the following relation: $$\mathrm{\Phi }=mod(T_{\mathrm{B}\mathrm{J}\mathrm{D}}^{\mathrm{A}\mathrm{P}\mathrm{O}}-M_0,P)/P.$$(22)

The terms Temp^RR and ${\text{Temp}}_{\text{err }}^{\text{RR }}$ represent the line-of-sight-velocity template (described by the Fourier series; see first row in Table 2) and a spline function that describes the scatter in the template¹⁰.

Following the steps above and applying one of the likelihood minimalization (or log-likelihood maximization) routines (e.g., SciPy and emcee Virtanen et al. 2020; Foreman-Mackey et al. 2013), we estimated υ_sys = 49.7 ± 1.6 km s⁻¹ for BH Aur¹¹. Our derived value based on APOGEE spectra is in agreement with the systemic velocity estimated by Gaia, υ_sys = 51.4 ± 1.5 km s⁻¹.

6.2 Testing the newly derived scaling relations and templates

In this subsection, we assess the accuracy of the derived line-of-sight velocity templates by estimating the systemic velocities for the large fraction of the Gaia RR Lyrae RVS dataset. In order to compile a sample of RR Lyrae stars with υ_sys to use as labels (i.e., the “ground truth”), we used only single-mode RR Lyrae stars and criterion on ipd_frac_multi_peak < 5 to avoid any blended stars. This resulted in a sample of 1083 RR Lyrae stars, and these stars were used to examine the precision of the new scaling relations and the derived templates and their impact on the systemic velocities for RR Lyrae stars with υ_los measurements in the Gaia catalog. For the estimation of υ_sys, we proceeded in the same fashion as in the previous subsection. We used the determined pulsation properties from the Gaia catalog and obtained the systemic velocities for the tested stars (Clementini et al. 2023).

Since some of the RR Lyrae variables in the Gaia catalog have spurious values for υ_los that would negatively affect the υ_sys determination, we imposed a criterion on the υ_los for a given star. First, for a given RR Lyrae star, we estimated the median υ_los value |υ_los|. The value of |υ_los| falls close to the actual υ_sys; thus, in the υ_sys determination, we only considered υ_los values falling within ±50 km s⁻¹ of the |υ_los|. This cut removed nearly all spurious measurements, and the interval of ±50 km s⁻¹ was motivated by the maximum line-of-sight velocity amplitudes (see Fig. 4), where all Amp_los are smaller than 90 km s⁻¹.

In Fig. 9, we depict a comparison between our derived values and the υ_sys measurements in the Gaia RR Lyrae catalog. We found an excellent match for our estimated values and a very low scatter, demonstrating a low RMSE equal to 5.7 km s⁻¹. Several of the outliers in the bottom panel of Fig. 9 (16 variables with differences above |20| km s⁻¹) can be explained by spurious Gaia υ_los measurements that although accompanied by a large ${\sigma }_{v_{\text{sys }}}$ in the Gaia RR Lyrae catalog, still contribute to and influence the Gaia υ_sys. Since we tried removing such measurements, our derived values differ from those published in the Gaia RR Lyrae catalog. An example of an RR Lyrae with spurious measurements is the source_id = 1826975050760883072 (HR Vul). This star has υ_los measurements ranging from 823 km s⁻¹ to −480 km s⁻¹, thus far exceeding the amplitude variation for a classical RR Lyrae star. For HR Vul, we found υ_sys = −8 ± 2 km s⁻¹, while the Gaia RR Lyrae catalog lists υ_sys = −60 ± 50 km s⁻¹. The general Gaia gaia_source catalog provides a value for HR Vul of υ_sys = −4.59 ± 4.35 km s⁻¹¹². This large discrepancy was caused by several spurious measurements of υ_los that were included in the estimation of υ_sys for the Gaia RR Lyrae catalog.

The Gaia RR Lyrae catalog provides numerous measurements of υ_los that, particularly for RR Lyrae stars in the Galactic bulge, are generally hard to come by. Therefore, we performed a more qualitative test of our method as well as a comparison with υ_sys values derived by a method used in the past for Galactic bulge RR Lyrae stars (using line-of-sight velocity template and scaling relation from Liu 1991). We proceeded in the same way as in the case above and only removed first-overtone pulsators since Liu (1991) is calibrated only on fundamental-mode RR Lyrae stars. To simulate possible observations for RR Lyrae stars toward the Galactic bulge, we used data published by the BRAVA-RR survey (Kunder et al. 2016), where most variables were observed with four epochs. We selected the same four epochs for individual stars and applied our methods and those of Liu (1991) to compare the two approaches.

The resulting comparison was made based on the RMSE between Gaia υ_sys values (used as reference values) and values derived using our new method and the procedure described by Liu (1991). Based on 884 tested RRab variables, we found a slightly lower RMSE for our newly derived approach, RMSE(Liu 1991) = 9.6 km s⁻¹ and RMSE^{this work} = 7.5 km s⁻¹. Furthermore, when we included only bright variables (mean apparent G-band magnitudes below 12.0 mag) that have very precise υ_los in the Gaia RR Lyrae catalog in the comparison, we found RMSE(Liu 1991) = 5.8 km s⁻¹ and RMSE^{this work} = 4.0 km s⁻¹ for the method by Liu (1991). Thus, for very precise data, our approach leads to an RMSE smaller by a factor of 1.5 than the RMSE from the method outlined by Liu (1991).

Furthermore, we explored how well both methods estimate uncertainties on υ_sys concerning the selected “ground truth” in Gaia υ_sys values. We first estimated the difference between the υ_sys values derived through two tested methods and the values from the Gaia RR Lyrae catalog, Δυ_sys. Then, for each method, we added uncertainties in quadrature with uncertainties values from Gaia, resulting in σ^tot. Under the assumption that the difference between the selected method and the Gaia values follows the normal distribution, we expected to find ≈68.2 percent of tested stars within one σ^tot of the difference distribution. Thus, from 884 tested RRab variables, 603 pulsators should be within one σ^tot of the Δυ_sys distribution. When compared with two tested methods, we found that our method and that of Liu (1991) yield 597 and 466 pulsators within one σ^tot, respectively. Therefore, our estimated uncertainties quite well represent the underlying error distribution. Viable room for further improvement in our approach would be to include a term for the uncertainty in the star’s ephemerides (e.g., an offset in phase Φ). This term would consider the possible time shifts between photometrically obtained ephemerides and spectral observations. Adding a term for the offset in E could potentially be helpful for RR Lyrae stars with modulated light curves due to the Blazhko effect, where additional modulation can cause shifts in the phase (e.g., RS Boo, see Fig. 6 in Jones et al. 1988). On the other hand, when we compared our results with those from the method of Liu (1991), we found that it underestimates the uncertainties of υ_sys. One of the conceivable reasons can be due to neglecting the scatter in the template, which is included in our procedure.

Fig. 9 Comparison between Gaia-determined systemic velocities (Clementini et al. 2023) and our determined υsys using newly derived line-of-sight velocity templates and scaling relations. The top panel shows a zoom-in of the stars with the difference between Gaia and our values below |20| km s−1, while the bottom panel shows the entire scale. The numbers in each panel denote the total number of stars displayed in a given plot.

7 Summary

In this study, we derived, for the first time, scaling relations between photometric and line-of-sight velocity amplitudes together with line-of-sight velocity templates for visual and near-infrared spectra (covering Ca T and H passband). The aforementioned tools can be used for estimating systemic velocities of single-mode RR Lyrae stars. The main purpose of this study was to derive tools that can be subsequently utilized in studies focused on the structure and kinematics of the Galactic bulge stellar component. In particular, our aim is to provide methods and tools for processing spectroscopic data for RR Lyrae stars in surveys such as the APOGEE and BRAVA-RR (Zasowski et al. 2017; Kunder et al. 2016).

To create scaling relations and line-of-sight velocity templates for the APOGEE survey, we utilized several RR Lyrae variables in the Solar neighborhood. Due to the subvisit spectra, we were able to cover the entire line-of-sight velocity curve for several RR Lyrae stars and construct scaling relations and templates for RRab and RRc subclasses. For scaling equations, we chose to use linear relations without the intercept in order to preserve the physical meaning of the relation. We tested and verified our developed approach for APOGEE spectra on one of the nearby RR Lyrae stars with a precisely determined systemic velocity in the literature.

For the Gaia spectra that cover Ca T, we proceeded in nearly the same way as in the APOGEE analysis. Thanks to the abundance and precise data provided by the Gaia RR Lyrae catalog, we found that the scaling relation for the fundamental-mode RR Lyrae stars is nonlinear and requires a polynomial of a higher degree to accurately describe the transformation of the photometric amplitudes into line-of-sight velocity amplitudes. We did not observe such behavior for first-overtone pulsators. The found nonlinearity was further explored using stellar pulsation models implemented in the MESA code. Using Fourier amplitudes, we were able to replicate the trend for spectroscopic and photometric data. It is important to note that we observed a shift toward larger amplitudes in the low-amplitude regime. The offset in low amplitudes cannot be explained by different physical parameters nor by a varying projection factor. To rectify this shift, we tested different sets of convective parameters, and we found that sets C and D that include turbulent pressure and turbulent kinetic energy flux (see Sect. 4 for details) significantly better describe the low-amplitude space. The effect may be caused by a shift of the IS toward lower effective temperatures. In the low-amplitude regime, the models are significantly cooler than those of sets A and B.

In our analysis, we examined the feasibility of using APOGEE and Gaia survey data interchangeably, focusing on RRab and RRc pulsators. This was driven by the scarcity of calibration RR Lyrae stars in APOGEE. Our comparisons of new and existing templates found a high degree of compatibility for both single-mode RR Lyrae types. Furthermore, we evaluated scaling relations for velocity amplitudes derived in our study and for those in the literature. We conclude that despite the differences between the surveys, they provide similar results, thus enhancing the use and interchangeability of Gaia and APOGEE to derive templates and scaling relations.

We have demonstrated that the templates and scaling relations we provided for RR Lyrae stars in our analysis are largely interchangeable with those found in the existing literature (e.g., Braga et al. 2021), offering a high degree of compatibility. Crucially, our approach distinguishes itself by offering a robust error treatment, notably enhancing the predictability and accuracy of line-of-sight amplitudes for low-amplitude RR Lyrae stars.

We also tested the newly obtained procedure on Gaia Ca T spectra, and our method accurately estimates systemic velocities for RR Lyrae stars in the Gaia RR Lyrae catalog (Clementini et al. 2023). Furthermore, we compared our approach for calculating systemic velocities with methods used in the past, and we found that our approach provides an approximately 50 percent higher precision in the determination of systemic velocities for observations similar to the BRAVA-RR survey strategy. In addition, our approach reliably estimates uncertainties in systemic velocities when compared to the residual distribution between our values and the original Gaia dataset.

Obtaining high-quality spectra of RR Lyrae stars at the distance of the bulge, even when using only moderate resolutions, is an expensive endeavor. Looking toward the future, Rubin-LSST will uncover hundreds of thousands of RR Lyrae stars out to distances more than 50 times the distance to the bulge, so an efficient way of obtaining line-of-sight velocity measurements of distant RR Lyrae stars will also be valuable. The line-of-sight velocity templates and scaling relations provided here allow the systemic velocity of an RR Lyrae star to be accurately found from a single line-of-sight measurement at any point during the pulsation cycle, assuming the stars amplitude, period, and time of maximum brightness are known. As such, our templates open up more precise avenues to carry out large kinematic studies of the bulge RR Lyrae star population as well as beyond. The derived methods for obtaining systemic velocities of RR Lyrae stars can be utilized by the upcoming large spectroscopic surveys such as 4MOST, WEAVE, and MOONS (de Jong et al. 2014; Dalton et al. 2014; Cirasuolo et al. 2014). We will be using our derived scaling relations and templates in forthcoming papers to examine the structure and kinematics of the Galactic bulge in both the northern and southern hemispheres.

Acknowledgements

R.S. was supported by the Polish National Science Centre, SONATA BIS grant, 2018/30/E/ST9/00598. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the US Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss4.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics I Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University. This research made use of the following Python packages: Astropy (Astropy Collaboration 2013, 2018), emcee (Foreman-Mackey et al. 2013), IPython (Pérez & Granger 2007), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), scikit-learn (Pedregosa et al. 2011), and SciPy (Virtanen et al. 2020).

Appendix A Additional tables

References

All Tables

All Figures

Fig. 1 Example of a single exposure spectrum (top panel) and line-of-sight velocities (bottom panel) used in this work. The top panel depicts a spectrum of TU-UMa taken at HJD = 2 457 484.880667 day. The bottom panel shows the line-of-sight velocities derived in this work (black dots) and the velocities extracted from APOGEE data products (red squares). The gray dots represent estimated velocities that did not fulfill the criteria in Eq. (1).

In the text

	Fig. 2 Scaling relations between the photometric and line-of-sight velocity amplitudes for the fundamental-mode RR Lyrae stars (left-hand panel) and first-overtone RR Lyrae variables (right-hand plot). Both RR Lyrae subclasses are depicted using blue points. The best-fit parameters are shown with a solid red line, while gray lines represent their variations.
In the text

Fig. 3 Line-of-sight velocity templates for the RRab (left-hand panel) and RRc (right-hand plots) pulsators. The black points in the top panels represent individual measurements of υlos, and the solid red lines depict the Fourier model for each normalized line-of-sight velocity curve. The bottom panels show the difference, Δ, between the model (Fourier fit) and the phased υlos data for both RR Lyrae subclasses.

In the text

	Fig. 4 Same as Fig. 2 but for line-of-sight velocity amplitudes determined based on Gaia calcium triplet spectra. In addition, the green solid, dashed, and dotted lines represent linear relations from Braga et al. (2021), Sesar (2012), and Liu (1991, using p-factor equal to 1.35), respectively. We note that the literature relations were scaled using Eqs. (11) and (13).
In the text

	Fig. 5 Same as Fig. 3 but for the line-of-sight velocity curves determined based on Gaia calcium triplet spectra.
In the text

	Fig. 6 Relation between the photometric and line-of-sight velocity Fourier amplitudes for RRab stars (gray dots) along with nonlinear pulsation models of [Fe/H] = −1.0 dex (top panel) and [Fe/H] = −2.0 dex (bottom panel) and different masses and luminosities, as indicated in the legend.
In the text

	Fig. 7 Relation between the photometric and line-of-sight velocity Fourier amplitudes for RRab stars (gray dots) along with four sequences of nonlinear pulsation models. We used the following fixed physical parameters (M = 0.65 M⊚, L = 45 L⊚, [Fe/H] = −2.0 dex, with varying effective temperature) and assumed different sets of convective parameters.
In the text

Fig. 8 Comparison between line-of-sight velocity templates derived in this work (top panels, based on APOGEE and Gaia data) and those derived in previous studies (bottom panels, together with APOGEE line-of-sight velocity template). For each line-of-sight velocity template derived in this work, we also included a shaded region (in top and bottom panels) depicting the measured scatter along the template (see Tables A.3 and A.5).

In the text

Fig. 9 Comparison between Gaia-determined systemic velocities (Clementini et al. 2023) and our determined υsys using newly derived line-of-sight velocity templates and scaling relations. The top panel shows a zoom-in of the stars with the difference between Gaia and our values below |20| km s−1, while the bottom panel shows the entire scale. The numbers in each panel denote the total number of stars displayed in a given plot.

In the text