A Data-driven Study of RR Lyrae Near-IR Light Curves: Principal Component Analysis, Robust Fits, and Metallicity Estimates

We have collected high-quality K_S and J-band photometry available from the literature for 101 RRLs. Table 1 summarizes the training set. The available data can be categorized into three main types: near-IR photometry taken with the aim of Baade-Wesselink analysis of RRLs (Barnes et al. 1992; Cacciari et al. 1992; Fernley et al. 1989, 1990; Jones et al. 1987a, 1987b, 1988, 1992; Liu & Janes 1989; Skillen et al. 1993); serendipitous observations of RRLs in the 2MASS (Szabó et al. 2014) and WFCAM (Ferreira Lopes et al. 2015) calibration fields; and the extensive J and K_S variability study of ω Centauri by Navarrete et al. (2015). There are other sources of near-IR time-series photometry available in the literature for RRLs (see, e.g., Angeloni et al. 2014), but we decided to utilize only photometry where the phase coverage and quality of the observations are adequate for the accurate determination of the near-IR light-curve shapes, at least in the K_S band.

Examining Table 1 reveals that the variables in the field of ω Centauri (Navarrete et al. 2015) dominate the sample. Although this could introduce a heavy bias toward a particular metallicity, ω Cen contains multiple stellar populations (e.g., Valcarce & Catelan 2011; Gratton et al. 2012, and references therein), with two contributing to the RRL sample (Sollima et al. 2006, and references therein). These two populations have [Fe/H] ∼ −1.2 and ∼−1.7. The variables utilized from the photometry of Ferreira Lopes et al. (2015) are all members of the globular cluster M3, with a metallicity of ∼−1.5 (Carretta et al. 2009). Furthermore, the field RRLs have a wider metallicity distribution, allowing us to cover the physical parameter space of RRLs.

As for the period distribution, the sample covers most of the possible period range of RRLs, from 0.39 to 1.01 days. While very metal-rich RRLs can have periods as short as about 0.35 days, their optical light-curve shapes are not drastically different from those with periods around 0.4 days. Therefore, we can surmise that the present data set can be considered representative of RRLs and their light-curve shapes.

To apply PCA to our data outlined in Section 2.1, we have to describe their phase-folded light-curve shape. As mentioned in Section 1, it can be hard to accurately represent the near-IR RRL light curves as a Fourier series. Furthermore, our light curves have vastly different numbers of data points: NSV 660 has almost 3000, while for the RRL found in ω Centauri, Navarrete et al. (2015) only acquired a total of 100 and 42 epochs in the K_S and J bands, respectively.

As an alternative, folded light curves can be described as the linear sum of a series of different basis functions of choice, such as Gaussians, which can be aligned to the phased light-curve points with ordinary least squares (OLS) regression. We have chosen to adopt the Gaussian sum fitting example of Figure 8.4 of Ivezić et al. (2014), to the K_S and J-band RRL light curves. As RRL light curves are strictly periodic, we have decided to replace the Gaussians with their periodic analog, the circular normal (also called Von Mises) distributions of the form

$$\begin{eqnarray}&&f(x)=\displaystyle \frac{{e}^{\kappa \cos (x-\mu )}}{2\pi {I}_{0}(\kappa )},\end{eqnarray} \tag{ 1 }$$

where μ is the measure of location (analogous to the mean of the Gaussian distribution), κ is the measure of concentration (where 1/κ is analogous to the variance, σ² of a Gaussian), and ${I}_{0}(\kappa )$ is the modified Bessel function of order 0. To model the light-curve shapes, we define the sum of 100 of these basis functions, distributed evenly between phases 0 and 1:

$$\begin{eqnarray}&&\mathrm{LC}=m+\displaystyle \sum _{i=0}^{99}{A}_{i}\displaystyle \frac{{e}^{\kappa \cos [2\pi (x-\tfrac{i}{100})]}}{2\pi {I}_{0}(\kappa )},\end{eqnarray} \tag{ 2 }$$

