Metallicity of Galactic RR Lyræ from Optical and Infrared Light Curves:
I. Period-Fourier-Metallicity Relations for Fundamental Mode RR Lyræ

As mentioned before, the sources selected for this analysis are chosen from both the high-resolution metallicity catalog of C20 and the full medium-resolution LAMOST DR6 and SDSS-SEGUE datasets, from which the spectrum selection criteria and $\Delta$S metallicity calibration of C20 have been applied. The resultant HR+$\Delta$S metallicity catalog is comprised of 8660 fundamental mode field RRLs that have also been cross-matched with the Gaia EDR3 database (Collaboration et al., 2020), as well as a number of other publicly available datasets, using an algorithm specifically developed for sparse catalogs (Marrese et al., 2019).

The HR+$\Delta$S metallicity catalog provides both a homogenized sample of RRab iron abundances gathered from various sources in literature (170 of which are derived from high-resolution spectra) and similarly homogeneous new $\Delta S$ metallicity estimates, as the $\Delta S$ calibration of C20 is based in part on the metallicity of 111 of the aforementioned HR RRab stars. For a complete and detailed description of the HR metallicity catalog’s demographics, the $\Delta S$ calibration, and the spectrum selection criterion, we refer the reader to the C20 paper.

All metallicities in this catalog are based upon the metallicity scale utilized by C20, which is based upon high-resolution spectroscopy, and the most updated iron line parameters, most of which have transition parameters derived in laboratory studies. This is the same scale utilized by For et al. (2011), Chadid et al. (2017), and Sneden et al. (2017), and will be the default scale used throughout this paper unless otherwise declared. Note, other literature field-RRL high-resolution works can be brought to the same metallicity scale with the addition of a simple offset. A full analysis of the offsets between various HR [Fe/H] scales is offered in C20. It is worth mentioning here that the often used Carretta et al. (2009, C09) [Fe/H] scale can be converted to this works scale with the addition of a small rigid shift of 0.08 dex.

In this work, we focus on the stars in the HR+$\Delta$S metallicity catalog that have a match in the ASAS-SN and WISE surveys. Figure 1 shows the distribution of period and metallicity for the subset of stars with available optical (ASAS-SN) or infrared (WISE) good quality light curves (i.e. passing the stringent photometric and Fourier decomposition criteria described in Section 3). The joint sample, also shown in the figure, comprises the smaller subset of RRLs with light curves available at both optical and infrared wavelengths. All samples cover the entire period range expected for RRab variables and are well representative of the metallicity of Galactic Halo RRLs, with a mean [Fe/H] abundance of $\approx-1.45$ in both the ASAS-SN and WISE samples. Figure 1 also shows that all histograms retain both the low metallicity $\textrm{[Fe/H]}\lesssim-2.2$ and the high metallicity $\textrm{[Fe/H]}\gtrsim-0.7$ tail present in the C20 catalog. For a discussion of the significance of populations in the Galactic Halo, we refer to Fabrizio et al. (2019)(hereafter F19). Here, we want to remark how the broad range in metallicity is an important feature of our samples, as it ensures broad leverage for accurate calibration of the metallicity slope in our period-$\phi_{31}$-[Fe/H] relation.

The optical ($V$-band) time series data for this analysis have been extracted from the ASAS-SN survey. The first telescope of ASAS-SN came online in 2013, and telescopes have gradually been added for a total of 24 telescopes scattered across the world as of early 2020. At its current capacity, ASAS-SN can survey the entire sky every night, providing high cadence $V$-band photometry, ideally suited to acquire long term, densely populated light curves of RRLs.

Out of the 6079 variables in the HR+$\Delta$S catalog with a match in ASAS-SN, we were able to extract good quality ASAS-SN light curves for 1980 stars using the procedure described in Section 3. Of these, 967 are in the joint sample, also having good quality light curves in the infrared. The left panel of Figure 2 shows the distribution of average $V$-band apparent magnitudes for both the entire optical dataset (solid black) and the stars in the joint sample (hatched). The cut-off at $V\simeq 17.4$ is due to the limiting magnitude of the ASAS-SN survey. The joint sample is instead truncated at $V\simeq 16$ due to the shallower photometric depth of the WISE survey. For a detailed description of the generation of the light curves from the ASAS-SN time-series, and calculation of average magnitudes for this dataset, we refer to Section 3.2.

