rrlfe: Software for Generating and Applying Metallicity Calibrations for RR Lyrae Variable Stars Across a Wide Range of Phases and Temperatures

RR Lyrae stars play a central role in tracing phase-space structures within the Milky Way because they are easy to identify, are relatively luminous, and are found in large numbers in the Galactic bulge, disk, and halo. In this work, we present a new set of spectroscopic metallicity calibrations that use the equivalent widths of the Ca II K and Balmer H-gamma and H-delta lines to calculate metallicity values from low-resolution spectra. This builds on an earlier calibration from Layden by extending the range of equivalent widths which map between Ca II K and the Balmer lines. We have developed the software rrlfe to apply this calibration to spectra in a consistent, reproducible, and extensible manner. This software is open-source and available to the community. The calibration can be updated with additional datasets in the future.


INTRODUCTION
RR Lyrae stars are old (≳10 Gyr), mostly iron-poor, helium-burning stars that exhibit regular luminosity variations as they pulsate.They are of central importance in tracing out phase-space structures in the Milky Way because they are numerous, relatively luminous (+0.4 ≲ ⟨M V ⟩ ≲ +0.9; Smith 1995), easy to identify based on their light-curve shapes, and are detectable deep into the Galaxy's baryonic halo, to ∼100 kpc from the Galactic centre and beyond.
The weak tidal gradients and shallow gravitational potential of the Galactic halo allow phase-space structures to persist for a few billion years, providing a fossilized eckhart.spalding@sydney.edu.aurecord of the Milky Way's assembly history (Wetterer & McGraw 1996;Amrose & Mckay 2001;Helmi 2008;Iorio & Belokurov 2021).Halo studies were used to investigate this history beginning with Eggen et al. (1962), who concluded that the Milky Way had undergone a rapid collapse that left behind a relatively homogeneous halo.This paradigm was challenged by Searle & Zinn (1978), who found that the heterogeneous metallicities and horizontal-branch morphologies of globular clusters suggested that the history of the halo involved a more protracted hierarchical accretion.Carollo et al. (2007), Carollo et al. (2010), Beers et al. (2012), and others since, have demonstrated the existence of at least two populations of halo stars, which they identified as the inner-halo population and the outer-halo population, with distinct spatial density profiles, stellar orbits, and stellar metallicities.They took this as evidence that the individual halo components likely were formed in fundamentally different ways, e.g., through successive dissipational (inner) and dissipationless (outer) mergers and tidal disruption of proto-Galactic clumps.The emerging picture of the Milky Way's outskirts (≳15 -20 kpc; Deason et al. 2011;Navarro et al. 2021) has indeed confirmed a "lumpy" appearance, with numerous halo over-densities, including some two dozen known stellar streams (Mateu et al. 2017;Bonaca & Hogg 2018;Shipp et al. 2018;Prudil et al. 2021) left over from the tidal disaggregation and accretion of dwarf galaxies.
In addition to being interesting in its own right as a record of the Milky Way's assembly history, halo substructure also provides a probe of ΛCDM cosmology at galactic (∼kpc) scales (Helmi 2008;Johnston et al. 2008;Cooper et al. 2010), where the physics of ΛCDM become complicated by stellar feedback in the form of heating, ionization, and metal enrichment (Springel et al. 2006).The physics to be gleaned may contribute to the resolution of longstanding conundrums such as the "missing satellites problem," whereby the dozens of known dwarf galaxies surrounding the Milky Way are much fewer than the hundreds expected from ΛCDM simulations (e.g., Bullock & Boylan-Kolchin 2017).Further progress on questions relevant to halo substructure and the Galactic potential will require tight constraints on the phase-space, age, and chemical information content of the diffuse halo and multiple stellar streams (e.g., Sanderson et al. 2017a;Ramos et al. 2020).
RR Lyrae stars are located at the intersection of the instability strip and the Horizontal Branch on the HR diagram.As lower-mass counterparts to Cepheid variables, RR Lyrae stars have a long history of use as standard candles for mapping the structure, kinematics, and formation history of the Milky Way (e.g., Vivas et al. 2008;Braga et al. 2015;Pietrukowicz et al. 2015;Minniti et al. 2016;Dong et al. 2017;Mateu et al. 2017).They include the subtype ab, which pulsates in the fundamental mode; subtype c, which pulsates in the first overtone mode; and subtype d, which pulsates in both modes.The ab and c subtypes (RRab and RRc stars) occupy distinct regions of period-amplitude space, whereby RRab stars are identifiable by qualitatively "shark-tooth" light curves, pronounced atmospheric shock waves, and cooler effective temperatures.RRc stars have a more sinusoidal light curves, and hotter effective temperatures (Catelan & Smith 2015).
The dominant drivers of RR Lyrae pulsations are the effect of the increasing opacity as a function of temperature at certain ionization zones, particularly that of the second ionization layer of He (the κ-mechanism) and the increased flow of heat into those same zones (the γ-mechanism).Within the instability strip, RR Lyrae stars rove back and forth on the HR diagram over a temperature range of T eff ∼ 6000 K to 7250 K (Catelan & Smith 2015).
In addition to their utility as standard candles, RR Lyrae stars have served as markers of metallicity.Given the long lifetimes of RR Lyrae stars, their iron abundances constrain the timescales, relative chronologies, and variegated parenthoods of stellar populations, in ways that are complementary to 6D phase-space information.
Within the temperature range of RR Lyrae stars, the Balmer lines become stronger with higher temperatures, as more electrons are available in the first excited states for the absorption of photons, and temperatures are not hot enough for the effects of ionization to dominate.Thus, the Balmer-line equivalent widths (EWs) are strongly correlated with effective temperature.Over the same temperature range, the Ca II lines weaken with increasing temperature as they ionize, but also possess some correlation with iron abundances (Gray et al. 2009;Fernández-Alvar et al. 2015).The inverse behavior of the changing Balmer and Ca II lines, and the partial correlation of Ca II K with iron, offers the possibility of calculating RR Lyrae iron abundances with multi-epoch spectroscopy.Preston (1959) was the first to note that the metallicity of RR Lyrae stars appeared to correlate with differences in spectral type, as determined from the Ca II K and hydrogen Balmer lines at minimum light.Later updates to the method of Preston (e.g., Butler 1975;Freeman & Rodgers 1975) led to the calibration of Layden (1994), henceforth L94, which was based on the pseudo-equivalent widths of the Ca II K and Balmer Hδ and Hγ lines.However, since the L94 calibrations were based on spectra at low-temperature phases (and thus narrow Balmer-line EWs), the calibration breaks down at hotter temperatures.Furthermore, the L94 calibrations with interstellar calcium correction were based on a model based on distance from the Galactic plane, symmetric in Galactocentric azimuth (Beers 1990).Since then, it has become increasingly clear that interstellar Ca II absorption varies significantly along different linesof-sight (e.g., Welty et al. 1996).Thus, the time is ripe for a new calibration which both maximizes the range of usable light-curve phases, takes a new look at interstellar Ca correction, and will be applicable to large spectroscopic datasets.In addition, none of the calibration stars used for setting iso-metallicities was more metal rich than [Fe/H] = −0.34,as calculated in L94.A calibration that includes more metal-rich stars would be applicable to many additional stars in the halo and disk systems.Given the range in effective temperatures of these stars, a new relation between EWs can also be made to allow for the non-linear behavior of the lines.
Fortuitously, RR Lyrae stars are plentiful in the Sloan Digital Sky Survey (SDSS) dataset (Sesar et al. 2007), which has yielded millions of low-resolution spectra since the first data release in 2003 (Abazajian et al. 2003).
Though not designed as a variable star survey (Castander 1998), a subset of perhaps thousands of SDSS spectra represent randomly phased, multi-epoch spectra of RR Lyrae stars.Given the survey sensitivity, RR Lyrae stars are detectable in SDSS data out to distances ≳100 kpc (Newberg 2003;Ivezić et al. 2005;Sesar 2011), and have been used since the early days of the SDSS to trace out the physical and radial velocity structure of halo over-densities (e.g., Ivezić et al. 2000;Sesar et al. 2007;Bell et al. 2008;Watkins et al. 2009;Sesar et al. 2010;Simion et al. 2014;Casetti-Dinescu et al. 2015).
The SEGUE Stellar Parameter Pipeline (SSPP) was designed to calculate, among other things, stellar metallicities based on SDSS spectra (Lee et al. 2008a,b;Allende Prieto et al. 2008;Smolinski et al. 2011).However, the SSPP assumes LTE conditions, which are not always a fair approximation in pulsating variable stars.
Photometric RR Lyrae metallicity calibrations also exist in the form of metallicity-period-ϕ 31 relations, where ϕ 31 is the Fourier phase parameter relation (Jurcsik & Kovács 1996;Sandage 2004;Ngeow et al. 2016;Hajdu et al. 2018;Li et al. 2023).Our spectroscopic calibration will help provide a basis for photometric metallicity calibrations, which will undergo renewed importance in the era of the photometric Legacy Survey of Space and Time (LSST) at the Vera C. Rubin Obser-vatory.After routine science observations begin in the near future, LSST will detect millions of RR Lyrae stars out to hundreds of kpc to even a few Mpc, to the farthest outskirts of the Galactic halo and even beyond to other members of the Local Group (Oluseyi et al. 2012;Baker & Willman 2015;Sanderson et al. 2017b).It will be impractical to obtain follow-up spectroscopy on such a large number of stars, especially since LSST goes down to r ∼ 24.5 for single visits and for some observations down to r ∼ 26 (Abell et al. 2009).Thus, a recalibration of photometric techniques will be important.
Here we present an updated metallicity calibration for RR Lyrae stars as generated from a basis set of synthetic spectra and empirical RRab spectra, and provide the software for applying this calibration to other stars and for revising this calibration in the future.To generate the calibration, our software normalizes the synthetic spectra, measures the EWs of the Ca II K, Hγ, and Hδ lines, and generates a functional fit.A final correction is applied, based on the comparison of retrieved values of [Fe/H] with a sample of RR Lyrae stars with a wide range of metallicities, as determined from highresolution spectroscopy from previous studies.
To obtain [Fe/H] estimates with our pipeline, we collected multi-epoch spectra for stars at R ∼ 2, 000 to mimic those in SDSS, as spectra for our program stars do not exist in SDSS as of DR14.We phase them using photometry from TESS, KELT, and other sources, so as to determine where in the phase the effects of shocks may cause the calibration to break down.The resulting fits can be applied to single-epoch RR Lyrae spectra from SDSS or other surveys.
The organization of our paper is as follows.We describe our methods in Sec. 2, and the reduction of observations in Sec. 3. Results are described in Sec. 4, which includes a test of our calibration on empirical spectra and the well-studied SDSS Stripe 82 in Sec.4.3.We provide a discussion of the implications in Sec. 5, and summarize in Sec. 6.