where A_i are the individual amplitudes of the circular normal distributions, and m is the intercept of the fit. Although we could use OLS to find the amplitudes of Equation (2), most of the light curves have less than 100 points available, leading to an underdetermined problem. In such cases, regularization can be introduced to penalize the magnitude of independent parameters (in our case, the amplitudes A_i), by modifying the loss function. We do so by utilizing the least absolute shrinkage and selection operator (LASSO or L1 regularization; Tibshirani 1996), which adds the sum of absolute values of the regression coefficients multiplied by a regularization parameter α to the loss function. We note that by utilizing LASSO, generally most coefficients end up being zero.⁶

We have utilized the linear regression routines of scikit-learn (Pedregosa et al. 2011) to fit the folded K_S and J light curves of each RRL with Equation (2) utilizing LASSO regularization. Our fit has two hyperparameters: κ and α. We have utilized both cross validation (leave-one-out for stars with few light-curve points, N-fold otherwise) and manual inspection of the resulting light-curve fits to determine the optimal values of these parameters. High values for the concentration parameter κ grants our model the ability to fit sharp features, such as the rising branches of certain variables, but can cause an overfit in phase ranges with few points. Conversely, at low values of κ the model cannot fit sharp features. We have found that a numerical value of κ = 6 is optimal for our model for both the K_S and J-band light curves of RRLs. As for the regularization parameter α, our cross validation resulted in different optimal values for different stars, typically in the range between 10⁻⁴ and 10⁻⁵. As visual inspection did not reveal significant differences when changing the regularization parameter between these two values, we have adopted an intermediate value of 10^−4.5 for all of our stars.

Figure 1 illustrates the quality of the light curves and their fits. During our fitting process, some outlying points have been removed manually and the periods of some variables have been revised, when tension existed between different photometric sources. Table 1 contains the revised periods for all variables. Figure 2 compares our fit with Fourier series of different orders. As can be seen, our method provides a better representation for variables with light-curve gaps.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The utilized light curves had been obtained in a variety of photometric systems, as detailed in Table 1. Therefore, the resulting light-curve fits were all transformed to the photometric system of VISTA. For variables not in the VISTA, WFCAM, or 2MASS systems, first they were transformed to the system of 2MASS utilizing the updated transformation formulae⁷ of Carpenter (2001). Then, the 2MASS and WFCAM photometry was transformed to the VISTA system with the help of the CASUVERS 1.4 transformations given by the Cambridge Astronomical Survey Unit (CASU).⁸

To apply PCA, we sample our light-curve fits on a grid of 100 phase points, evenly distributed phases from 0.0 to 0.99. In the analysis of pulsating variables, it is customary to set the light-curve maxima at phase 0.0; however, inspecting the K_S-band light curve examples of Figure 1 reveals that the maxima of RRLs in the K_S band are ill defined: the timing of the maxima of the light curve depends heavily on the strength of the bump on the rising branch. Therefore, we have chosen to align our light curves by the much sharper feature of the light-curve minima (similar to the case of eclipsing binaries). Furthermore, in PCA, the sample values (the magnitudes in our case) are usually normalized to a mean of 0 and scatter of 1 along each dimension (phase). As our goal is to describe the light-curve shapes of RRLs in the K_S band as a linear combination of PCs, we have chosen to normalize each light curve independently to have a mean of 0 and scatter of 1. These aligned, normalized input light curves can be seen in Figure 3.

Zoom In Zoom Out Reset image size

We carry out PCA by utilizing singular value decomposition, adopted from the PCA module of scikit-learn (Pedregosa et al. 2011). Figure 4 shows the first six PCs, according to our decomposition. As we have chosen not to normalize in each phase point, the first PC contains the average light-curve shape of the normalized light curves. Including further components to describe the light curves modifies this average shape, and this can be easily understood in the context of the individual light curves, for PCs of low order: the second component can make the bump at the end of the rising branch (around phase 0.15) more or less pronounced, while the third component is important to reproduce the double-peaked light curves displayed by some of the variables.

Zoom In Zoom Out Reset image size