We obtained archival infrared photometric time-series from the WISE mission $W1$ and $W2$ bands (3.4 and 4.6 $\micron$ respectively). The primary WISE mission surveyed the entire sky in four infrared bands every six months from January 2010 to the end of the post cryogenic phase in February 2011. The spacecraft was then reactivated in September 2013 in just the $W1$ and $W2$ bands as the NEOWISE mission. Still having sensitivity similar to the full cryogenic phase, NEOWISE, as of early 2020, has yielded 12 additional full-sky survey epochs (having a minimum of 12 measurements per survey epoch). With a minimum of $\sim$156 individual epochs ($\sim$12 from WISE and 144 from NEOWISE) for each position in the sky, collected over a baseline of almost ten years, the combined WISE and NEOWISE photometry provide the most comprehensive full-sky catalog in the infrared. From hereafter in the paper, the combination of photometry from both the primary WISE and ongoing NEOWISE mission will be referred for brevity as WISE.

Out of 6264 sources in the HR+$\Delta$S with a WISE catalog match, we were able to generate 1083 good quality light curves in at least one of the two adopted WISE bands (1083 and 707 stars in the $W1$ and $W2$ bands, respectively). Figure 2 (right panel) shows the WISE distribution of the average $W1$ apparent magnitude, for both the entire infrared dataset (solid grey) and those stars in the joint sample (hatched). Although the WISE sensitivity limit is $\sim 16$ mag in the two bands of interest, the figure shows a sharp drop at $W1\simeq 14.7$; fainter RRLs, while detected in the WISE and NEOWISE catalogs, tend to have noisy light curves that are then rejected by the quality control procedures described in Section 3. The histogram for the $W2$-band is similar, although the cut-off in the magnitude distribution happens at $W2\simeq 14.1$ due to the lower sensitivity in this band, resulting in noisier light curves. The magnitude range and other characteristics in the joint sample subset of WISE bands are nearly the same as the complete infrared dataset due to the large overlap with available ASAS-SN data.

Lastly, we select a dataset separate from those used in deriving our period-$\phi_{31}$-[Fe/H] relations to use as an independent check of this work (see Section 4.2). This dataset is comprised of eight globular clusters (GCs) homogeneously spread between [Fe/H] = $-1.1$ and $-2.3$ dex with a sizable number of RRab stars. These clusters are listed in Table 2. The GC dataset comes mainly from the homogeneous photometric database of P. B. Stetson¹¹1https://www.canfar.net/storage/list/STETSON/homogeneous/Latest_photometry_for_targets_with_at_least_BVI For additional information on the PBS database, readers may also contact Peter B. Stetson directly at [email protected] (hereafter, PBS), except for NGC 3201 data, which comes from Piersimoni et al. (2002). The data in the PBS database was collected from ground-based telescopes using archival data from 1984 to present. The main telescopes and cameras which contributed the most to the $V$-band data used in this work include the following: AAO LCOGT 1m (CCD), CTIO 0.9m (Tek2K), CTIO 1m (Y4KCam), CTIO 1.3m (ANDICAM), CTIO LCOGT 1m (CCD), Cerro Pachón SOAR 4.1m (SOI), La Silla NTT 3.6m (EMMI TK2048EB), La Silla ESO/MPI 2.2m (WFI), LCO Warsaw 1.3m (8k-MOSAIC), Maunakea CFHT 3.6m (CFH12K), ORM La Palma JKT 1m (EEV7), ORM La Palma INT 2.5m (WFC), SAAO LCOGT 1m (CCD). For further information about the summary of observing runs for each cluster, the bands observed, and the imagers and telescopes used, we direct the reader to the PBS database (previously mentioned).

The spectroscopic metallicities of these GCs (second column in Table 2) are listed in the scale of C09, and the third column of Table 2 shows the total number of fundamental mode RRLs in each GC according to Clement’s catalog²²2http://www.astro.utoronto.ca/~cclement/read.html (Clement et al., 2001). In parentheses, we list the number of RRab actually available in the PBS photometry. In addition, the average number of epochs per $V$-band light curve is shown in the fourth column.