Synthetic Spectra
We computed a grid of synthetic spectra which covers the range in expected physical parameters for RR Lyrae stars.The grid spans 5750 K ≤ T eff ≤ 7750 K, in increments of 250 K, 2.0 ≤ log g ≤ 3.0, in increments of 0.5 dex, and −2.5 ≤ [Fe/H] ≤ 0.0, in increments of 0.5 dex.The synthetic spectra were generated using Castelli & Kurucz (2003) Atlas9-ODF, LTE model atmospheres with microturbulence of 2.0 km s −1 , and the spectral synthesis code SPECTRUM v2.76 (Gray & Corbally 1994).The synthetic spectra were generated at high resolution and Gaussian smoothed to a resolution of 3.2 Å with a dispersion of 1.4 Å per pixel.Since we are using the Ca II K line EW as a proxy for Fe abundance, we scaled the calcium abundance at low metallicity using the relation from Eqn. 2 in Allende Prieto et al. (2006).The spectra span a range of 3900 Å ≤ λ ≤ 5000 Å, which included the lines of Ca II K and H, and the Balmer lines of Hβ, Hγ, Hδ, and Hϵ.These spectra were then processed by our custom pipeline rrlfe, which we describe in Sec.3.1.1.

Literature Metallicities
Estimates of [Fe/H] based on synthetic spectra must be compared with high-resolution studies of actual stars.We selected a set of 19 program stars for comparison and observation, including both RRab and RRc stars, most of which have metallicities calculated from highresolution studies.These literature sources of [Fe/H] are listed in the appendix in Table A1.We derived a net [Fe/H] value from multiple appearances of each star in the literature, following the general method of Chadid et al. (2017), which we summarize as follows.
Metallicity values from a given high-resolution study are plotted against those of L94, for stars shared between both works.L94 included only RRab stars and was not a high-resolution study, but involved over 300 RR Lyrae stars, increasing the probability of sample overlap.The offset between a line of best fit at an arbitrarily-chosen abscissa value of [Fe/H] L94 = −1.25, and the ordinate of the linear fit to the metallicities of 26 program stars of Chadid et al. (2017), is added in to shift all the highresolution [Fe/H] values from a given study.After this is done for each of the high-resolution studies, this acts to remove systematic differences between them and maps them onto a common scale, while preserving the intrinsic scatter in metallicity among stars within each study as compared to L94.Fig. 1 shows these shifted literature values against L94.
We considered using a second-order polynomial or a piecewise-linear function to fit the metallicities from high spectral-resolution studies as a function of the L94 values.However, we found that the scatter at the metalpoor end did not provide better fits, at the price of an additional degree of freedom.This left us with a net literature metallicity value [Fe/H] lit based on the linear fit [Fe/H] lit = A [Fe/H] L94 + B, with A = 1.040 ± 0.026 and B = −0.060± 0.035.Metallicities mapped in this way, along with their literature values, are tabulated in Table A3.It should be noted that this mapping is generated using only RRab stars.To our knowledge, the largest single compendium of RRc metallicities is Kemper (1982), which used the ∆S method on image We therefore refrain from making another metallicity mapping for RRc stars alone, and apply any resulting calibration to them as-is.