The power of the PCA lies in the fact that generally the first few PCs, together with their amplitudes, are sufficient to describe the original input data, and the rest of the PCs can be discarded. The number of significant PCs can be decided by examining the fraction of the variance explained by the components. As is obvious, the first PC dominates the explained variance, because we normalized each light curve individually instead of normalizing the magnitudes along each dimension (phase). If we look at the variance explained by the rest of the PCs, they explain 3.68, 1.39, 0.96, 0.72, etc. percent of the total variance, or 88.8, 3.4, 2.3, 1.7, etc. percent of the residual variance, when the variance explained by the first PC is subtracted from the total variance. On the basis of these values, we deem the first four PCs are sufficient for describing the K_S-band light-curve shapes of RRab stars.

We applied the PCA method on the normalized light curves of our sample of variables to emphasize the light-curve shape differences among RRLs in the K_S band, instead of the different pulsation amplitudes of each variable. Consequently, the individual amplitudes (also called eigenvalues) of the first four PCs, u_1j, u_2j, u_3j and u_4j, where j = 1..101 is the index of the variables, carry no amplitude information on the original light curves. By multiplying these amplitudes with the normalization constants used to normalize each of the light curves, or by utilizing OLS to directly align the PCs to the light curves, we can determine the PC amplitudes U_1j, U_2j, U_3j and U_4j for each star in our sample. The sum of the PCs multiplied with these amplitudes recover the original light-curve shapes (and amplitudes) with high accuracy.

The distribution of the amplitudes U_ij is illustrated in Figure 5. The shapes of the light curves of RRLs, represented by the amplitudes of these first few PCs, potentially contain information on the physical properties of the variables themselves, when the pulsation periods of the stars are taken into account. The first amplitude, U₁ (from now on, we omit the j index for simplicity), can be viewed as an analog to the total amplitude in amplitude-period (also called Bailey) diagrams. As described before, additional PCs modify the light-curve shape. Keeping this in mind, it is immediately obvious that, at a given period, the amplitudes U₁ and U₂ separate two sequences, which we can associate with the Oosterhoff type I and II groups of variables (Oo; Oosterhoff 1939; see Catelan & Smith 2015, for a recent review and references). As these groups at least partly correlate with metallicity, we are going to explore the possibility of estimating the metallicity of individual RRLs when their K_S-band light curves are described by the PC amplitudes U_i, in Section 4.

Zoom In Zoom Out Reset image size

The light-curve shapes and amplitudes of pulsating variables vary among bandpasses, due to the difference in the relative contribution of the change of radius and photospheric temperature during the pulsation cycle to the emitted flux at different wavelengths (Catelan & Smith 2015, and references therein). Comparing the K_S and J-band light curves in Figure 1 clearly demonstrates this.

During the course of the VVV survey, each field has been observed only a few times in the J-band. To determine accurate mean J magnitudes for the RRLs, it is necessary to describe the light-curve shapes of the stars in the J-band, i.e., the deviation of the J-band magnitude from its average in each phase of the pulsation cycle. Usually, the difference between the light-curve shapes is ignored, and the K_S-band light-curve shape is used to estimate the average J-band magnitude. However, this method, depending on the light-curve phases of the observations, introduces additional scatter in the derived magnitudes.

As the PC amplitudes U_i provide a concise description of the K_S-band light-curve shapes, we can evaluate whether the J-band light-curve shapes can be approximated with their use. Besides these amplitudes, we consider the period as a possible additional parameter, to assess its effect on variables that otherwise possess the same light-curve shapes in the K_S band.

The ω Centauri photometry of Navarrete et al. (2015) contains only 42 epochs in the J-band; therefore, depending on their periods, some variables have gaps in their J-band light curves at critical phases. A few other stars possess photometric anomalies in their light curves, preventing us from obtaining a good light-curve fit for them in this band. We decided to omit these stars, marked in Table 1, and continue with the analysis of the remaining 87 RRLs. These J light curves are first sampled in the same 100 phase points, as has been done for the K_S-band light curves in Section 2.3, where phase 0.0 is the phase of the K_S-band light-curve minimum for each of the variables.⁹ These curves are aligned to have a mean of 0, but in contrast to the PCA analysis of the K_S-band light curves, they are not normalized.