The large temporal baseline of the ASAS-SN and WISE surveys makes our time-series very sensitive to even small errors in the nominal period of the stars. Even inaccuracies of $10^{-6}$ days, over the 10 years span of WISE and its extended mission, result in a $\sim 0.7$% phase shift when the photometric time series of a 0.55 day period of our typical RRab star is phased, readily detectable in our datasets. To avoid this issue, we have re-derived the periods of all our stars on the basis of their $V$ and $W1$ band³³3We did not use the $W2$-band for period determination due to its larger photometric error, often leading to less accurate periods. photometry, using the Lomb-Scargle method (Lomb, 1976; Scargle, 1982) within a search window from 0.2 to 1.5 days.

Typical re-derived periods differ from the ones in the literature by $10^{-5}$ days or less. Light curves with period discrepancy larger than $10^{-5}$ days between the $V$ and $W1$ datasets, or with the nominal period from literature, were inspected manually for further quality checks and possible rejection. Approximately 5% of derived periods fell in this category and required manual inspection. The inability to find a reliable period was also used as a rejection criterion for an individual light curve. Note that this quality assurance process causes the removal of some high amplitude Blazhko stars (Blažko, 1907), due to their naturally larger dispersion in the light curve causing a high “false alarm” probability in the period determination returned by the Lomb-Scargle routine. This removal is intended, as the photometric modulation of Blazhko stars can affect the Fourier decomposition of their light curves, complicating the dependence of [Fe/H] from period and $\phi_{31}$ (see e.g. Skarka 2014).

Rather than directly calculate the Fourier decomposition of the phased light curve, we elected to smooth it first to better deal with the uneven sampling of the observing epochs, and primarily to allow the removal of data points with excessive noise. The smoothing was performed using a Gaussian locally-weighted regression smoothing algorithm (GLOESS, Persson et al. 2004). This method places the phased light curve on an interpolated grid, to which a second-degree polynomial is locally fit to the photometric data. Fitting weights depend on both the photometric error and the Gaussian distance of the phased photometry from the interpolation point. The procedure was repeated multiple times in order to apply an iterative sigma clipping rejection scheme, designed for noisy data point removal. Only time series with a minimum of 30 valid photometric data points after sigma clipping were retained for processing, with the rest excluded from further consideration. A full description of our implementation of the GLOESS method is available in Neeley et al. (2015).

Figure 3 shows example light curves of the phased data (red points) for both a typical short and long period RRL star (DM Cyg and NSVS 13688631, respectively). Data points automatically rejected by the GLOESS fitting procedure are removed. From the GLOESS light curve, we have measured the mean magnitude (calculated as the mean of the smoothed light curve in flux units, converted back to magnitude), the amplitude, and the epoch of maximum (separately in each band). The mean magnitude uncertainty was calculated as the sum in quadrature of the photometric uncertainty and the uncertainty in the fit (see Neeley et al. 2015 for details). The error associated with the amplitude was instead defined as the standard deviation of the data points residuals with respect to the smoothed light curve. Stars with an amplitude equal to less than three times the amplitude error were found to correspond to excessively noisy light curves and were excluded from further analysis. Note that this process also excluded large modulation Blazhko stars, which would appear as noisy light curves with large uncertainty in amplitude.

Fourier decomposition is a widely adopted tool to quantify the shape of a light curve, as the lower order terms in the expansion are usually sufficient to fully characterize its shape (see e.g. Simon & Lee 1981). As justified above, for this work we elected to perform Fourier expansion of the evenly spaced smoothed light curves, rather than the individual data points, in the form:

We determined the Fourier coefficients with a weighted least-squares fit of Equation 1 to the smoothed GLOESS light curve and its locally calculated error. The locally calculated error was defined as the local photometric scatter around the smoothed light curve, estimated as the weighted standard deviation of the residual with the data convolved with the GLOESS smoothing kernel. We found that a fifth-order ($n=5$) Fourier expansion was sufficient to reproduce the shape of the light curve in each band. Examples are shown in Figure 3 for two typical RRLs. The Fourier decomposition fit (black solid line) is plotted on top of the actual photometric data in red.