Photometry
The very pulsation that causes effective temperature changes in RRab sends shocks through the photosphere, which can perturb the line profiles (Chadid & Gillet 1996;Fokin et al. 1999;Chadid et al. 2008;Preston 2011;Chadid et al. 2017).To put a constraint on this region of the phase in our program stars, as well as to make a general prediction of the phase for planning purposes, we initially used data from the website of the American Association of Variable Star Observers (AAVSO) (Kafka 2021).To minimize the effect of drift in predicted phase we obtained follow-up observations from December 2012 to January 2014, from three locations: the MacAdam Student Observatory (MSO) on the University of Kentucky campus in Lexington, Kentucky, the University of Louisville's Moore Observatory (ULMO) near Louisville, Kentucky, and in one case, the University of Southern Queensland's Mt.Kent Observatory, near the town of Toowoomba in eastern Australia.
All three telescopes were 20-inch PlaneWave Corrected Dall-Kirkham f/6.8 models (CDK20).Filters used were V , G, and in one case, R. The Moore data was reduced using AstroImageJ (Collins et al. 2017).Fourth-order polynomials were fit to the relevant parts of the light curves in order to determine an epoch of maximum in JD and calculate the phases of the spectroscopy.
The photometry for a given star from MSO, LMSO, and Mt.Kent were from a single night to capture the maximum light, with a dense sampling cadence of ≈ 20 sec to 2.5 min, so as to establish the time at which the pulsation phase was zero.Pixel scales at these telescopes were ≈ 0.4 -1"/pixel.
To make the most precise determinations of phases following acquisition of the stars' spectra, we also obtained photometry from the Kilodegree Extremely Little Telescope (KELT) and the space-based Transiting Exoplanet Survey Satellite (TESS) (Pepper et al. 2007;Ricker et al. 2015).The KELT and TESS observations use wide fields of view of 26 • × 26 • and 24 • × 24 • from a single camera, respectively, larger pixel scales of ∼23"/pixel and ∼21"/pixel, and observing baselines stretching over months or years, at typical sampling cadences during observing windows of ≈ 2-30 min (Pepper et al. 2012).The long temporal baselines of KELT and TESS, and the extra benefit of the excellent photometric precision of TESS, enabled the dense phase sampling required for establishing highly precise pulsation periods.All told, 13 stars were observed at McDonald Observatory in December 2012, and 6 in July 2013, totaling 18 different stars, including RU Psc, which was observed during both runs.Of these stars, 13 were RRab stars.A total of 169 useable spectra were reduced in the analysis (see Fig. 2).The stars whose periods were covered as much as possible were VX Her (an RRab, with 26 spectra) and RU Psc (an RRc, with 22 spectra).
Raw spectra from December 2012 ranged over wavelengths of 3900 Å to 5000 Å, and in July 2013, 3610 Å to 4960 Å. Spectra were later truncated to the range 3911 Å to 4950 Å, to limit distortions of the Ca II K and Hβ absorption lines during the normalization process.Integration times ranged from 40 to 1800 seconds, and always remained less than 10% of the period in order to avoid smearing of the phase.Of the 169 spectra, the integration times of six of them represented less than 1% of the period, 52 represented between 1% and 2%, 86 between 2% and 5%, and 24 between 5% and 9%.In a few instances, integration was paused due to a passing  Pancino et al. (2015).Red lines are spline fits from the folding process, and vertical blue lines indicate the phase epochs of our spectra.Vertical axes are scaled to fit the spline fits, leaving some photometric points outside the display region, so as to make it easier to see the curves.RU Psc, RW Ari, and V535 Mon have poorly defined periods (e.g., Wils 2008;A. Odell, pers. comm.).Photometry for those stars were limited to short baselines at times near or around the time the spectra were taken.Offsets in the BH Peg subplot may be due to astrophysical Blazhko changes (period ∼10s to ∼100 days; Smith 1995; Skarka 2014) over the 6-year span of the AAVSO observations, or changing photometric zero points in the observations.Overlaid plots at the bottom centre and right exclude RU Psc and V535 Mon.(See Appendix B and Sec. 2.4.) cloud, which may have increased the total elapsed time by a minute or two.