The training set of J-band light curves is shown in Figure 6. We are searching for the ideal set of variables to describe these light-curve shapes, e.g., the deviations from the mean in each phase. We can consider this as a separate problem at each light-curve phase for which we are looking for a linear model description, as a combination of yet to be determined parameters as the basis. The connection between the J magnitudes in different light-curve phases is that these linear models should depend on the same parameters. Our candidate parameters are the period P, the four PC amplitudes of the K_S-band light curves U₁, U₂, U₃, and U₄, as well as their polynomial combinations up to the third order (P², PU₁, PU₂, PU₃, PU₄, ${U}_{1}^{2}$, etc.).

Zoom In Zoom Out Reset image size

We performed an exhaustive search for the best parameter combination, where we considered various permutations containing up to six candidate parameters. As our sample size is fairly small, we have chosen not to utilize a hold-out set, nor N-Fold Cross Validation for the validation of our results, as the exact (random) choice of the hold-out set or the N-Folds could lead to heavily biased results (for example, if all the short period variables are selected to be in one fold). Therefore, for each candidate parameter set, we have opted to perform Leave-One-Out Cross Validation (LOOCV) the following way:

• 1.  
  we evaluate 87 separate sub-cases, where we omit the light curve of each of the 87 variables;
• 2.  
  in each case, we optimize a linear solution using OLS for the J-band magnitudes, using the candidate parameter set as basis;
• 3.  
  with these solutions, we approximate the J-band light curve of the omitted variable, and calculate the residuals; and
• 4.  
  for each candidate parameters set, we assess the total root mean square error (RMSE) as a performance metric.

We have found that out of all the possible parameter combinations, by only considering the PC amplitudes of the K_S-band light curves U₁, U₂, U₃, and U₄, we reach a total RMSE of our linear prediction of just 0.012 mag. Although this value can still be lowered by including combinations with the period, such as PU₁, we noticed that this mainly decreases the errors for variables coming from ω Centauri, where multiple stars have very similar periods and K_S and J light-curve shapes, while the errors for other RRLs, especially the short period ones, increase drastically. Figure 7 illustrates the residuals both as a function of the pulsation phase (top panel) and for individual stars from the LOOCV (bottom panel).

Zoom In Zoom Out Reset image size

Based on our analysis, we conclude that the J-band light-curve shapes of RRLs can be estimated using their K_S-band PCA amplitudes with high precision. Figure 8 illustrates this process by comparing the K_S-band PCs derived in Section 2.3 (left-hand panel) with the coefficients found by OLS for the prediction of the J-band light-curve shape (right-hand panel), when all of the 87 RRLs with good J light curves are considered.

Zoom In Zoom Out Reset image size

We employed our light-curve fitting method described in Section 3 to determine the K_S-band light-curve parameters of RRab variables listed in the OGLE-IV bulge RR Lyrae catalog (Soszyński et al. 2014). Toward this end, we made use of the data processed by the VISTA Data Flow System (VDFS; Emerson et al. 2004), provided by the CASU. The aperture photometry of detector frame stacks (pawprints) was extracted using the Starlink Tables Infrastructure Library (STIL; Taylor 2006) at the coordinates of the OGLE RRLs. This process has resulted in near-IR light curves for 24217 RRLs.

The VDFS provides photometry in circular apertures of different radii, but due to the source crowding in the Galactic bulge fields, the adopted radii have to be chosen carefully. Generally, smaller apertures provide better photometry for dimmer objects, but in relatively uncrowded cases, this does not always hold true. Furthermore, variables found in overlapping regions of tiles and pawprints can have different offsets between them, but generally the best aperture also minimizes this offset. Therefore, for each star, the optimal aperture has to be chosen on a case-by-case basis.

We utilized our fitting procedure described in Section 3 on the five smallest apertures provided by the VDFS for each OGLE RRL. It has been found that the final value of the Huber cost function divided by the number of data points after the 3.5σ clip is a good indicator of the quality of the photometry obtained with a given aperture with respect to the others. Hence, in this section we are going to utilize the PC amplitudes determined for each star using the aperture for which this value is the smallest.