Simon & Lee (1981) first showed that certain combinations of Fourier coefficients were directly related to some physical parameters of pulsating stars. These coefficients are typically defined either as linear combinations of Fourier phases:

In Appendix A, we discuss the correlation between various combinations of Fourier parameters among themselves, and with period, and how they help discriminate sources with different metallicities. Our analysis confirms the conclusions of early studies such as JK96, suggesting that a relation between period and $\phi_{31}$ is a good indicator of metallicity in the $V$-band. We found this to be true also in the infrared. Uncertainties in the calculated $\phi_{31}$ values are several orders of magnitude less than our best uncertainties in metallicity and are therefore deemed negligible throughout the rest of this work.

Figure 4 shows that the $W1$ and $W2$ bands produce indistinguishable values of $\phi_{31}$ due to the light curves in these bands being nearly identical. This is demonstrated quantitatively by Equation 4, which shows that the best-fit slope between $\phi_{31(W1)}$ and $\phi_{31(W2)}$ is within the errors close to unity, with a dispersion of 0.30.

We take advantage of this strong correlation by averaging, whenever possible, the $\phi_{31}$ parameters calculated for the two bands, falling back on $\phi_{31(W1)}$ or $\phi_{31(W2)}$ when only one is available. This essentially doubles the signal of the WISE light curves by combining the data of independent measurements taken at different wavelengths.

The metallicity of each star in either our ASAS-SN or WISE datasets can be effectively represented by a plane in the period, $\phi_{31}$, and [Fe/H] space, as is demonstrated in Appendix B by performing Principal Component Analysis (PCA) on each dataset. To determine the orientation of this plane, we adopted the Orthogonal Distance Regression (ODR) routine part of the SciPy package⁴⁴4https://docs.scipy.org/doc/scipy/reference/odr.html, which utilizes a modified trust-region Levenberg-Marquardt-type algorithm (Boggs & Rogers, 1990) to estimate the best fitting parameters. We chose ODR because it can be used similarly to PCA with both techniques minimizing the perpendicular distance to the fit with no differentiation between dependent and independent variables. This allows ODR to produce an optimal fit despite the correlation between variables we found in Appendix B (especially the well-known correlation between period and metallicity, see e.g. F19). These correlations, in the case of ordinary least-squares or similar fitting algorithms, would result in a trend of the fitted metallicity residuals with respect to the other variables.

We performed the ODR fit on both our ASAS-SN and WISE calibration sets using the equation:

Our best fit period-$\phi_{31}$-[Fe/H] relation based on the ASAS-SN $V$-band light curves is:

The period-$\phi_{31}$-[Fe/H] relation, based on the average $\phi_{31}$ parameters of the WISE light curves (averaged when possible between the $W1$ and $W2$ bands, as described in Section 3.3), is instead:

Parameter uncertainties have been checked with bootstrap re-sampling and are consistent with those via ODR. The errors from bootstrap re-sampling are those listed in the above relations. The Root Mean Square (RMS) of the two relations (0.41 and 0.50 dex for the ASAS-SN and WISE samples, respectively) are similar, showing that indeed accurate photometric metallicities can be obtained from infrared light curves. The RMS values are also comparable to the dispersion that we found with a non-parametric regression scheme, based on the $k$-NN method, of the same data (0.33 and 0.40 dex respectively, see Appendix C). This shows that our ODR fits provide an accurate description of the dependence of [Fe/H] from period and $\phi_{31}$, with uncertainty only slightly larger than the data’s intrinsic scatter. Table 3 shows the derived photometric properties for both the $V$-band and mid-IR (WISE) calibration datasets, including the period, $\phi_{31}$ value, and photometric metallicity in each band. Following the light curve fitting, quality control process, and the final plane fit described above, we were left with 1980 variables with good quality ASAS-SN light curves, and 1083 variables with good WISE (in at least one of the two $W1$ and $W2$ bands) light curves.