Calculation of Periods and Phases
After the spectroscopy were gathered at McDonald Observatory, we re-calculated periods for most of the stars using phase-folded photometry from KELT and TESS.The heavy sampling of photometry in these long- term surveys maximizes the precision of the periods, which minimizes the error in the calculated phase based on a given epoch-of-maximum.Two of our program stars were unavailable in KELT or TESS: for VX Her we obtained photometry from the AAVSO, and for BH Peg we used the period in Monson et al. (2017).
The KELT light curves used dates in terrestrial time, so the dates were converted to barycentric Julian dates (BJD) using astropy (Astropy Collaboration et al. 2013, 2018).TESS light curves were downloaded from the online Mikulski Archive for Space Telescopes (MAST), and were already in BJD.
The VARTOOLS light-curve analysis program (Hartman & Bakos 2016) was used to provide a common interface as we tried several different period-finding algorithms.We chose to use the Fast-χ 2 algorithm (Palmer 2009), limiting our search to between 0.1 and 0.9 days, with the number of harmonics set to 3 (fundamental plus the first two harmonics).The periodogram was over-sampled by a factor of 10 and the two most significant peaks were analyzed.For many stars in the KELT footprint, separate curves taken from both east and west of the meridian were available.In such a case, the period we adopted for the rest of the analysis corresponded to the best periodogram peak in the best-fitting light curve (East or West) as determined by eye.We compared these to the published periods from the General Catalog of Variable Stars (GCVS, Samus' et al. 2017) to check for consistency (Fig. 3).
For BH Peg, we used the period as determined by Monson et al. (2017), who used archival data with a baseline of thirty years, together with robotic observations from the Three-hundred MilliMeter Telescope (TMMT) at Las Campanas Observatory.The data were phase folded across multiple basis wavelengths, and transformed onto a basis filter set of Johnson B, V , Kron-Cousins I C , 2MASS J, H, K S , and Spitzer [3.6] and [4.5] µm channels.We used the period from Monson et al. (2017) to phasefold the AAVSO photometry, which we limited for self consistency to observations contributed from 2015 to 2021 by one observer, G. Samolyk.We found that the folded light curve was the most consistent of any produced with other literature period values, or period values found from our own phase folding of sparsely sampled AAVSO data.It became clear that not all of our program stars have well-defined periods.The KELT light curves for RU Psc and V535 Mon showed significant variations in the pulsation period, the envelope amplitude of the pulsations characteristic of Blazkho variations, or both, to the extent that they could not be phased with the entirety of the available photometry.We therefore restricted the photometry to a short baseline, closest in time to when the spectroscopy was collected.For RU Psc we used the East KELT light curve limited to BJDs 2456187.75 to 2456663.58.For V535 Mon we limited the West light curve to BJDs 2456281.76 to 2456304.98.In addition, RW Ari is known to have an irregular pulsation (A.Odell, pers.comm.).
To determine the epoch of maximum brightness, we used the derived periods to phase the light-curve data.The phased light curves were then binned into 100 phase bins and the median of each bin was determined.In the cases of V535 Mon, RU Psc, and BH Peg, 20, 50, and 65 bins, respectively, were used due to a smaller number of data points.The binned data was then fit with a 3rd order univariate spline from the SciPy interpolate package, and the brightest point on the spline curve was taken to be the maximum.The phase of maximum was then converted back into an epoch and used to re-phase the light curve, so that maximum light occurred at zero phase.The periods of our program stars were used to The sources of photometry and final periods are listed in Table 1.

Reduction of synthetic spectra
We wrote a pipeline in Python to take the synthetic spectra and generate a higher-order version of the calibration of L94, using the Python version of Robospect to measure EWs (Waters & Hollek 2013).
Our expanded version of the L94 relation is where W (K) is the EW of the Ca II K line and W (H) is the EW in angstroms of the Balmer lines.Coefficients {a, b, c, d} are used in L94, and our addition of {f, g, h, k} is explained in the Appendix B. The Hβ line EWs were discarded from the analysis because they introduced considerable scatter in both the synthetic and empirical spectra, apparently due to metal contamination, as was also noted by L94 (Fig. 4).
We checked for consistency in the EWs returned by Robospect and manual use of IRAF (Tody 1986).1With IRAF we tested window sizes of 10, 28, and 56 Å on the synthetic spectra, and found results between Robospect and IRAF to be consistent with window sizes set to 56 Å.
The steps of our pipeline workflow are as follows: 1. Normalize the spectra 2. Measure EWs and errors for Ca II K, Hδ, and Hγ using Robospect CCD readouts of the spectra taken at McDonald Observatory were reduced with IRAF.Finding a solution to the shorter-wavelength end of the spectrum proved challenging, due to the presence of no more than two reliable calibration-lamp emission lines at wavelengths less than about 4150 Å (which includes the Hδ and Ca II K absorption lines in the object spectra), and a complete lack of lines below 3950 Å (which includes Ca II K).
Spectra were written out to ASCII files, and cosmic rays were manually removed by deleting data lines in the ASCII file.Fig. 5 shows example spectra taken of X Ari, including the time at which a shock travels through the photosphere.For each star, we assigned the metallicity found from the mapping of [Fe/H] from high-resolution studies as described in Sec.2.2.

Contamination from Interstellar Calcium
The presence of interstellar calcium along the line-ofsight to an observed RR Lyrae can add to the total EW of the measured Ca II K line, and systematically increase the star's computed metal abundance.To compensate for this contamination, L94 applied the interstellar calcium EW model of Beers (1990), which is a function of a star's Galactic latitude and vertical height above the Galactic plane.This model allowed for the subtraction of the EW of the interstellar component from the EW of the net Ca II K line.
However, SDSS spectroscopy since then has revealed highly textured interstellar Ca II K absorption (Murga et al. 2015) within a large observational footprint that contains some, but not all, of our program stars.Even if the line-of-sight ISM calcium absorption to a single star can be quantified precisely-and the resolution of the map in Murga et al. (2015) is about 3.7 • ×3.7 • -the saturated Ca II K stellar line will experience shifts of up to ≈ 70 km s −1 , in addition to the intrinsic radial velocity of the star.Thus, the ISM will contaminate the stellar Ca II K line different amounts at different times.
To see whether it is worth modeling these effects, we investigated the effect of the interstellar calcium on our most metal-poor program star, X Ari.The interstellar Ca II K line is readily seen redward of its stellar counterpart in the high-resolution spectrum in Fig. 6.At this resolution, the two lines can easily be deblended.When using low-resolution spectroscopy (R ≈ 2000), however, the contribution of the interstellar component is blended into the line and cannot be directly removed.
We measured the EW of the interstellar Ca II K line in the high-resolution spectroscopy from Nemec et al. (2013) to be 80 m Å. Inclusion of this contaminant increases the net Ca II K EW by 4.6%.After smoothing the spectral resolution to that of R = 2000, the increase of the EW (which includes the interstellar component) as measured with a Gaussian fit, over the highresolution, deblended stellar EW, is 1.2%.This is well below the random fitting error of the low-resolution line.
We also used the maps of Murga et al. (2015) to find upper bounds to ISM Ca II K absorption close to the lines-of-sight towards our other program stars.The maps were generated from the integrated absorption along lines-of-sight to distant quasars and galaxies, and thus represent an over-estimate of the absorption to our stars, which are at 0.25 kpc < d < 2.1 kpc (Gaia Collaboration et al. 2021).Fig. 7 compares Ca II K EWs of our stars with values of ISM absorption, which we take to be the maximum of the ISM absorption within an angle of 7.4 • (two resolution elements) of the star.For the stars which appear in or near the footprint, the upper bounds of the ISM are never larger than the smallest Ca II K EWs, and are an average of 0.18 Å of the minimum Ca II K EW for each star, and 0.08 Å of the maximum.
Based on the high-resolution spectrum of X Ari and the upper bounds on the ISM absorption, we conclude that the measured Ca II K EWs of RR Lyrae variables with [Fe/H] > −2.5 will be essentially unaffected by  interstellar calcium, due to the significantly larger EW of the stellar line, regardless of pulsation phase.The overall small and variable contribution to the stellar EW as measured from low-resolution spectroscopy justifies our decision to not attempt to correct for the interstellar component in our calibration.

Coefficients
Fig. 8 shows a plot of the measured EWs from the synthetic spectra and iso-metallicity contours from the fit.The original L94 calibration implied a Ca II K EW of zero at a Balmer-line EW of ≈ 12 Å for [Fe/H] = −1.Our calibration is comparatively 'shallower,' and can handle a larger range of Balmer-line EWs.Errors in the EWs as measured by Robospect were larger than would be expected based on the scatter evident in Fig. 8, so we fit a linear scaling relation between the Robospect errors and the standard deviation in measured EWs from the star VX Her, for which we obtained particularly heavy phase coverage with our spectroscopy (Fig. 2).We subdivided the spectra of VX Her into phase bins of ∆ϕ = 0.25, found the standard deviations therein, and compared it to the mean of Robospect errors for the same lines.Based on this exercise, we implement in our pipeline a rescaling of the errors returned by Robospect, and consider this the true error in EW.

Tests on McDonald Observatory Data
This calibration was tested on McDonald Observatory spectra of both RRab and RRc stars.Out of 169 spectra, 9 had line-fit failures.Fig. 9 shows the measured EWs, and Fig. 10 shows the results of the metallicity estimates, compared to what one would expect from high-resolution studies.After visual inspection, the badphase region was found to be so narrow in phase given our sampling that it was disregarded.
There is an over-estimation of [Fe/H] at the metalpoor end, and an under-estimation at the metal-rich end.To counteract this, we implemented a final correction in the pipeline to effectively map the best-fitting line onto the one-to-one line.

Tests on SDSS Spectra
This calibration was tested on single-epoch SDSS spectra of unique stars that appear in the footprint of the SDSS Stripe 82 or the photometric Catalina Sky Sur-vey (Larson et al. 2003;Drake et al. 2013;Abbas et al. 2014).Stripe 82 is a subsection of the SDSS footprint along the celestial equator that was imaged 70 to 90 times between 1998 and 2007 (Jiang et al. 2014).After manually removing cosmic rays, Stripe 82 spectra were processed in the same way as the McDonald spectra.
Fig. 11 shows that the SDSS spectra have a clear peak in at metallicity [Fe/H] ≈ −1.6, with a metalpoor tail which is thicker than the metal-rich tail, possibly because of a lower-metallicity stellar component at [Fe/H] ≈ −2.2.This is consistent with the metallicities assigned to an inner halo (∼ 10-15 kpc from the Sun) and an outer halo (beyond ∼ 15 kpc from the Sun), respectively (Carollo et al. 2007;Carollo et al. 2010;Beers et al. 2012).Fig. 12 shows comparisons of our estimates of [Fe/H] based on SDSS spectra, against the SSPP pipeline and values from other recent studies.In the case of the SSPP, the persistent offset may be in part due to the fact that SSPP metallicities depend more heavily on the Ca II K line as the metallicity decreases; and as the metallicity increases more techniques can be applied the find the abundances.

Effective Temperature
We also considered estimates of T eff by using the net Balmer-line width.This relation is highly linear, and a simple fit leads to the relation: where A = 192.8± 1.6 K/ Å and B = 5434 ± 11 K.We apply this to synthetic spectra in Fig. 13.The esti-  Red squares are average values taken across all the spectra for that star, and are the same as the retrieved values listed in the appendix in Table A3.
mated temperatures are sufficiently accurate that "temperature curves" can be obtained, as shown in Fig. 14.

DISCUSSION
The calibration we present in this work, and the pipeline used to generate it, is available as a software package (Sec.7).It can be applied to SDSS survey spectroscopy, and to refine photometric calibrations such as those that will make use of the LSST dataset.It would also be of interest to find where consistency among these different methods fails.
One advantage of our methodology is that the lowresolution input spectra can be restricted to a range of ∼ 3910 -4950 Å, or a bandwidth of ∆λ ≈ 0.1µm.Our methodology is also not limited to RR Lyrae stars.It can also be used to generate and apply calibrations to main-sequence non-pulsating stars, until the Balmer lines fade in strength at later types.This is not only relevant for studies of Milky Way structure, but also investigations further afield, such as correlations between exoplanet demographics and the metallicities of their host stars.The latter provide clues to the physics of the planet-formation process (Kolecki et al. 2021;Boley et al. 2021;Ghezzi et al. 2021), and gleaning metallicities from low-resolution host-star spectra is particularly them all up with high-resolution spectroscopy with large telescopes.
There are a number of possible avenues to improve our calibration in the future.One would be to include more synthetic spectra in the metal-poor ([Fe/H] < −2.5) and metal-rich ([Fe/H] > −1.0) regimes.This will provide a better lever arm for making refinements based on comparisons with the SSPP pipeline.As a secondary priority, it would be useful to obtain empirical, low-resolution spectra that cover more of the pulsation phase of a given sample of stars.In this work we found consistent metallicities across almost all phases tested, but with more phase coverage the phase region where this calibration breaks down could be constrained further.
To expand the sample of input stars with highresolution spectroscopy for mapping literature [Fe/H] values, one may consider stars accessible from the Southern Hemisphere, so as to overlap with the 11 field RRab stars of For et al. (2011) (see also Chadid & Preston 2013).Their stars have unmatched phase coverage with high-resolution (R ≈ 27, 000) spectroscopy, and the calculation by For et al. (2011) of stellar parameters like log g would allow us to try to find correlations between them and the metallicities generated by our calibrations.
In this work we have also obtained low-resolution spectra of RRc stars and applied our calibration thereto, but a better comparison of metallicity estimates is needed for more RRc stars with a set of [Fe/H] values from high-resolution spectroscopy.A larger basis set of RRc  stars would also allow this subtype to be folded in to the generation of a calibration, or enable a separate RRc calibration.High-priority RRc star targets for lowresolution spectroscopy include those which already appear in high-resolution studies: U Com, T Sex, DH Peg, YZ Cap, and eight southern RRc stars in the All Sky Automated Survey (ASAS) (Govea et al. 2014).Target lists that expand beyond these, and which overlap among multiple studies, will be useful for making a mapping specific to RRc stars analogous to our mapping of RRab stars to those in L94.

SUMMARY
In this work, we have established a metallicity calibration for low-or mid-resolution spectroscopy of RR Lyrae variables from a spectroscopic basis set including both RRab and RRc sub-types.This includes a final correction based on comparison with estimates based on high-resolution spectroscopy in the literature.We make our "living" code available, which can be updated with new data, including for RRc stars.
Acknowledgements: E.S. acknowledges support from Huffaker Travel Scholarships, and department Chair Sumit Das for providing additional travel funding, under the auspices of the University of Kentucky Physics and Astronomy Department.T.C.B. acknowledges partial support for this work from grant PHY 14-30152; Physics Frontier Center/JINA Center for the Evolution of the Elements (JINA-CEE), awarded by the U.S. National Science Foundation.The work of V.M.P. is supported by NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.Y.S.L. acknowledges support from the National Research Foundation (NRF) of Korea grant funded by the Ministry of Science and ICT (NRF-2021R1A2C1008679).Y.S.L. also gratefully acknowledges partial support for his visit to the University of Notre Dame from OISE-1927130: The International Research Network for Nuclear Astrophysics (IReNA), awarded by the US National Science Foundation.We also acknowledge Rhodes Hart, who obtained a light curve for one of our program stars from Mt. Kent Observatory in Australia, which is administered by the University of Southern Queensland.
People who provided invaluable assistance for the observations included Kyle McCarthy and Tim Knauer at the MacAdam Student Observatory, and David Doss and John Kuehne at McDonald Observatory.We also thank the University of Arizona's and CyVerse's computing staff members Blake Joyce and Julian Pistorius for their computational assistance.The logo for rrlfe was created by Anna McElhannon.
Finally, we acknowledge Jim Nemec, Thomas Gomez, and George Wallerstein for kindly providing their highresolution spectroscopy to constrain line-of-sight interstellar calcium absorption; Chris Waters, for his development of the Python version of Robospect; and Andrew Odell for pointing out the irregular pulsation of RW Ari.
We acknowledge with thanks the variable star observations from the AAVSO International Database contributed by observers worldwide and used in this research.In particular, we are indebted to observer Gerard Samolyk for his observations of BH Peg and VX Her.Other observers who contributed to the heavy sampling of the VX Her light curve were Teofilo Arranz, Jerry Bialozynski, Neil Butterworth, Shawn Dvorak, Massimiliano Martignoni, Kenneth Menzies, Tom Polakis, John Ritzel, John Ruthroff, Richard Sabo, Neil Simmons, and Denis St-Gelais.This paper also makes use of data from KELT (DOI 10.26133/NEA8).Early work on KELT-North was supported by NASA Grant NNG04GO70G.Work on KELT-North was partially supported by NSF CAREER Grant AST-1056524 to S. Gaudi, and from the Vanderbilt Office of the Provost through the Vanderbilt Initiative in Data-intensive Astrophysics.
Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England.The SDSS Web Site is http://www.sdss.org/.
The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions.The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington.
This paper includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST).Funding for the TESS mission is provided by NASA's Science Mission directorate.This research has also made use of the SIM-BAD database, operated at CDS, Strasbourg, France (Wenger et al. 2000).
This research was supported in part by the Notre Dame Center for Research Computing's computing clusters.Some computations were also performed on Cy-Verse Atmosphere cyberinfrastructure, which is supported by the National Science Foundation under Award Numbers DBI-0735191 and DBI-1265383, URL: www.cyverse.org.Jeff Chilcote also kindly allowed the use of his planetfinder cluster.
Author contributions: E.S. led the acquisition and analysis of the spectroscopy and photometry from Mc-