Due to the large number of RRLs with both OGLE-IV and VVV photometry, we decided to search for correlations between the abundances determined on the former, and the light-curve shapes of the latter. As the main bulge population of RRLs have spectroscopic iron abundances of [Fe/H] ∼ −1 (Walker & Terndrup 1991), this component is expected to be dominant in the calculated photometric metallicities as well.

The OGLE-IV (or OGLE-III, if the former was not available) I-band light curves of RRLs in the common OGLE/VVV sample were fit with a sixth-order Fourier series to calculate the Fourier parameters necessary for the metal abundance formulae of Smolec (2005).

As our goal is to provide a method for estimating metallicity based on the near-IR light-curve parameters, we eliminate all variables where either estimate is uncertain.

By far, the most common reason for the elimination of variables was the presence of the Blazhko effect, which leads to distorted light-curve shapes in RRLs, which are never equivalent to those of non-modulated RRLs (Jurcsik et al. 2002). Prudil & Skarka (2017) studied the Blazhko effect on a subsample of bulge RRLs in the OGLE sample, and found an incidence rate of 40%. We have inspected each I-band light curve and its fit visually to reveal the presence of the Blazhko effect. In dubious cases, we have inspected the Discrete Fourier Transform of the residual light curve for signs of the characteristic side peaks of modulation in the Fourier domain. We opted to eliminate all variables from our selection where inspection of the light curve gave hints of the Blazhko effect, which accounted for approximately the same fraction of stars as found by Prudil & Skarka (2017). Furthermore, this inspection revealed some dubious RRL light curves, as well as many light curves where the Fourier parameters cannot be determined with high accuracy (faint variables, too few light-curve points, gaps in the folded light curve, etc.), which were also removed in order to preserve the purity of the sample.

Besides requiring good-quality light curves in the I-band, the K_S-band data were also subjected to the following quality cuts: we discarded all variables where either the first PC amplitude, U₁, or the PC amplitude ratios ${U}_{2}/{U}_{1}$, ${U}_{3}/{U}_{1}$, or ${U}_{4}/{U}_{1}$ exceeded the limit described in Section 3, as these indicate problematic photometry in the VVV data, such as severe blending.

We utilized Equations (2) and (3) of Smolec (2005) to calculate the metal abundances of the 6215 remaining variables. The top panel of Figure 10 reveals a striking systematic on top of the overall linear trend in the photometric metallicity calculated with Equation (2) of Smolec (2005) as a function of the pulsation period. In the main locus of stars, denoting the Oosterhoff I population of the bulge RRLs, longer-period RRLs seem to have systematically higher metallicities than their shorter-period counterparts. Moreover, a similar trend can be seen on the bottom panel: the lower-amplitude stars have systematically higher metallicities.

Zoom In Zoom Out Reset image size

In globular clusters, for both Oosterhoff groups of fundamental-mode variables, the amplitude decreases with increasing periods. It is hard to conjure up a scenario where systematically higher-metallicity RRLs would only populate the lower-amplitude, longer-period part of the main locus on this diagram, while the lower-metallicity variables occupy the higher-amplitude, shorter-period part. Therefore, we conclude that the metal abundance formula described by Equation (2) of Smolec (2005) suffers from a systematic bias as a function of amplitude/period. This finding is confirmed by a demonstrably monometallic data set: on both panels of Figure 10, empty and filled circles illustrate the I-band photometric metallicities of the Oosterhoff I and II variables in the globular cluster M3, which is well known to be monometallic (Kraft et al. 1992; Cohen & Meléndez 2005), highlighting the systematic offset between these two groups of objects as well.