Figures 5 and 6 (for the ASAS-SN and WISE samples, respectively) demonstrate in graphical form the ability of our fits to provide photometric metallicities in agreement with the spectroscopic values in the HR+$\Delta$S calibration sample. The top left panel shows the distribution of the sources in the $\phi_{31}$ vs. period plane, superimposed with the relations in Equation 6 and 7 calculated for fixed values of [Fe/H]. Both data and fit lines are color-binned by metallicity. Note the good agreement between the distribution of the spectroscopic metallicities of individual sources with the locus corresponding to the same metallicity bin defined by the spectroscopic relations. The histogram in the bottom left panel shows how the distribution of the photometric [Fe/H] from our fits closely reproduces the metallicity distribution, from spectroscopy, of the calibration sample. The two plots on the right instead show the distribution of the residuals between the photometric and spectroscopic metallicities: the top panel confirms that our choice of fitting the data with an ODR method indeed avoids residual trends, while the bottom panel shows a close-to-Gaussian distribution of the residuals.

The analysis in Section 3.4 shows that optical and infrared period-$\phi_{31}$-[Fe/H] relations provide photometric metallicities of comparable accuracy. However, we still need to test if Equation 6 and 7 provide consistent values of [Fe/H] for individual stars. Figure 7 shows that this is indeed the case: we find an excellent agreement between the photometric metallicities derived at the two wavelengths ranges for the joint sample on a per-star basis. The dispersion between the two datasets, of 0.44 dex, is comparable with the RMS of the individual fits as well as with the dispersion of the residuals from the $k$-NN method in Appendix C.

This tight correlation between optical and infrared photometric metallicities, with no apparent trends, led us to consider whether averaging the two photometric metallicities would yield a value significantly closer to the spectroscopic metallicity from the HR+$\Delta$S sample. Furthermore, besides the obvious advantage of increasing the statistics of the photometric measurements, by combining optical and infrared data we can probe the effects that different aspects of stellar atmospheres have on the RRL light curves acquired in these two separate wavelength ranges (the optical light curves are dominated by temperature, while infrared emission follows more closely radius variations).

Based on these arguments, we average together the photometric metallicities derived from the optical and infrared relationships for each star in the joint sample. Figure 8 shows that the average photometric metallicity (black histogram in left panel) closely mirrors the spectroscopic metallicity values from the HR+$\Delta$S calibration sample (grey histogram) for the entire metallicity range. The right panel shows that the residuals between the average photometric [Fe/H] and the spectroscopic metallicities fall fairly symmetrically around zero.

The RMS dispersion between the two sets of values is $\sim$0.37 dex: smaller than the individual error in the optical ($\sim$0.41 dex) or infrared ($\sim$0.50 dex) sample alone. This RMS value is near the dispersion we found with the k-NN method in Appendix C (0.33 and 0.40 dex for the ASAS-SN and WISE datasets, respectively), which shows that combining these datasets allows one to approach the local scatter in the individual datasets.

In order to test the $V$-band relation obtained in Section 3.3 on an independent sample, we selected a list of eight Galactic GCs homogeneously spread between [Fe/H] = $-1.1$ and $-2.3$ dex. The selected GCs are described in Section 2.4 and listed in Table 2. For each cluster, we chose the RRab stars with the best sampled light curves in the $V$-band, avoiding those RRab stars that suffer the Blazhko effect because of the modulation of the amplitude and the shape of their light curves. A Fourier decomposition was performed on each light curve to obtain their $\phi_{31}$ parameter. The period-$\phi_{31}$-[Fe/H] relation (Equation 6) was then applied to estimate the metallicity of each cluster star.

Figure 9 shows the spectroscopic [Fe/H] versus the mean photometric [Fe/H] values calculated for each GC with our relation. For consistency with the field RRLs described in Section 2.1, the spectroscopic metallicities of these GCs (second column in Table 2) from C09 were converted into the scale of this paper with the addition of an offset of $0.08$ dex, as noted in C20 and described in Section 2.1. The error bars correspond to the spectroscopic and photometric metallicity uncertainties. The former comes from Carretta et al. (2009), while the latter was assessed as the standard error of the mean for star photometric metallicities in each cluster. The good performance obtained using the period-$\phi_{31}$-[Fe/H] relation derived in this work is clearly noticeable in Figure 9. There are no signals of possible systematic effects, and the predicted [Fe/H] using the period-$\phi_{31}$-[Fe/H] relation obtained in this work is within $\pm$0.09 dex of the spectroscopic [Fe/H]. In fact, this accuracy is similar to the metallicity uncertainties measured from the high-resolution spectroscopy on individual stars.

