Metallicity of Galactic RR Lyrae from Optical and Infrared Light Curves. I. Period–Fourier–Metallicity Relations for Fundamental-mode RR Lyrae

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

The HR+ΔS metallicity catalog provides both a homogenized sample of RRab iron abundances gathered from various sources in literature (170 of which are derived from HR spectra) and similarly homogeneous new ΔS metallicity estimates, as the ΔS calibration of C21 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 ΔS calibration, and the spectrum selection criterion, we refer the reader to the C21 paper.

All metallicities in this catalog are based on the metallicity scale utilized by C21, which is based on HR 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 that other literature field RRL HR 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 C21. It is worth mentioning here that the often-used (Carretta et al. 2009, hereafter C09) [Fe/H] scale can be converted to this work's scale with the addition of a small rigid shift of 0.08 dex.

In this work, we focus on the stars in the HR+Δ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 ≈−1.45 in both the ASAS-SN and WISE samples. Figure 1 also shows that all histograms retain both the low-metallicity [Fe/H] ≲ −2.2 and the high-metallicity [Fe/H] ≳ −0.7 tail present in the C21 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 that 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–ϕ₃₁–[Fe/H] relation.

Zoom In Zoom Out Reset image size

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+Δ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 data set (solid black) and the stars in the joint sample (hatched). The cutoff at V ≃ 17.4 is due to the limiting magnitude of the ASAS-SN survey. The joint sample is instead truncated at V ≃ 16 owing 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 data set, we refer to Section 3.2.

Zoom In Zoom Out Reset image size

We obtained archival infrared photometric time series from the WISE mission W1 and W2 bands (3.4 and 4.6 μm, respectively). The primary WISE mission surveyed the entire sky in four infrared bands every 6 months from 2010 January to the end of the post-cryogenic phase in 2011 February. The spacecraft was then reactivated in 2013 September 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 ∼156 individual epochs (∼12 from WISE and 144 from NEOWISE) for each position in the sky, collected over a baseline of almost 10 yr, the combined WISE and NEOWISE photometry provides the most comprehensive full-sky catalog in the infrared. Hereafter in the paper, the combination of photometry from both the primary WISE and ongoing NEOWISE mission will be referred to for brevity as WISE.

Out of 6264 sources in the HR+Δ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 data set (solid gray) and those stars in the joint sample (hatched). Although the WISE sensitivity limit is ∼16 mag in the two bands of interest, the figure shows a sharp drop at W1 ≃ 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 cutoff in the magnitude distribution happens at W2 ≃ 14.1 owing 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 data set owing to the large overlap with available ASAS-SN data.