The second formula given by Smolec (2005) in the form of Equation (3) gives much more consistent abundance estimates as a function of the amplitude. However, to compare the behavior of the calculated abundances of the two Oosterhoff groups of variables, we have to separate them. This has been done using the period–amplitude diagram, utilizing a cut that is dependent on the calculated iron abundance, as illustrated by Figure 11. Comparing the top and middle panels of Figure 12 indicates that there is still a slight trend for the Oosterhoff I variables as a function of the amplitude, as well as an offset with respect to the Oo II stars, which have a calculated average [Fe/H] of −1.05 dex across the whole amplitude range. We correct for this offset using the ridge of Oo I RRLs, derived from a KDE of amplitude bins of the estimated abundances of Oo I stars. The resulting rectified Oo I metallicity estimates are consistent with the Oo II variables, as well as across the whole range of RRL amplitudes (bottom panel of Figure 12).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We do note that not only the iron abundance formulae of Smolec (2005) suffer from obvious biases when applied on a population of monometallic RRLs, nor is the problem limited to the I-band. As a stopgap measure (until new, carefully calibrated photometric abundance formulae become available) we make available a routine on GitHub¹¹ that implements the procedure outlined here to rectify the [Fe/H] values of Oo I RRLs calculated from Equation (3) of Smolec (2005), given the period, as well as the I-band Fourier parameters A₁, A₂ and ${\phi }_{31}$.

To assess whether it is possible to determine the iron abundance of RRLs from the K_S band light-curve shapes, we made use of the K_S-band light-curve parameters of OGLE RRLs in the VVV fields determined in Section 4.1, combined with the photometric metallicities determined in Section 4.2. As the RRL sample on which Smolec (2005) calibrated the relation used to determine the initial [Fe/H] of the RRLs in question only extends down to [Fe/H] = −1.7, using even the rectified abundance estimates below this threshold is uncertain. Nevertheless, because the typical abundance of Oo II globular clusters bearing RRLs is about [Fe/H] ∼ −2.2 dex, we are only discarding the RRLs below this limit. This final cut results in a final sample of 6193 RRLs.

Figure 13 shows the optical photometric metallicity as a function of the period, and the first two PC amplitudes. Some general trends are apparent; for example, as expected from the form of Equation (3) of Smolec (2005), longer-period RRLs have higher iron abundances, but even at the same period and U₁, stars of different U₂ have systematically different iron abundances.

Zoom In Zoom Out Reset image size

As additional features, we also consider Fourier parameters for this regression problem. During the light-curve fitting, the individual PCs are represented as Fourier series (Equation (3)), and the total light curve is given as the sum of these series multiplied by the PC amplitudes (Equation (4)). Therefore, traditional Fourier parameters, such as the Fourier amplitudes and epoch-independent phase differences (Simon & Lee 1981), can be calculated straightforwardly. Of these, we have decided to use the amplitude of the Fourier harmonics A₁, A₂, and A₃, as well as the epoch-independent phase differences ϕ₂₁ and ϕ₃₁. Together with the pulsation period, the total number of independent variables is 10.

We explored the regression routines implemented in scikit-learn (Pedregosa et al. 2011) with the aim of finding a relation to determine the iron abundance from some combination of the considered parameters. We found that regardless of the parameters, the underlying relation is nonlinear and multimodal due to the systematic offset between the Oo I and II variables. Additionally, the heavy excess of bulge RRLs with metallicities around [Fe/H] = −1 dex biases any kind of regression without resampling or weighting the input data.

The best results were achieved using the MLPRegressor routine, implementing a multilayer perceptron regressor, a type of neural network. Artificial neural networks are inspired by and modeled after biological neural networks, and have been used in many branches of science as well as in commercial applications for their ability to learn highly nonlinear relations inherent in many types of data (see Haykin 2011 for a general description of neural networks).

Neural networks have many different parameters that have to be decided (hyperparameters) before the training of the network. As these can vastly influence their performance, we decided to train the MLPRegressor with a grid of different hyperparameter values, and use cross validation to determine the best combination of hyperparameters. The grid of considered hyperparameter values is documented in Table 2. This grid only contains possible neural network architectures up to two hidden layers of neurons; although adding more hidden layers increases the flexibility of the neural network, doing so did not significantly improve the performance in our tests. Besides the parameters of the network, we considered different combinations of the 10 independent variables as an additional hyperparameter.