In this section, we compare our optical period-$\phi_{31}$-[Fe/H] relation (Equation 6) with previous relations found in the literature for similar wavelength ranges. In particular, we focus on the relations found in JK96, N13, MV16, and IB20. Note that a similar comparison for our mid-infrared relation (Equation 7) is not possible, since we could not find any previous work studying the relation between metallicity and Fourier parameters at wavelengths longer than the $I$-band (Smolec, 2005). The results are shown in Figure 10, which plots the [Fe/H] abundances calculated with the periods and $\phi_{31}$ parameters from the ASAS-SN sample, using our relation vs. the relations published in the literature.

The top left panel shows the comparison with the formula found in JK96 (their Equation 3). Their relation was derived using a total of 81 field RRab with $V$-band photometry from heterogeneous observations at various sites. Due to lack of phase coverage or excessive noise, a direct Fourier fit in JK96 was often not available, and individual polynomial fits or small parabolas were fit to light curve segments as needed. Since JK96 adopted metallicities based on the high-dispersion spectroscopy scale of Jurcsik (1995), for consistency, we first converted the metallicities derived with their relation to the C09 scale, using the relation from Kapakos et al. (2011): $\rm{[Fe/H]_{C09}}=1.001\rm{[Fe/H]_{JK96}}-0.112$. A scale offset of $0.08$ was then added to transform the C09 scale to that adopted by this work’s HR+$\Delta$S calibration sample (see Section 2.1). We found agreement within the calibration range of JK96 (red horizontal lines), with an RMS value of 0.17 dex, which is smaller than the nominal uncertainty in the calculated [Fe/H] abundances, and only a mildly apparent trend.

N13 introduced a quadratic period-$\phi_{31}$-[Fe/H] relation (their Equation 2), calibrated using stars observed in the Kepler photometric band during the first 970 days of the Kepler mission. Since Kepler, during its nominal mission, surveyed the same region of the sky nearly continuously, each star resulted in having $\sim$350,000 data points spread over $\sim$2,500 pulsation cycles, yielding the best sampled light curves (and Fourier decomposition) currently available for any sample of RRLs. However, due to the fixed field of view and shallow depth of Kepler’s field, their calibration dataset only had a total of 26 RRab stars with accurate metallicity measurements, nine of which were Blazhko. In order to use the relation of N13, derived in the Kepler photometric system ($Kp$), with our data, we had to convert the $V$-band $\phi_{31}$ value to the Kepler system by adding the systematic offset derived by Nemec et al. (2011):