B. FUNCTION FITTING
The mapped metallicities for our program stars (see Sec. 2.2 and Table A3) and the Ca II K and Balmer EWs were fed into the MCMC package emcee to fit a variation of the L94 relation: We fit an expansion of the relation up to third order.Using the shorthand F for [Fe/H], H for W(H), and K for W(K): with coefficients {a, b, c, d} all non-zero, and all possible zero/non-zero permutations of the extra six coefficients {f, g, h, k, m, n}.The number of permutations is the total of all choices of "6 pick k, for k from zero to 6": where is the number of free parameters, n is the number of data points, and L is the model likelihood function.
We have taken L ∝ exp(− 1 2 χ 2 ), where χ 2 is the sum of the weighted squared residuals between the measured equivalent width of the Ca II K line and that retrieved by using the model Eqn.B2 with the best-fitting coefficients.The sum is taken over spectra i, weighted by the error σ K0 in the Ca II K line returned by Robospect: where K 0 is measured, and K m,i is the model value.This value depends on the Balmer-line EW and metallicity, or x i x i x i = {H i , F i }; θ θ θ is the vector of coefficients as defined in Eqn.B2.
We fit models representing every permutation in Eqn.B3 and calculated the BIC.The difference ∆BIC between one model and L94 compares the fits while penalizing extra model parameters.models that converged, and were an improvement on L94.We found the best-fitting model to be Eqn.B2 with m = n = 0 and all other coefficients non-zero.
The solution to Eqn.B2 with m = n = 0 for the metallicity F is simply the the quadratic formula: where (B7) Eqn.B6 has two solutions, of which the '+' solution in the numerator is the physical one.
2.4.Spectroscopy The spectroscopy was gathered with 2.08 m Otto Struve Telescope at McDonald Observatory in Texas, using the low-to-moderate resolution Electronic Spectrograph Number 2 (ES2), with R ∼ 1, 000.The two observing runs took place on UT 2012 December 21-28, and UT 2013 July 20-26.