Lastly, we select a data set separate from those used in deriving our period–ϕ₃₁–[Fe/H] relations to use as an independent check of this work (see Section 4.2). This data set is composed of eight 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 data set comes mainly from the homogeneous photometric database of P. B. Stetson ¹⁴ (hereafter PBS), except for NGC 3201 data, which comes from Piersimoni et al. (2002). The data in the PBS database were collected from ground-based telescopes using archival data from 1984 to the present. The main telescopes and cameras that contributed the most to the V-band data used in this work include the following: AAO LCOGT 1 m (CCD), CTIO 0.9 m (Tek2K), CTIO 1 m (Y4KCam), CTIO 1.3 m (ANDICAM), CTIO LCOGT 1 m (CCD), Cerro Pachón SOAR 4.1 m (SOI), La Silla NTT 3.6 m (EMMI TK2048EB), La Silla ESO/MPI 2.2 m (WFI), LCO Warsaw 1.3 m (8k-MOSAIC), Maunakea CFHT 3.6 m (CFH12K), ORM La Palma JKT 1 m (EEV7), ORM La Palma INT 2.5 m (WFC), SAAO LCOGT 1 m (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 ¹⁵ (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⁻⁶ days, over the 10 yr span of WISE and its extended mission, result in a ∼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 data sets. To avoid this issue, we have rederived the periods of all our stars on the basis of their V- and W1-band ¹⁶ photometry, using the Lomb–Scargle method (Lomb 1976; Scargle 1982) within a search window from 0.2 to 1.5 days.

Typical rederived periods differ from the ones in the literature by 10⁻⁵ days or less. Light curves with period discrepancy larger than 10⁻⁵ days between the V and W1 data sets, or with the nominal period from the 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 ϕ₃₁ (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' residual 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

$$\begin{eqnarray}&&m({\rm{\Phi }})={A}_{0}+\sum _{i=1}^{n}{A}_{i}\sin [2\pi i({\rm{\Phi }}+{{\rm{\Phi }}}_{0})+{\phi }_{i}],\end{eqnarray} \tag{ 1 }$$

where m(Φ) is the observed magnitude for either the ASAS-SN or WISE bands, A₀ is the mean magnitude, n is the order of the expansion, Φ is the phase from the GLOESS light curve varying from 0 to 1, Φ₀ is the phase that corresponds to the time of maximum light T₀, and A_i and ϕ_i are the ith-order Fourier amplitude and phase coefficients, respectively.

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,

$$\begin{eqnarray}&&{\phi }_{{ij}}=j\cdot {\phi }_{i}-i\cdot {\phi }_{j},\end{eqnarray} \tag{ 2 }$$

where ϕ_ij is cyclic in nature and ranges from 0 to 2π, or as ratios of the Fourier amplitudes,

$$\begin{eqnarray}&&{R}_{{ij}}=\displaystyle \frac{{A}_{i}}{{A}_{j}}.\end{eqnarray} \tag{ 3 }$$

In Appendix A, we discuss the correlation between various combinations of Fourier parameters among themselves, as well as 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 ϕ₃₁ is a good indicator of metallicity in the V band. We found this to be true also in the infrared. Uncertainties in the calculated ϕ₃₁ 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 ϕ₃₁ owing 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 ϕ_31(W1) and ϕ_31(W2) is within the errors close to unity, with a dispersion of 0.30:

$$\begin{eqnarray}&&{\phi }_{31(W2)}=(1.011\pm 0.023)\cdot {\phi }_{31(W1)}+(0.066\pm 0.043).\end{eqnarray} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

We take advantage of this strong correlation by averaging, whenever possible, the ϕ₃₁ parameters calculated for the two bands, falling back on ϕ_31(W1) or ϕ_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 data sets can be effectively represented by a plane in the period, ϕ₃₁, and [Fe/H] space, as is demonstrated in Appendix B by performing principal component analysis (PCA) on each data set. To determine the orientation of this plane, we adopted the orthogonal distance regression (ODR) routine part of the SciPy package, ¹⁷ 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

$$\begin{eqnarray}&&[\mathrm{Fe}/{\rm{H}}]=a+b\cdot (P-{P}_{0})+c\cdot ({\phi }_{31}-{\phi }_{{31}_{0}}),\end{eqnarray} \tag{ 5 }$$

where P₀ and ${\phi }_{{31}_{0}}$ are pivot offsets necessary to add an extra element of robustness in the fitting procedure and reduce the fitting parameter uncertainties. We chose the pivot offsets to be near the mean period and ϕ₃₁ value of each data set, equal to P₀ = 0.58 days for both data sets, and ${\phi }_{{31}_{0}}=5.25$ and 1.90 rad for ASAS-SN and WISE, respectively. To ensure an accurate fit, two rounds of fitting were performed with an intermediate 4σ clipping between the fitted and calibration [Fe/H] (removing ≲1% of the stars in each calibration data set).

Our best-fit period–ϕ₃₁–[Fe/H] relation based on the ASAS-SN V-band light curves is

$$\begin{eqnarray}\begin{array}{rcl}[\mathrm{Fe}/{\rm{H}}] & = & (-1.22\pm 0.01)+(-7.60\pm 0.24)\cdot (P-0.58)\\ & & +\ (1.42\pm 0.05)\cdot ({\phi }_{31}-5.25)\ ;\ \mathrm{rms}=0.41,\end{array}\end{eqnarray} \tag{ 6 }$$

where, due to the 2π periodic ambiguity in the ϕ₃₁ coefficient, some ϕ₃₁ values required adding 2π to their phase in order to lie closer to the mean ϕ₃₁ value, as suggested by JK96.

The period–ϕ₃₁–[Fe/H] relation, based on the average ϕ₃₁ parameters of the WISE light curves (averaged when possible between the W1 and W2 bands, as described in Section 3.3), is instead

$$\begin{eqnarray}\begin{array}{rcl}[\mathrm{Fe}/{\rm{H}}] & = & (-1.47\pm 0.02)+(-8.33\pm 0.34)\cdot (P-0.58)\\ & & +\ (0.92\pm 0.05)\cdot ({\phi }_{31}-1.90)\ ;\ \mathrm{rms}=0.50.\end{array}\end{eqnarray} \tag{ 7 }$$

Parameter uncertainties have been checked with bootstrap resampling and are consistent with those via ODR. The errors from bootstrap resampling are those listed in the above relations. The rms values 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 nonparametric 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] on period and ϕ₃₁, 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 data sets, including the period, ϕ₃₁ 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+ΔS calibration sample. The top left panel shows the distribution of the sources in the ϕ₃₁ versus period plane, superimposed with the relations in Equations (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 Note the good agreement between the distribution of the spectroscopic metallicities of individual sources with the locus of the photometric relations corresponding to the same metallicity bin. 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 panels 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The analysis in Section 3.4 shows that optical and infrared period–ϕ₃₁–[Fe/H] relations provide photometric metallicities of comparable accuracy. However, we still need to test whether Equations (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 wavelength ranges for the joint sample on a per-star basis. The dispersion between the two data sets, of 0.44 dex, is comparable to the rms of the individual fits, as well as to the dispersion of the residuals from the k-NN method in Appendix C.

Zoom In Zoom Out Reset image size

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+Δ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+ΔS calibration sample (gray 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.

Zoom In Zoom Out Reset image size

The rms dispersion between the two sets of values is ∼0.37 dex—smaller than the individual error in the optical (∼0.41 dex) or infrared (∼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 data sets, respectively), which shows that combining these data sets allows one to approach the local scatter in the individual data sets.

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 ϕ₃₁ parameter. The period–ϕ₃₁–[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 C21 and described in Section 2.1. The error bars correspond to the spectroscopic and photometric metallicity uncertainties. The former comes from C09, 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–ϕ₃₁–[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–ϕ₃₁–[Fe/H] relation obtained in this work is within ± 0.09 dex of the spectroscopic [Fe/H]. In fact, this accuracy is similar to the metallicity uncertainties measured from the HR spectroscopy on individual stars.

Zoom In Zoom Out Reset image size

In this section, we compare our optical period–ϕ₃₁–[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 IB21. 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 ϕ₃₁ parameters from the ASAS-SN sample, using our relation versus the relations published in the literature.

Zoom In Zoom Out Reset image size

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 stars with V-band photometry from heterogeneous observations at various sites. Due to the 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): [Fe/H]_C09 = 1.001[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+Δ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–ϕ₃₁–[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 ∼350,000 data points spread over ∼2500 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 data set 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 ϕ₃₁ value to the Kepler system by adding the systematic offset derived by Nemec et al. (2011):

$$\begin{eqnarray}&&{\phi }_{31}(V)={\phi }_{31}({Kp})-(0.151\pm 0.026).\end{eqnarray} \tag{ 8 }$$

Furthermore, we converted the [Fe/H] abundances provided by N13 from their adopted metallicity scale C09 to the scale in this work, 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 dex ≲ [Fe/H] ≤ 0.03 dex). For lower metallicities, however, the two relations rapidly diverge (rms ≃ 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 [Fe/H] ≲ −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 GCs binned by period, with the addition of eight 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 [Fe/H] ≃ −0.99 since there was a discontinuity in the data beyond which only had two variable stars.

Finally, IB21 introduced a G-band period–ϕ₃₁–[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 first had to convert the V-band ϕ₃₁ value of our data set to the G-band system, using the relation derived by Clementini et al. (2016):

$$\begin{eqnarray}\begin{array}{rcl}{\phi }_{31}(G) & = & (0.104\pm 0.020)+(1.000\pm 0.008){\phi }_{31}(V);\\ \mathrm{rms} & = & 0.055.\end{array}\end{eqnarray} \tag{ 9 }$$

Furthermore, an additional π offset had to be subtracted from ϕ₃₁ to set the coefficient on the same scale as IB21. Metallicity abundances from IB21 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: [Fe/H]_C09 = 1.105[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 IB21 relation are in the range of −2.53 ≤ [Fe/H] ≤ 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–ϕ₃₁–metallicity relations, both derived in this paper and from the literature, in providing [Fe/H] abundances, we compare their predictions with metallicities measured from HR spectroscopy in the HR+Δ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 C21 with metallicity from an HR 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.

Zoom In Zoom Out Reset image size

The figure shows that all relations have similar residuals (rms ≃ 0.3–0.4 dex). The scatter in our relations (top row), in particular, is consistent with the expected uncertainties of our photometric metallicities. There is marginal evidence for a systematic shift downward in panels of Figure 11 that we attribute to the minor differences between the ΔS and HR metallicities as shown in Figure 9 of C21. All relations from the literature show a trend of the residuals with HR 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 showing a significantly larger scatter (rms ≃ 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 data set's 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 IB21, which underestimates the metallicity of all stars with [Fe/H] ≲ −2.0 and overestimates the metallicity for [Fe/H] ≳ 0.7 dex, which was already partially noted by IB21 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 data set 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 data set (as was done in this section). We consider the larger dispersion in our relations as denoting a more realistic representation of the intrinsic scatter found in the period–ϕ₃₁–[Fe/H] relation, which is not apparent in smaller calibration samples, usually covering a smaller metallicity range. Our calibration sample is 76× larger than N13 (26 RRab) and 24× larger than JK96 (81 RRab). We could decrease the dispersion of this work's 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.