Due to the heavily biased nature of the bulge photometric metallicity distribution (most stars have abundances [Fe/H] ∼ −1), we implemented a special sampling method. We grouped observations in ten 0.2 dex wide bins (Table 3). Because there are only a few variables with [Fe/H] above 0 and below −2 dex, these were merged with neighboring bins. In each bin, 10 stars were selected randomly to be part of the validation sample. Of the remaining stars, 80 variables were selected randomly with replacement to be part of the training sample, resulting in a training set of 800 data points. As there are less than 90 stars in the most metal-poor and metal-rich bins, some variables at the extremes of the [Fe/H] distributions are selected multiple times due to the selection with replacement. The variables that were not selected for training were added to the cross-validation sample. In case of an uneven distribution of stars between the two extremes in a metallicity bin, the bin was split into sub-bins to guarantee a more even sampling as a function of metallicity.

This separation into training and cross-validation samples guarantees that the training set is large enough for the training of neural networks; the training set is balanced on the total range of possible abundances; by repeating this separation, variables in bins with few stars get selected for the cross-validation sample at least a few times.

This separation was repeated 40 times, resulting in 40 training and cross-validation samples. The neural networks were trained on all 40 training samples with all combinations of the possible hyperparameters detailed in Table 2. In all cases, the predictions of the trained neural networks were calculated for the corresponding cross-validation samples. Finally, for each star and for each hyperparameter combination, the predictions were averaged over the 40 repeats. Every star appeared at least six times in the cross-validation samples.

Due to the imbalanced data set, neural networks with the highest R² score (also called the coefficient of determination) calculated on the complete sample of averaged predictions is heavily biased toward solutions predicting [Fe/H] ∼ −1, irrespective of the values of the independent parameters. Therefore, we have instead selected the best hyperparameter combination by using the sum of the R² scores calculated separately for each of the abundance bins for the averaged cross-validation predictions. The cross-validation result with the highest score is illustrated by the top panel of Figure 14, with features P, A₁, A₂, A₃, ϕ₂₁, and ϕ₃₁; hyperparameters α = 10⁰; two hidden layers with 100 and 20 neurons; activation function relu; and solver lbfgs.

Zoom In Zoom Out Reset image size

After determining the optimal hyperparameters, 100 new training samples were drawn by randomly selecting 80 stars with replacement from each metallicity bin, but without withholding any stars for cross validation. The MLPRegressor was trained on each of these samples, and the predicted metallicity is the average of these 100 trained regressors. The code estimating the iron abundance using these neural networks from the K_S-band light-curve parameters is available on GitHub.¹²

To provide an unbiased evaluation of the performance of a regression model, it is important to test the performance on a test data set not used during the fit for either training or validation. In this case, the model was used to predict the iron abundance from the shape parameters determined directly from the PCA analysis for the variables in Table 1. Then, these predictions can be compared to independent determinations of their abundances from the literature.

The cross-validation results in the top panel of Figure 14 indicate that the method described in Section 4.3 gives reasonable photometric abundance estimates with a scatter of ∼0.22 dex in the −1.7 < [Fe/H] < 0.0 abundance range. The estimated iron abundance is biased for stars below −1.7 toward higher values, probably caused by the combination of overestimation of the I-band photometric abundances, as well as contamination of the parameter space of low-abundance stars with those of higher estimated abundances, due to the higher relative photometric errors in low-amplitude, long-period stars (see Figure 13).

The middle panel of Figure 14 compares the photometric abundances of 17 stars from our Table 1 to the abundances listed in Table 1 of Jurcsik & Kovács (1996). Additionally, the photometric abundance of NSV 660 is compared with the [Fe/H] = −1.31 value determined from the SDSS spectra of the SEGUE Stellar Parameter Pipeline (Lee et al. 2011; Szabó et al. 2014). The two values for this sample of stars generally agree. We note, however, that there is a hint that the abundances of high-metallicity RRLs are slightly underestimated. Furthermore, the [Fe/H] of the most metal-poor star is overestimated, in concordance with the cross-validation results. The abundances of two stars, RR Cet and RR Leo, are overestimated by 0.5 and 0.6 dex, respectively. As there is a single source of photometry for both stars, we believe that the original K_S-band photometry of both stars is affected by systematic trends on their rising branches, causing the estimated light-curve shapes to be deformed, resulting in the overestimation of their abundances.