Furthermore, we converted the [Fe/H] abundances provided by N13 from their adopted metallicity scale C09 to this work scale, using the previously noted scale offset of 0.08 dex. The results are presented in the top right panel of Figure 10 and show an excellent agreement for higher metallicities in the N13 calibration range ($-1.5\lesssim\textrm{[Fe/H]}\leq 0.03$ dex). For lower metallicities, however, the two relations rapidly diverge (RMS $\simeq 0.30$ dex if calculated over their entire calibration range), caused by the higher-order term in N13, and possibly due to the scarcity of calibrators in their samples for low [Fe/H] (they only have one RRL with $\textrm{[Fe/H]}\lesssim-2.0$ dex.

Next, we compare our relationship with the one derived by MV16 (end of Section 2 in their paper), calibrated using a sample of 381 RRLs in globular clusters binned by period, with the addition of 8 field RRLs chosen to extend the metallicity range of the sample. The C09 metallicity scale used by MV16 has also been converted to that of this work to allow a comparison with our sample. The result is presented in the bottom-right panel of Figure 10, showing a general agreement with our fit within the calibration range of MV16, albeit with a small slope with respect to our relation. Note that we set the lower limit of the calibration range plotted to $\textrm{[Fe/H]}\simeq-0.99$ since there was a discontinuity in the data beyond which only had two variable stars.

Finally, IB20 introduced a $G$-band period-$\phi_{31}$-[Fe/H] relation (their equation 3), based on Gaia DR2 light curves. The relation was calibrated with a sample of 84 stars with known spectroscopic metallicity. To allow comparison with our metallicities, we had first to convert the $V$-band $\phi_{31}$ value of our dataset to the $G$-band system, using the relation derived by Clementini et al. 2016:

Furthermore, an additional $\pi$ offset had to be subtracted from $\phi_{31}$ to set the coefficient on the same scale as IB20. Metallicity abundances from IB20 were on the scale of Zinn & West (1984, ZW84) and were thus converted using the following equation from C09 with a subsequent scale conversion to that of this work: $\rm{[Fe/H]_{C09}}=1.105\rm{[Fe/H]_{ZW84}}+0.160$.

The results are shown in the bottom right panel of Figure 10 and exhibit a clear trend between the two relations for the entire range of metallicity. It should be noted that although the stars used for the calibration of the IB20 relation range from ($-2.53\leq\textrm{[Fe/H]}\leq 0.33$), the fit lacks a significant number of calibrators at the low and high metallicity ends. In comparison to this work, their relation tends to overestimate the metallicity at the metal-poor end and underestimate the metallicity at the metal-rich end.

To assess the reliability of the period-$\phi_{31}$-metallicity relations, both derived in this paper and from literature, in providing [Fe/H] abundances, we compare their predictions with metallicities measured from high-resolution spectroscopy in the HR+$\Delta$S sample. The results are shown in Figure 11, plotting the difference between spectroscopic and photometric [Fe/H] abundance for each of the stars in C20 with metallicity from a high-resolution spectrum. The RMS scatter of each relation, calculated with respect to zero residuals over the entire spectroscopic metallicity range, is indicated for each photometric relation.

The figure shows that all relations have similar residuals ($\textrm{RMS}\simeq 0.3$-$0.4$ dex). The scatter in our relations (top row), in particular, are consistent with the expected uncertainties of our photometric metallicities. There is marginal evidence for a systematic shift downwards in panels of Figure 11 that we attribute to the minor differences between the $\Delta$S and HR metallicities as shown in Figure 9 of C20. All relations from literature show a trend of the residuals with high-resolution spectroscopic metallicity due to their tendency to overestimate [Fe/H] abundances at the low-metallicity range and underestimate high-metallicities, with the exception of N13 that shows a significantly larger scatter ($\textrm{RMS}\simeq 0.57$ dex). This apparent trend may result from an inadequate sampling of the metal-rich and metal-poor tails or possibly non-LTE effects in each of the calibration datasets’ metallicities that are not fully taken into account. More investigation is needed in order to clarify the mechanisms behind this trend. This trend is largest in the IB20, which underestimates the metallicity of all stars with $\textrm{[Fe/H]}\lesssim-2.0$ and overestimates the metallicity for $\textrm{[Fe/H]}\gtrsim 0.7$ dex, which was already partially noted by IB20 within their calibration range.

It is worth noting that the RMS values presented in Figure 11 are generally larger than those quoted inside these respective works of literature. The RMS values quoted in other papers are a measure of their fit’s dispersion with respect to their calibration dataset and should not be construed as a measurement of accuracy between two relations. A comparison between works needs to be made on an identical validation dataset (as was done in this section). We attribute the larger dispersion in our relations as denoting a more realistic representation of the intrinsic scatter found in the period-$\phi_{31}$-[Fe/H] relation, which is not apparent in smaller calibration samples, usually covering a smaller metallicity range. Our calibration sample is 76$\times$ larger than N13 (26 RRab) and 24$\times$ larger than JK96 (81 RRab). We could decrease the dispersion of this works’ relations by artificially cutting the tails of the metallicity distribution (very metal-poor/metal-rich). However, this approach would introduce severe systematic drifts in the metallicity estimates — i.e. the same problem affecting other current calibrations — which is what we have strived to resolve in this work.