Figure 2 .
Figure 2. The phase-folded light curves of our program stars, in the style of Fig. 1 ofPancino et al. (2015).Red lines are spline fits from the folding process, and vertical blue lines indicate the phase epochs of our spectra.Vertical axes are scaled to fit the spline fits, leaving some photometric points outside the display region, so as to make it easier to see the curves.RU Psc, RW Ari, and V535 Mon have poorly defined periods (e.g.,Wils 2008; A. Odell, pers.comm.).Photometry for those stars were limited to short baselines at times near or around the time the spectra were taken.Offsets in the BH Peg subplot may be due to astrophysical Blazhko changes (period ∼10s to ∼100 days; Smith 1995; Skarka 2014) over the 6-year span of the AAVSO observations, or changing photometric zero points in the observations.Overlaid plots at the bottom centre and right exclude RU Psc and V535 Mon.(See Appendix B and Sec.2.4.)

BFigure 3 .
Figure 3. Periods we use for our program stars, compared to those in the GCVS 5 catalog of variable stars.The dashed line is one-to-one in the top plot, and zero in the bottom plot.Vertical lines connect the same star to guide the eye.The bottom plot is in units of 10 −4 day, or 8.64 sec.The residual of the period of V535 Mon is outside the plot, at −70.734 × 10 −4 day = 611 sec.