The bottom panel of Figure 14 compares the K_S-band photometric abundances of the RRLs in ω Cen with their spectroscopic measurements by Sollima et al. (2006). The photometric abundances for five variables from the globular cluster M3 are also compared to the average cluster abundance, [Fe/H] = −1.5 dex, given by Carretta et al. (2009). We do note that so far all abundances were on the metallicity scale defined by Jurcsik (1995). To compare the K_S-band photometric abundance estimates of the ω Cen and M3 variables in a consistent way to their spectroscopic determinations, we convert them to the metallicity scale established by Carretta et al. (2009). The conversion between the two scales can be determined by comparing the common clusters with measured abundances in Table 1 of Jurcsik (1995, J95) and Table A1 of Carretta et al. (2009, C09) (neglecting the cluster NGC 5927, where the difference between the two sources is almost 0.8 dex):

$$\begin{eqnarray}&&{[\mathrm{Fe}/{\rm{H}}]}_{{\rm{C}}09}=1.044{[\mathrm{Fe}/{\rm{H}}]}_{{\rm{J}}95}-0.037,\end{eqnarray} \tag{ 6 }$$

with a residual scatter of 0.1 dex. Furthermore, the iron abundances of Sollima et al. (2006) are also converted to the Carretta et al. (2009) scale by shifting them by −0.02 dex, following Braga et al. (2016).

The photometric abundances of ω Cen stars reproduce the bimodal distribution of iron abundances (Sollima et al. 2006). When taking into account the uncertainties of individual measurements made by Sollima et al. (2006), there is no offset for the metal-rich group at −1.2 dex. For the metal-poor group of variables, there is a hint that our method systematically overestimates abundances by about 0.1 dex, similarly to the results on the top and middle panels of the Figure 14. As for the M3 RRLs, their average K_S-band photometric abundance estimate is [Fe/H] = −1.33 dex, only 0.17 dex higher than the cluster metallicity given by Carretta et al. (2009).

During the calibration of the abundance prediction, we discarded Blazhko variables, as their modulated light curves lead to systematic biases in their calculated photometric abundances (see, e.g., Figure 4 of Jurcsik & Kovács 1996). Nevertheless, a significant fraction of the total population of RRLs suffers from this light-curve modulation. The nature of the Blazhko effect in the near-IR has not been established until very recently, due to the lack of well-populated light curves covering sufficiently long time spans. In Jurcsik et al. (2018), we show that the modulation is indeed present in the K_S-band light curves, but with diminished amplitudes. We have no abundance estimates for the 22 Blazhko RRLs studied in Jurcsik et al. (2018), but we surmise that the majority of them must come from the bulge population of RRLs. We calculated the abundances based on their individual K_S-band mean light curves, for which we get a mean of [Fe/H] = −1.11 dex with a standard deviation of 0.28 dex. These values are in agreement with those of the original bulge sample, where the mean and standard deviation are [Fe/H] = −1.06 dex and 0.26 dex, respectively. Therefore, we surmise that the iron abundances of RRLs should not show systematic biases when estimated with the parameters of the complete, long-term near-IR light curves, even when the Blazhko effect is present.

In summary, the K_S-band light-curve parameters allow the estimation of the iron abundance of RRLs to an accuracy of 0.20–0.25 dex, with slight hints of systematic trends (at the <0.05 dex level) in the range of −1.7 < [Fe/H] < 0 (underestimation for high abundances, overestimation for low ones), but abundances lower than this range are systematically overestimated. To establish better relations, one will need high-quality K_S-band observations of a large sample of RRLs in globular clusters covering a wide metallicity range, and/or K_S-band observations of the hundreds of bright RRLs in the Solar neighborhood, and/or spectroscopic iron abundances for the OGLE field RRLs, with emphasis on the high and low-abundance extremes.