Figure 4 .
Figure 4.A comparison of the EWs of hydrogen lines with that of Hδ (filled data points), showing the strong linear trend between Hγ and H δ .For clarity, Hϵ and H β have been offset, and error bars have been removed.

3.
Generate a net Balmer EW from a simple average between the Hγ and Hδ EWs 4. Feed the EWs, errors, and input metallicities into an MCMC that samples the posteriors of the coefficients and their errors as per the expanded version of the relation in L94, using the emcee package (Foreman-Mackey et al. 2013) 5. Save the sampled posteriors for application to empirical spectra 3.1.2.Reduction of empirical spectra

Figure 5 .
Figure5.Normalized spectra of X Ari at different phases.Spectra are offset for clarity, with the spectrum at phase 0.4 underplotted in grey for comparison.The pulsation shock at late phases leads to a filling-in of lines.

Figure 6 .
Figure 6.High-resolution spectroscopy from Nemec et al. (2013) of the metal-poor star X Ari.Lines are identified with the model atmosphere program SPECTRUM (Gray 1999).The fortuitous pulsation phase of the star allows separation between the stellar and interstellar Ca II K line.Data kindly provided by J. Nemec.Also see Sec. 3.2.

Figure 8 .
Figure 8.The space of Ca II K and Balmer EWs of the synthetic spectra, with iso-metallicity contours from the final calibration.Points with error bars show the median sizes of the error bars at selected metallicities.The S/F ratio indicates the number of spectra with successful and failed line fits.

Figure 9 .
Figure9.The space of Ca II K and Balmer EWs for the program stars, with the same contours and relative [Fe/H] marker size scale as in Fig.8.The "looping" cycle in the plane is illustrated for six stars, where phase increases along the direction of the arrows.

Figure 10 .
Figure10.Program RRab star [Fe/H]s as estimated by our calibration, compared to their values based on highresolution spectroscopy.Each data point corresponds to one spectrum.The size of the dots indicate the absolute value of distance in phase from the maximum of the light curve (i.e., ϕ = {0.2,0.7} maps to |∆ϕ| = {0.2,0.3}.The dashed line is one-to-one or zero, and the solid line is the best fit.Red squares are average values taken across all the spectra for that star, and are the same as the retrieved values listed in the appendix in TableA3.

Figure 11 .
Figure 11.Histogram of estimated [Fe/H] of RR Lyrae stars in Stripe 82 and the Catalina footprint.

Figure 12 .Figure 13 .
Figure 12.Estimates of [Fe/H] of RR Lyraes from SDSS spectra, as compared with the SSPP pipeline, Liu et al. (2020), and Li et al. (2023).Left column: histograms; right column: binned S/N of the SDSS spectra upon which our estimates are based.S/N≥10 for all spectra, and the scale is truncated at 50.Dashed lines are one-to-one, and solid lines are best fits weighted by S/N.The SSPP values are from spectra acquired separately, or coadded into one, with no phase information.Liu et al. (2020) used synthetic spectral template matching to estimate [Fe/H].When more than one spectrum was available for a star they used observed spectra at phases unaffected by atmospheric shocks, based on the EW of the Ca II K line.Li et al. (2023) used photometric relations based on the metallicities of Liu et al. (2020).

Facilities
include the null set of extra coefficients, which is equivalent to the L94 relation.The validity of each model was tested with the Bayesian Information Criterion (BIC)(Schwarz et al. 1978) defined as: BIC ≡ ln(\) − 2 ln(L), (B4)

Table 1 .
Photometric Periods of Program Stars aFor compactness of notation, the uncertainties in parentheses correspond to the last two significant digits (e.g., 0.354311(17) means 0.354311±1.7×10−5 ).b Blazhko star (e.g., Skarka 2013) c Period is known to vary * The period of this star is poorly defined over long time baselines.The value tabulated here is based on photometry taken close in time to when spectra of the star were obtained at McDonald Observatory.
Observables of our program stars, as returned by rrlfe.Solid lines are RRab stars, and dashed lines are RRc stars.Qualitatively, the retrieved T eff , which is directly proportional to the Balmer EWs, resembles RR Lyrae photometric light curves.

Table A1
lists the literature sources of [Fe/H] based on high-resolution spectroscopy for our program stars, upon which we make a basis set to map our own estimated values.In addition, we list L94 upon which a consistent mapping between literature sources is based.TableA3lists both the mapped and estimated values for each star.

Table A1 .
Literature [Fe/H] Sources FeII Used spectra near light curve minimum.We use abundances which are the average of those of the Fe I and Fe II lines.See C17 Table3.

Table A2 continued
Table A2 (continued)The first error term is the average error in [Fe/H] across all spectra for that star.The second is the scatter in [Fe/H] between spectra of the same star.

Table A3 continued
TableA3(continued) The first error term is the average error in [Fe/H] across all spectra for that star.The second is the scatter in [Fe/H] between spectra of the same star.