© ESO 2019

1. Introduction

RR Lyrae (RRLs), are very accurate distance indicators and solid tracers of old (age >10 Gyr) stellar populations. The near-infrared (NIR) period-luminosity (PL) relations of RRLs will be the first calibrator of the extragalactic distance scale based on population II stars (Beaton et al. 2016), which will provide an independent estimate of H₀. The NIR bands, when compared with the optical bands, present several advantages. These become even more compelling for variable stars such as RRLs. It is therefore mandatory to fully exploit the advantages brought up by the NIR bands.

i) The NIR bands are less prone to uncertainties in reddening corrections and are less affected by the occurrence of differential reddening. Indeed, the K band is one order of magnitude less affected than the visual band. For this reason, the highly reddened regions of the Galactic centre and of the inner bulge can only be investigated effectively in NIR bands. At these Galactic latitudes, the absorption in the K band becomes of the order of 2.5–3.0 mag (Gonzalez et al. 2012), meaning ∼25–30 mag in the V band. This is well beyond the capabilities of current and near-future optical observing facilities.

ii) The luminosity variation in the optical bands is dominated by variations in effective temperature, while in the NIR bands it is dominated by variations in stellar radius (Madore et al. 2013). This means that the NIR light curves are minimally affected by non-linear phenomena like shock formation and propagation, which cause the appearance of either bumps and/or dips along the light curves. Moreover, the luminosity amplitudes steadily decrease when moving from the optical to the NIR bands and approach an almost constant value for wavelengths equal to or longer than 2.2 μm (Madore et al. 2013). Indeed, the ratio in luminosity amplitudes A(K_s)/A(V) and A[3.6]/A(V) attain values ranging from 0.22 to 0.41 (RRc and RRab, respectively, Braga et al. 2018) and from 0.18 to 0.22 (RRc and RRab, respectively, Neeley et al. 2015).

iii) The typical sawtooth shape of the light curves of RRab in the optical is less sharp in the NIR, where light curves become more symmetrical. This means that, even with a modest number of phase points (eight-twelve, depending on the photometric cadence) the light curve can be well fitted.

The NIR bands, together with these “intrinsic advantages” also bring up several “extrinsic advantages” concerning the RRL distance scale.

i) Solid theoretical (Bono et al. 2001; Marconi et al. 2015) and empirical (Longmore et al. 1986; Bono et al. 2003) evidence indicates that RRLs obey well-defined period-luminosity-metallicity (PLZ) relations in the NIR bands. The slope of the relation becomes steeper and its standard deviation decreases when moving towards longer wavelengths. The RRLs in the optical bands also obey mean magnitude-metallicity (MZ, Sandage 1981a,b) relations, but these are affected by non-linearity and evolutionary effects, and are less precise than the PLZ relations in the NIR bands (Caputo et al. 2000).

ii) In the case that both optical and NIR bands are available, one can adopt the newly developed algorithm REDIME (Bono et al. 2019). REDIME is capable of providing homogeneous and simultaneous estimates of metal content, distance and reddening.

There are, however, some disadvantages to the use of the NIR bands.

i) The identification and the characterization of RRLs is more difficult in the NIR bands, due to the decrease in luminosity amplitude and the less characteristic shape of the light curve when moving from shorter to longer period variables.

ii) Accurate and deep NIR photometry is never trivial, in particular in crowded stellar fields. NIR observations are quite demanding of telescope time, since specific observing strategies must be devised to properly subtract the sky background. This means that NIR observations typically have shallower limiting magnitudes and longer observing runs when compared with optical bands. A practical example of this disadvantage is the comparison between the OGLE and VVV surveys in the Bulge. While the first covers a larger sky area and provides time series with thousands of phase points, the second achieved ∼100 phase points per time series in a smaller area, although being much more capable of piercing the dust in the Galactic plane. A very interesting and promising approach to overcome several of the limitations affecting the NIR bands is to use observing facilities that are assisted by an adaptive optics system. However, these complex detectors have a quite limited field of view, typically of the order of one arcminute or even smaller. This means that they can hardly be adopted for a photometric survey and/or for an efficient detection and characterization of variable stars.

iii) To overcome possible non-linear effects in cameras and/or the saturation of bright stars, and to improve sky subtraction, the NIR images are collected as series of short-exposure images, arranged in specific dithering patterns. Several approaches have been suggested in the literature to perform PSF photometry of NIR images, all of which present both pros and cons.

These limitations of NIR photometry are the main reasons for the development of NIR light-curve templates, both now and in the past. The first concepts of template light curves were provided by Freedman (1988), based on the BVI light curves of classical Cepheids (CCs) in the Local Group galaxy IC1613. However, the first NIR light-curve templates for RRLs were provided more than twenty years ago, in a seminal investigation by Jones et al. (1996, henceforth, J96). Of most importance in this context is that once the period of an RRL is known, preferentially from optical data, together with its optical amplitude and its epoch of maximum light (anchor point), a template provides the opportunity to estimate its mean magnitude on the basis of a single NIR measurement. However, the J96 templates were provided only for the K-band. Furthermore, owing to the limited number of NIR measurements available at that time it was only based on 17 RRab and 4 RRc. J96 divided RRab variables into four subgroups and kept the RRc variables within a single group, therefore obtaining four and one light-curve templates, respectively. However, the bins in luminosity amplitude adopted to split the fundamental pulsators into different sub-groups did not overlap one another (see Fig. 1). It is also worth mentioning that, in the Bailey diagram (luminosity amplitude versus period), the trend of both RRc and RRab luminosity amplitudes is not linear over their typical period range (Cacciari et al. 2005; Kunder et al. 2013). This means that two variables that have the same amplitude might have different periods. Therefore, the use of the luminosity amplitude to discriminate RRLs with different light curve shapes might also have the disadvantage of being affected by degeneracy.

Fig. 1. Bailey diagram, V-band amplitude versus period, for ω Cen RRLs. Blue circles mark RRc variables, while red squares mark RRab variables. The ranges in period for the RRab light-curve templates are indicated by vertical black dashed lines. The blue striped areas show the ranges in amplitude adopted for the light-curve templates by J96. Note that they provided thresholds in the B band, but here they have been rescaled by 1.25, i.e., the typical amplitude ratio (A(B)/A(V)) for RRab variables (Braga et al. 2016).

To overcome some of these intrinsic limitations of the J96 NIR light-curve templates, new approaches have been proposed recently in the literature. It has been suggested by Freedman & Madore (2010) that accurate optical bands for CCs can be transformed into the NIR bands using only a few measurements. The same approach was also applied to RRLs by Beaton et al. (2016), the experiment was limited to a single RRL as a preliminary result of their ongoing investigation based on HST data. The key advantage of this method is that it requires no knowledge of the epoch of maximum light to phase the NIR measurements. More recently, Hajdu et al. (2018) suggested an interesting new method to use a well-sampled K_s-band light curve to estimate the J- and the H-band mean magnitude of an RRL from single-epoch measurements. They used data from the VISTA Variables in the Vía Láctea (VVV) survey and decomposed the K_s-band light curves of 101 RRab variables into orthogonal principal components. Their method also provides estimates of photometric metallicities.

Light-curve templates of RRLs have also been developed in the visual bands. Layden (1998) obtained six V-band light-curve templates, but they were limited to RRab variables. The adopted sample of 103 field RRLs was divided according to the shape of the light curve (Bailey types a and b, plus the phase range of the rising branch). They were used to simultaneously estimate mean magnitude and luminosity amplitude. More recently, optical (ugriz) light-curve templates of RRLs were derived by Sesar et al. (2010) from SDSS photometry of 379 RRab and 104 RRc. They provided 22 RRab templates and two RRc templates for the five ugriz SDSS bands. They found evidence that the shape of the RRab light curves steadily changes when moving from the blue to the red edge of the instability strip, while RRc light curves are dichotomous. They claim that this evidence might suggest the possible occurrence of second-overtone RRLs. However, theoretical models and spectroscopic measurements indicate that shorter period RRc variables are, on average, more metal-rich than the bulk of field RRc variables (Bono et al. 1997; Sneden et al. 2017), providing an alternative explanation to the hypothesis of second overtone RRLs. It is worth mentioning in passing that the light-curve templates by Sesar et al. (2010) were mainly developed for RRL identification – especially within the upcoming LSST survey – rather than to determine their mean magnitudes.

The main aim of this investigation is to provide new NIR light-curve templates for RRLs based on a detailed optical and NIR data set that our group collected for RRLs in the Galactic globular clusters (GGC) ω Cen and M 4, supplemented by literature photoelectric photometry of Milky Way RRLs.

The structure of the paper is as follows. In Sect. 2, we describe the optical and the NIR photometric data sets adopted for the current analysis. In Sect. 3 we deal with the NIR light-curve templates and, in particular, with the criteria adopted to select the period bins and the normalization of the light curves. The analytical form of the light-curve templates are discussed in Sect. 3 together with a detailed discussion of the adopted anchor point to phase NIR measurements. Section 5 is dedicated to the validation of the templates. The validation is based on ω Cen data and OGLE+VVV (Udalski et al. 1992; Minniti et al. 2010) data and it was performed for one- and three-phase point light curves. In Sect. 6 we apply the new NIR templates to the J and K_s light curves of RRLs in the extragalactic GC Reticulum and provide a new true distance modulus determination. We summarize our results in Sect. 7 and briefly outline future developments of the current project.

2. Optical and near-infrared data sets

2.1. Collection and homogenisation of the data sets

We use proprietary – still unpublished until this work – optical and NIR PSF-reduced photometry of RRLs in M 4 (Stetson et al. 2014) and in ω Cen (Braga et al. 2016, 2018). Optical data are in the (Landolt 1983, 1992) system, and NIR data are in the 2MASS photometric system (Skrutskie et al. 2006). We note that the NIR data were binned by epoch, therefore each phase point is actually an average of three to five phase points belonging to the same dithering sequence. The binning process is described in detail in Braga et al. (2018). More insights on the data (telescopes, cameras, and reduction) can be found in Stetson et al. (2014), Braga et al. (2016), and Braga et al. (2018).

These data were supplemented with i) relatively old optical and NIR photoelectric photometry of 26 Milky Way (MW) field RRLs, mostly collected to perform Baade-Wesselink (BW) analysis (Carney & Latham 1984; Cacciari et al. 1987; Jones et al. 1987, 1988a,b, 1992; Barnes et al. 1988; Liu & Janes 1989; Skillen et al. 1989; Clementini et al. 1990; Fernley et al. 1990; Barnes et al. 1992; Cacciari et al. 1992; Skillen et al. 1993a,b), the BW sample, and ii) optical data from long-term photometric surveys (ASAS: Pojmanski 1997; NSVS: Woźniak et al. 2004). We note that the photoelectric data were not available in machine-readable format, therefore we have digitized the tables available in the original papers. Moreover, to deal with such inhomogeneous data sets, we transformed all the photoelectric NIR data to the 2MASS system. We have used the transformations by Carpenter (2001) to convert the magnitudes from the CIT system (Jones et al. 1987, 1988a,b, 1992; Liu & Janes 1989; Barnes et al. 1992), UKIRT system (Skillen et al. 1989), SAAO system (Fernley et al. 1990; Skillen et al. 1993a), and ESO system (Cacciari et al. 1992) to the 2MASS system. We note that more recent transformations between the SAAO and 2MASS system are available (Koen et al. 2007), but these would require measurements in the H-band which are not available for the Fernley et al. (1990) data. On the other hand, the optical photoelectric data are all in the Johnson system. However, we only use these optical data to derive the epoch of the mean magnitude on the rising branch (t_ris), independently for each variable. Therefore, they were not homogenized with the CCD data in the Landolt system.

The key advantage of ω Cen RRLs is that this stellar system contains almost 200 RRLs and they cover a range in metal content of at least one dex (Rey et al. 2000; Sollima et al. 2006). Moreover, it is the only cluster – with the exception of the peculiar metal-rich clusters NGC 6388 (Pritzl et al. 2002) and NGC 6441 (Pritzl et al. 2001) – hosting a sizable sample of long-period (P >  0.7 days) RRLs.

To make a homogeneous optical and NIR data set available to the entire astronomical community, Tables 1 and 2 give, respectively, the UBVRI and JHK_s light curves of 233 RRLs in ω Cen and M 4; Table 2, also provides JHK_s light curves for 21 RRLs in the BW sample. Table 1 is based on the optical data collected during ∼20-year-long campaigns (Stetson et al. 2014; Braga et al. 2016) and are calibrated to the Landolt (UBV) and Kron-Cousins (RI) photometric system. We note that Table 1 contains also literature data (Sturch 1978; Kaluzny et al. 1997, 2004 plus the CATALINA Drake et al. 2009 and ASAS-SN surveys Shappee et al. 2014; Kochanek et al. 2017) which we used to supplement our photometry (more details in Sect. 3.1 of Braga et al. 2016). Table 2 includes objects for which either we collected NIR time series data during ten-year-long observation campaigns (Stetson et al. 2014; Braga et al. 2018) or NIR photometry was available in the literature. The NIR measurmenets listed in Table 2 are in the 2MASS photometric system. We note that the fraction of objects adopted for the NIR light-curve templates is 57% of the total number of objects listed in Table 2.

2.2. Sample selection

To derive accurate and precise NIR light-curve templates we selected from the initial RRL sample the variables satisfying the following criteria.

• 1)

  At least ten phase points in J, H or K_s.

• 2)

  An accurate estimate of t_ris (see Appendix A for the calculation of t_ris), that is, the epoch to which the template is anchored.

• 3)

  A small dispersion (σ ≲ 0.1) of the phase points along the normalized light curve. To derive the light-curve templates, all the light curves were divided by their amplitude (see Sect. 4). The variables with limited photometric accuracy are more likely to increase the dispersion of the normalized light curve, and in turn of the light-curve template. Our approach was conservative: we only included variables with a “clean” trend in the normalized light-curve fit.

• 4)

  Special care was taken to include variables that trace the shape of the light curve of both RRab and RRc when moving from shorter to longer period RRLs. This means the occurrence of either dips just before the phase of maximum light and/or bumps just before the phase of minimum light.

Once we had applied these selection criteria we were left with a sub-sample of 94 RRab and 51 RRc variables. The excluded variables are marked with an asterisk in Table 2. In the following, the selected objects, belonging to ω Cen, to M 4, or to the BW sample are called “template data sample” (TDS); their pulsation properties are listed in Table 3. The reader interested in a more detailed discussion of the approach adopted to derive periods, mean magnitudes, amplitudes, and their uncertainties is referred to Stetson et al. (2014), Braga et al. (2016) and Braga et al. (2018). The photometric properties of field RRLs were derived using the PLOESS polynomial fit (Braga et al. 2018). We note that the number of variables with accurate light curves in all three filters is limited. More specifically, the light-curve templates rely on a number of variables ranging from 142 for the J band to 101 for the H band and 112 for the K_s band. The difference among the three bands is mainly caused by the paucity of H-band data for field and M 4 RRLs. Moreover, the H- and K_s-band light curves have luminosity amplitudes that are half the J-band amplitudes. This means that the photometric scatter in the normalized light curves appears larger.

3. Near-infrared light-curve templates

3.1. Selection of the period bins

We defined the template bins according to the pulsation period of the variable. The reasons are manifold. i) The period is a solid observable, since it can also be firmly estimated for variables showing multi-periodicity (Blazhko, mixed-mode). The same statement does not apply to the luminosity amplitude adopted by J96. ii) The period range covered by TDS variables (0.28–0.47 days for RRcs and 0.39–0.87 days for RRabs) is much larger than the RRL sample adopted by J96 (0.25–0.34 days for RRc and 0.39–0.66 days for RRabs). iii) The optical luminosity amplitude is not a linear function of the period (Cacciari et al. 2005; Kunder et al. 2013). Data plotted in the Bailey diagram (period versus luminosity amplitudes, Fig. 1) clearly show that RRab and RRc variables with similar amplitudes can have significantly different pulsation periods. iv) The period is tightly correlated with the intrinsic parameters (stellar mass, luminosity, effective temperature) of the variable (Bono & Stellingwerf 1994).

We have checked that, for RRc variables, one template bin is enough because the shape of the light curve in the NIR bands is almost sinusoidal over the whole period range. On the other hand, the RRab variables were divided into three period bins for following reasons. i) To improve the sampling along the light-curve template we required at least ten variables per bin for each band, limiting the number of possible period bins. ii) The RRab variables display in the optical Bailey diagram a parabolic trend (Cacciari et al. 2005) when moving from shorter to longer periods. Data plotted in Fig. 1 show that the maximum is located at approximately 0.55 days. iii) We also decided to cut the sample at 0.7 days, because empirical evidence indicates that a transition – both in Blazhko properties and in the optical-to-NIR amplitude ratios – takes place across this boundary (Prudil & Skarka 2017; Braga et al. 2018).

This means that the RRab variables were split into short (RRab1, P ≤ 0.55 days), medium (RRab2, 0.55 <  P <  0.70 days) and long (RRab3, P ≥ 0.70 days) period bins, while the RRc constitute a single period bin (0.28 <  P <  0.47 days). It is worth mentioning that we could have extended the period range of the RRab3 template up to 0.9 days, by including ωCen-V91 and ωCen-V150. However, both variables have light curves with a significantly different shape when compared to the other RRLs in RRab3 sub-sample. More data is required to establish whether RRLs with periods longer than ∼0.87 days require a separated template bin.

Finally we mention that the number of phase points per template bin is 1226 (J), 698 (H), 959 (K_s) for the RRc template; 931 (J), 478 (H), 1125 (K_s) for the RRab1 template; 1662 (J), 995 (H), 1709 (K_s) for the RRab2 template; 440 (J), 284 (H), 512 (K_s) for the RRab3 template. The current data set is more than six times larger than that adopted by J96, and more than 2.5 times larger when considering only the K_s-band data.

3.2. Normalization of the light curves

The NIR light-curve templates we are developing provide the mean magnitude ⟨X⟩ of an RRL with an accuracy of the order of a few hundredths of a magnitude provided that the following data are available:

• i)

  the epoch (t) and the magnitude (X_t) of a phase point;

• ii)

  the period of the variable (P);

• iii)

  the luminosity amplitude in either the V or B band (A(V) or A(B));

• iv)

  the epoch of the anchor point along the light curve. In this investigation, we adopt the epoch of the mean magnitude on the rising branch (t_ris). Inno et al. (2015) demonstrate that for CCs, t_ris is a more precise anchor point than the more commonly-used epoch of the maximum light (t_max).

We performed a number of simulations using optical and NIR light curves for which both t_ris and t_max were available and we find that the former is better defined when moving from the blue to the red edge of the RRL instability strip. The reasons for which t_ris is better defined than t_max are twofold. Firstly, large-amplitude RRab variables characterized by a “sawtooth” light curve show a cuspy maximum. This means that the phases across maximum light occur during a short time range, so an accurate estimate of the epoch of maximum light requires high time resolution. Secondly, some RRc variables display a well-defined dip just before maximum light (U Com, Bono et al. 2000). To properly identify and separate the two maxima, high time resolution is also required for these short-period variables.

The mean NIR magnitude, ⟨X⟩, of a variable for which the aforementioned parameters are available can be estimated by using the following relation:

$$\begin{array}{c}\hfill ⟨X⟩=X_t-A(X)·T({\varphi }_t)\end{array}$$(1)

where ${\varphi }_t={\displaystyle \frac{t-t_{\mathrm{ris}}}{P}}$ is the difference in phase between the NIR phase point that was observed and the epoch, t_ris, of the anchor point, while A(X) is the luminosity amplitude in the X band. We note that the latter is typically unknown, but it can be estimated from the optical amplitude and empirical NIR-over-optical amplitude ratios (Braga et al. 2018). We note also that the light-curve templates must be normalized.

To generate the normalized light-curve templates, we adopted the magnitudes m_ijk of the TDS variables, marked with an asterisk in Table 2 and listed in Table 3, where i indicates the ith phase point of the empirical light curve, j indicates the band (1 for J, 2 for H, and 3 for K_s), and k indicates the kth RRL in the TDS sample. We have transformed all the empirical m_ijk measurements into normalized magnitudes M_ijk by subtracting from each kth RRL its mean magnitude in the jth band (see Table 3 ⟨m_jk⟩) and by dividing for the jth band amplitude (see Table 3A_jk) according to the following relation:

$$\begin{array}{c}\hfill M_{\mathit{ijk}}={\displaystyle \frac{m_{\mathit{ijk}}-⟨m_{\mathit{jk}}⟩}{A_{\mathit{jk}}}}·\end{array}$$(2)

Figures 2 and 3 show the final normalized light curves as a function of the pulsation phase for the TDS sample.

Fig. 2. Left panels: from left to right the different panels display the Fourier fits of the normalized JHKs light curves. Top to bottom panels: four (RRc, RRab1, RRab2, and RRab3) different ranges in period. Right panels: same as the left, but for the residuals of the normalized light curves with the Fourier fits. The median and standard deviation of the median are labelled in red.

Fig. 3. Left panels: from left to right the different panels display the PEGASUS fits of the normalized JHKs light curves. From top to bottom the panels show the four (RRc, RRab1, RRab2, and RRab3) different ranges in period. Right panels: same as the left, but for the residuals of the normalized light curves with the PEGASUS fits. The median and standard deviation of the median are labelled in red.

4. Analytical fits to the light-curve templates

Once the normalized light curves for the three NIR bands and for the different period bins were derived, we performed an analytical fit of the light-curve templates. We adopted two different fitting functions: Fourier series (Sect. 4.1) and periodic Gaussians (PEGASUS, Sect. 4.2). We chose the quoted functions for several, key reasons. Firstly, they are analytical functions. This means that they can be easily used as fitting functions, on top of the classical use as templates. Secondly, Fourier series are very easily fitted when the number of phase points is similar to those in our database and the coefficients are quite robust. Thirdly, Fourier and PEGASUS series were already successfully used in an analogous paper for NIR light-curve template of CCs (Inno et al. 2015). Finally, PEGASUS series are a sound complement to Fourier series since the former is not affected by the presence of secondary bumps (ripples) affecting the latter, (see Sects. 4.1 and 4.2).

4.1. Fourier fit

We have fitted the normalized light curves with Fourier series of the ith order

$$\begin{array}{c}\hfill F(\varphi )=A_0+{\mathrm{\Sigma }}_i{A}_{i}cos(2\pi i\varphi -{\varphi }_i)\end{array}$$(3)

with i ranging from two to seven. The red lines plotted in the left panels of Fig. 2 show the individual fits for the three different bands and for the four light-curve templates. The coefficients and the standard deviations of the different fits are listed in Table 4.

We note that the agreement between the analytical fits and observations is, within the errors, quite good over the entire pulsation cycle. In particular, the fits properly represent the data across the phases of minimum light in which the variation of the luminosity is more cuspy. Interestingly enough, we found that the residuals between the normalized light curves and the Fourier fits plotted in the right panels of the same figure are vanishing. They are typically smaller than the fourth decimal place. Moreover and even more importantly, the residuals do not show any phase dependence within the standard deviation (dashed red lines) of the analytical fits. In this context it is worth mentioning that the light-curve templates derived by J96 were obtained using second-order Fourier fits for the RRc variables and sixth-order Fourier fits for the RRab variables. We used different orders for almost all the period bins, however, we adopted the sixth order for the fit of the RRab3 K_s-band templates. This template includes roughly the same number of variables as the RRab1 template by J96 (A(B) < 1.0 mag), however the coefficients of the fit are significantly different.

4.2. PEGASUS fit

We also performed an independent fit of the normalized light curves using a series of periodic Gaussians, presented in Inno et al. (2015) with i ranging from two to six.

$$\begin{array}{c}\hfill P(\varphi )=A_0+{\mathrm{\Sigma }}_i{A}_{i}exp(-sin({\displaystyle \frac{\pi (\varphi -{\varphi }_i)}{{\sigma }_i}}{)}^2)·\end{array}$$(4)

Data plotted in the left panels of Fig. 3 show that PEGASUS fits follow the variation of the normalized light curves quite well over the entire pulsation cycle. This applies not only to the RRc, but also to the RRab light-curve templates. The main difference between the fits based either on PEGASUS or on Fourier series is that the former display a smoother variation over the entire pulsation cycle, while the latter show several small bumps/ripples. The NIR light curves with accurate photometry and very well sampled light curves do not display these bumps. This suggests that the bumps/ripples are spurious variations of the order of a few thousandths of a magnitude among the different variables included in a given period bin.

The residuals between the normalized light curves and the PEGASUS fits are plotted in the right panels of the same figure. They are of the order of a few thousandths, that is, slightly larger than the residuals of the Fourier fits. The difference is mainly due to the smoothness of the PEGASUS fits compared with the Fourier fits. Coefficients and the standard deviations of the fits are listed in Table 5.

4.3. Phases of maximum and minimum along the light-curve template

Although there are legitimate reasons supporting the idea that t_ris is easier to derive than the epoch of maximum light, t_max, and it provides a more precise epoch of reference, we are aware that all the recent surveys adopt t_max as the reference epoch for RRLs and other variable stars. For this reason we also provide the phases of both maximum and minimum (ϕ_max and ϕ_min) of the current light-curve templates (see Table 6). These pulsation phases – which can be considered typical – provide the opportunity to use the current templates to estimate the mean magnitude of variables for which only t_max and/or t_min is available in the literature.

5. Validation of the light-curve templates

5.1. Validation based on ω Cen RR Lyrae

To validate the light-curve templates, we need optical and NIR light curves of RRLs from which we can derive accurate estimates of their photometric properties (mean magnitudes, amplitudes and t_ris). However, to perform an independent check we cannot use RRLs in the TDS (Table 3). Therefore we defined a template validation sample (TVS) including four ωCen RRLs: ωCen-V20 (RRc), ωCen-V57 (RRab1), ωCen-V107 (RRab2) and ωCen-V124 (RRab3). The selection of the four TVS RRLs was based on the following criteria: i) – the TVS RRLs have well-sampled J-, H- and K_s-band light curves and cover the four-light-curve templates we develop; ii) – the estimate of epoch of reference (t_ris, t_max) is very accurate.

We estimated the mean NIR (⟨JHK_s⟩_best) magnitudes of the TVS RRLs by fitting the light curves in flux and then converting the mean flux to mean magnitude. To estimate the mean NIR magnitude (⟨JHK_s⟩_templ) with the light-curve template we followed two different paths based either on single phase point (Sect. 5.1.1) or on three independent phase points (Sect. 5.1.2). The key idea is to estimate the accuracy of the light-curve templates from the difference Δ⟨JHK_s⟩ between the measured (⟨JHK_s⟩_best) and the estimated (⟨JHK_s⟩_templ) mean magnitudes. The mean NIR magnitudes will be estimated from the Fourier and PEGASUS fits for both the single-phase point and the triple-phase points method. To discriminate among them we add suffixes to the subscript of the mean magnitudes ⟨JHK_s⟩_{templ[P/F][1/3]}, where [P/F] indicates that we used either the PEGASUS or the Fourier fit, and [1/3] indicates that we used either the single-phase point or the triple-phase point approach. Finally, to provide a more quantitative comparison with the light-curve template available in the literature we also fit the TVS RRLs with the J96 templates.

5.1.1. Light-curve templates applied to a single phase point

We extracted 100 phase points (ϕ_i,JHK_s(i), where i runs from 1 to 100) starting from an evenly-spaced grid of phases ϕ_i = [0.00, 0.01, ... 0.99]. For each ϕ_i, we generated a random magnitude JHK_s(i) = JHK_s(fit(ϕ_i)) + rσ. The two components of this extracted light curve are i)JHK_s(fit(ϕ_i)), which is the value of the fit of the light curve at the phase ϕ_i, and ii)rσ, which simulates random noise: σ is the standard deviation of the phase points around the fit and r is a random number extracted from a normal distribution.

We also derived, by applying the template with Eq. (1), 100 estimates of ⟨JHK_s⟩_templ(i), one for each extracted phase point. Subsequently, we estimated the median and the standard deviation of the median over the 100 ⟨JHK_s⟩_templ(i) extractions. Figures 4, 5 and 6 display the extracted phase points and the fits based on the light-curve templates in the J, H and K_s bands.

Fig. 4. Panels a and b: Black crosses show the randomly extracted Ks-band phase points over the light curve of the RRc variable ωCen-V83. Grey dashed lines display the fit of the light-curve template (Fourier, a; PEGASUS, b) to the extracted phase points. The ID of the RRL is labelled. Panels c and d: Same as panels a and b, but for the RRab variable ωCen-V107. The RRab1 light-curve templates were adopted. Panels e and f: Same as panels a and b, but for the RRab variable ωCen-V125. The RRab2 light-curve templates were adopted. Panels g and h: Same as panels a and b, but for the RRab variable ωCen-V15. The RRab3 light-curve templates were adopted.

Fig. 5. Same as Fig. 4, but for the H band light-curve templates.

Fig. 6. Left (a,d,g,j) and middle (b,e,h,k) panels are the same as Fig. 4, but for the Ks band light-curve templates. The right panels (c,f,i,l) display the fit based on the Ks band the J96 light-curve template.

The estimates of ⟨JHK_s⟩_templ – using the Fourier, PEGASUS and J96 templates – of the TVS RRLs are listed in Table 10. The same Table also gives the difference in magnitude (Δ⟨JHK_s⟩) among the different fits.

It is worth noting (see Table 8, Cols. 2–4) that the mean of the residuals with respect to the measured magnitudes is at most 0.010 mag for all the templates. In all cases, the standard deviations are larger than the residuals, meaning that the latter can be considered null within the dispersion. The largest residuals are found in the H band for the RRab1 template: the mean magnitudes estimated from the templates are ∼0.01 mag fainter than the measured mean magnitude. This happens because the fit of the H-band light curve has minor deviations from the light-curve template, and the extracted single phase points follow these deviations. We note that, in performing this test, we are maximizing the uncertainty, since the error on the individual phase points is estimated as a Gaussian distribution with a σ equal to the standard deviation of the analytical fit. Indeed, we find that when using the individual measurements the residuals are systematically smaller.

The comparison between the new and old K_s-band templates indicates that the former are on average better than the latter. Indeed, the residuals in the longest period bin (RRab3) of the new templates are one order of magnitude smaller than for the J96 template (–0.001 [Fourier]/0.000 [PEGASUS] mag vs. –0.011 mag). We note, however, that the standard deviations are of the same order of magnitude of the difference in offset between our templates and those of J96. Moreover, the standard deviation of the current RRab1 period bin is more than a factor of two smaller than for the J96 template (0.016 [Fourier]/0.016 [PEGASUS] mag vs 0.038 mag). A glance at the data plotted in the right column of Fig. 6, and in particular in the panels d to f, clearly shows the difference.

5.1.2. Light-curve templates applied to three phase points

The application of the NIR light-curve templates to individual NIR measurements does require the knowledge of three parameters: i) the period, ii) the luminosity amplitude, and iii) the epoch of the anchor point (t_ris). The third parameter poses a severe limitation for RRLs because their periods range from a quarter of a day to less than one day. This means that either the pulsation period and the epoch of the anchor point have been estimated with very high accuracy (∼one part per million) or the separation between the time at which the optical and the NIR photometry were collected must be shorter than a few years.

To overcome this limitation we decided to perform a number of tests by assuming that three independent NIR measurements were available. The advantage of this approach is that the light-curve template is used as a fitting function. The approach is quite simple and follows the following steps: i) an estimate of the NIR luminosity amplitude A(NIR) using the optical to NIR amplitude ratios available in the literature; ii) a least-squares fit of a light curve including at least three phase points, minimizing the χ² of two parameters: a shift in phase (Δϕ) and a shift in magnitude (Δmag). The functions to be minimized are:

$$\begin{array}{c}\hfill F(\varphi ;\mathrm{\Delta }\varphi ,\mathrm{\Delta }\mathrm{mag})=\mathrm{\Delta }\mathrm{mag}+A(NIR)·(A_0^F+\\ \hfill {\mathrm{\Sigma }}_i{A}_i^{F}cos(2\pi i\varphi -{\varphi }_i-\mathrm{\Delta }\varphi ))\end{array}$$(5)

and

$$\begin{array}{cc}\hfill P(\varphi ;\mathrm{\Delta }\varphi ,\mathrm{\Delta }\mathrm{mag})=& \mathrm{\Delta }\mathrm{mag}+A(NIR)·(A_0^P+{\mathrm{\Sigma }}_i{A}_i^{P}exp\hfill \\ \hfill & (-sin({\displaystyle \frac{\pi (\varphi -{\varphi }_i-\mathrm{\Delta }\varphi )}{{\sigma }_i^P}}{)}^2)),\hfill \end{array}$$(6)

for the Fourier and PEGASUS templates, respectively. To further investigate the difference between new and old light-curve templates, the same minimization was also performed using the J96 templates:

$$\begin{array}{cc}\hfill J(\varphi ;\mathrm{\Delta }\varphi ,\mathrm{\Delta }\mathrm{mag})=& \mathrm{\Delta }\mathrm{mag}+A(NIR)·(A_0^J+{\mathrm{\Sigma }}_i{A}_i^{J}cos\hfill \\ \hfill & (2\pi i\varphi -{\varphi }_i-\mathrm{\Delta }\varphi )).\hfill \end{array}$$(7)

To validate the templates with this approach, we generated 100 triplets of phase points (ϕ_ij,JHK_s(ij), where i runs from 1 to 100 and j from 1 to 3). The phases are randomly extracted from a uniform distribution between 0 and 1. The extracted magnitudes, JHK_s(ij), were treated following the approach discussed in Sect. 5.1.1.

Once the 100 three-phase point light curves were generated, we performed the fits using Eqs. (5), (6) and (7). The individual K_s-band fits are displayed in Fig. 7. We computed 100 estimates of the mean ⟨JHK_s⟩_templ(i) magnitudes as the integral in flux over the template fits. The final mean magnitude (⟨JHK_s⟩_templ) and its uncertainty were determined as the median and the standard deviation of the median over the 100 random estimates of ⟨JHK_s⟩_templ(i) (see Sect. 5.1.1). Table 8 also shows the magnitude differences Δ⟨JHK_s⟩ between the template estimates of the mean magnitudes ⟨JHK_s⟩_templ(i) and the best estimates of the mean magnitudes based on the fit of the light curve (⟨JHK_s⟩_best).

Fig. 7. Top-left: Fourier, PEGASUS and J96 template fits applied to RRc variable ωCen-V83. Grey crosses show the randomly extracted phase points. The red crosses display the three phase points of a single random extraction. The horizontal red line shows the mean magnitude of the variable based on the flux average fit of the empirical data. The black dashed curve and the horizontal black long-dashed line show the fit with the Fourier template and its mean magnitude. The magenta dashed curve and the horizontal magenta long-dashed line display the same, but for the PEGASUS fit. The blue dotted curve and the horizontal blue dotted line display the same, but for the J96 fit. Top-right: Same as the top-left, but for the RRab variable ωCen-V107 (RRab1). Bottom-left: Same as the top-left, but for the RRab variable ωCen-V125 (RRab2). Bottom-right: Same as the top-left, but for the RRab variable ωCen-V15 (RRab3). We note that, in these panels, the light curves are not phased using tris, but to an arbitrary epoch (HJD = 2 350 000) to underline that the three phase points method is independent of the reference epoch.

Data plotted in Fig. 7 show that the residuals are similar to the fits based on a single phase point. Indeed, the residuals are, within the standard deviations, zero. However, the standard deviations of the template fits based on three phase points are larger than those based on a single phase point. The difference is mainly caused by the fact that the three randomly-selected phase points span, in some of the extractions, a very small range in pulsation phase (Δϕ ≤ 0.05, see Fig. 8). This is also the reason why the residuals are correlated with the difference in phase between the two closest points in phase (Δϕ).

Fig. 8. Left: Same as Fig. 7, but for random extraction on the light curve of the RRab variable ωCen-V125 (RRab2). In this specific case the three randomly extracted phase points are close in phase. The J96 fit does not provide an accurate estimate of the mean magnitude. Right: Same as the left, but for an extraction where the three randomly extracted phase points have a smaller difference in phase. In this case the light-curve templates do not provide an accurate estimate of the Ks band mean magnitude.

The current findings indicate that the light-curve templates used as fitting curves provide accurate mean magnitudes when i) the distance between the phase points is at least 0.1 pulsation cycles. Otherwise, we suggest averaging the two close phase points. ii) the number of available phase points is modest, that is, larger than two, but smaller than a dozen. Classical analytical fits (e.g. Fourier, Spline, PLOESS, or PEGASUS) become more accurate for a larger number of measurements.

5.2. Validation based on OGLE + VVV RR Lyrae

An independent path to validate the current light-curve templates is offered by the two different long-term photometric surveys collecting time-series data in the optical (OGLE, Udalski et al. 1992) and the NIR (VVV, Minniti et al. 2010) of a significant fraction of the Galactic bulge. The photometric catalogues provided by these surveys can be simultaneously used to validate the K_s-band templates. We note that we were not able to validate the J- and H-band templates because the VVV survey only collected K_s-band time series. Indeed, the current NIR photometric surveys either collected time series data only for limited sky regions (2MASS, Quillen et al. 2014) or they are not yet available (VMC, Cioni et al. 2011).

The validation relies on the OGLE-IV catalogue of 38,257 Bulge RRLs Soszyński et al. (2014). Using a searching radius of 2″, we found 2,517 matches in the VVV point source catalogue. We used a very small searching radius because this provides a faster selection of the good matches. Obviously, the completeness is modest, but the validation only requires a few variables per period bin. Among them we selected 22 RRLs and the criteria we adopted for the selection are the following: i) good coverage of the K_s-band light curve, and in fact they all have at least 38 phase points (80% of them have at least 49 phase points); ii) good coverage of both the V- and the I-band light curve to provide accurate estimates of the luminosity amplitudes (A(V), A(I)) and of the epochs of the mean magnitudes on the rising branch (t_ris(V), t_ris(I)). Optical and NIR light curve phasing was performed using the pulsation period provided by OGLE. The distribution of these variables among the different period bins is the following: RRc (six), RRab1 (five), RRab2 (five), and RRab3 (six). These variables were called the “bulge template validation sample” (BTVS) and their pulsation properties are listed in Table 9.

The validation with the BTVS RRLs follows the approach adopted for the ω Cen RRLs (see Sect. 5.1). The key idea is to compare the mean magnitude estimated by using the template (⟨K_s⟩_templ) with the mean magnitude evaluated by using the K_s band measurements (⟨K_s⟩_best). For these objects we compare eight independent estimates of ⟨K_s⟩_templ, because we apply Fourier and PEGASUS fits to the light curve parameters based on the V- and on the I-band data. Moreover, the validation is applied to both single phase points and triple phase points. We add suffixes to the subscript of ⟨K_s⟩_{templ[P/F][V/I][1/3]}. where [P/F] indicates that we used either the PEGASUS or the Fourier fit, [V/I] indicates that we used either the V- or the I-band data, and [1/3] indicates that we used either the single phase point or the triple phase points.

The two methods are identical to those described in Sects. 5.1.1 and 5.1.2. The only difference is that in this case we have more than one RRL per template bin. Therefore, we also estimated the median difference ⟨Δ⟨K_s⟩⟩ for all the RRLs in the period bin. The results are listed in Table 10. Fig. 9 displays the fits to four BTVS RRLs, one for each template bin.

Fig. 9. Panels a–c: black crosses represent the randomly extracted Ks-band phase points over the light curve of the OGLE RRc variable 11254. Grey dashed lines display the template fits to the individual phase points. The Fourier (panel a), PEGASUS (panel b), and J96 (panel c) RRc light-curve templates are also displayed. The ID of the RRL is labelled. Panels d–f: same as panels a–c, but for the OGLE RRab variable 14578. The RRab1 light-curve templates are displayed. Panels g–i: same as panels a–c, but for the OGLE RRab variable 14806. The RRab2 light-curve templates are displayed. Panels j–l: same as panels a–c, but for the OGLE RRab variable 14958. The RRab3 light-curve templates are displayed.

In this context it is worth mentioning that two (OGLE ID: 34618, 11992) out of the 22 BTVS RRLs, both belonging to the RRab2 period bin, are Blazhko RRLs. The amplitude modulation is 0.2 mag in the I band and 0.3 mag in the V band. The Blazhko modulation does not significantly affect the mean magnitude (Δ⟨K_s⟩) for two main reasons. Firstly, the OGLE data are well sampled and we could estimate the average amplitude over the Blazhko cycle. Secondly, the Blazhko variables with extreme amplitude modulation, that is, 0.5 mag in V, and sampled only across the phases of the maximum will be affected by an error of the order of 0.4 mag in V amplitude. The impact of this amplitude uncertainty on the mean magnitude estimated by using the template is minimal, indeed it is of the order of 0.002 mag in the J band and even smaller for the other bands. We note, however, that this limitation becomes severe for the J96 RRab templates, because the different light-curve templates are based on the luminosity amplitude. The use of a wrong template causes a systematic error in the mean magnitude of the order of a few hundredths of a magnitude.

6. Application of the new light-curve templates to Reticulum RRLs

Reticulum is an extragalactic globular cluster associated with the halo of the Large Magellanic Cloud (LMC). It hosts a sizable sample of RRLs (32 in total as by Walker 1992) and it is an interesting workbench, because the J96 light curve templates were adopted by Dall’Ora et al. (2004) to derive the mean K_s-band magnitudes of 30 RRLs that were observed with SOFI at NTT. However, the mean J-band magnitudes were estimated as the mean of the measurements. The number of measurements was limited, typically 46 unbinned phase points, which means on average ten binned phase points (see below). This means that the classical analytical fits (spline, Fourier series) could be applied. Moreover, the J-band light curve templates were not available. For these reasons, the authors focused their cluster distance determinations only on the K_s-band PL relation. The new light-curve templates will be used to provide new J- and K_s-band mean magnitudes, new NIR PL relations and, in turn, new cluster distance determinations.

6.1. Phasing of the data and application of the light-curve templates

We plan to use the photometric data collected by Dall’Ora et al. (2004), but we derive new NIR (JK_s) curves. In particular, we plan to take advantage of the new pulsation periods and epoch of the anchor point recently provided by (Kuehn et al. 2013). Moreover, the SOFI JK_s-band data were binned using the same approach adopted in Braga et al. (2018). The data collected in one dither pattern were binned into a single phase point using a time interval of 108 s. The binned J- and K_s-band light curves have a number of phase points ranging from ten to fourteen. The J- and K_s-band light curves of three variables, V10 (RRc), V19 (RRab1) and V5 (RRab2), are displayed in Fig. 10 together with the template fits (black dashed lines) and the mean magnitude (green solid line).

Fig. 10. Top: J (left panel) and Ks (right panel) band light curve for the RRc variable V10. The red diamonds display the binned phase points. The black dashed curves show the light-curve template applied to the individual binned phase points. The thick green line displays the mean magnitude listed in Table 11. Middle: same as the top panel, but for the RRab variable V19. Bottom: same as top panel, but for the RRab variable V5.

We have folded the light curves with the periods published by Kuehn et al. (2013). However, the decimal places provided in their Tables 1 and 2 are limited and for seven RRLs (V3, V4, V11, V15, V24, V28, and V32, using the new notation introduced by Kuehn et al. 2013), the folded light curves show significant phase drifts. Therefore, for these RRLs we estimated our own periods, based on their V-band light curves (see Table 11).

Subsequently, we estimated t_ris from the V-band light curves provided by Kuehn et al. (2013). We fit the optical light curves using the PLOESS method described in Braga et al. (2018). We find that the difference between our V-band mean magnitudes and those provided by Kuehn et al. (2013) is negligible, with a mean of 0.003 mag, a standard deviation of 0.012 mag and a maximum difference of 0.035 mag. On the basis of the new periods and of the new epochs (t_ris), we folded the NIR light curves.

We note that Reticulum hosts six mixed-mode RRLs (RRd) and we have NIR data for five of them (except V32). We do not provide templates for this type of variable, but since the dominant mode is the first overtone, we decided to apply the RRc light-curve template to these variables.

To apply the template, we need an estimate of the optical amplitudes of the RRLs and of the NIR-to-optical amplitude ratios (Braga et al. 2018) to rescale the template function. We decided to adopt our own V-band amplitudes – estimated from the PLOESS fits derived in Sect. 6.1 – because they differ from those published by Kuehn et al. (2013). The mean difference ΔA(V) = A(V)_our − A(V)_K13 is –0.08 mag, with a standard deviation of 0.07 mag and a maximum difference of –0.32 mag. We obtained smaller luminosity amplitudes because, for Blazhko and RRd variables, we did not fit the brightest or faintest envelopes of the data (Kuehn et al. 2013) since we are interested in the application of the template to determine their NIR mean magnitudes.

We then applied to each phase point of the NIR binned light curve both the PEGASUS and the Fourier light-curve templates. This means that we estimated two mean magnitudes (⟨J⟩_i, ⟨K_s⟩_i) per phase point, where i indicates the ith phase point. Interestingly enough, the Fourier and the PEGASUS templates provide, within the photometric uncertainty of the individual phase points, similar estimates of both ⟨J⟩_i and ⟨K_s⟩_i. The final values of ⟨J⟩ and ⟨K_s⟩ are the medians of all the ⟨J⟩_i and ⟨K_s⟩_i. These values are listed in Table 11, together with their standard deviations.

6.2. New empirical J and K_s PL relations and the distance to Reticulum

We derived the PL relations in the J and K_s band after correcting the NIR mean magnitudes for reddening. Following the same arguments of Muraveva et al. (2018a), we adopted the cluster reddening (E(B − V) = 0.03 ± 0.02 mag) originally derived by Walker (1992). We also adopted R_V = 3.1 and the optical-to-NIR extinction ratios by Cardelli et al. (1989). We note that in the current PL relations the periods of RRc and RRd variables were “fundamentalized”, meaning that we adopted log P_F = log P_FO+0.128 (Kuehn et al. 2013). We obtained the following PL relations, where J₀ and K_s0 indicate the un-reddened magnitudes:

$$\begin{array}{c}\hfill J_0=(17.78\pm 0.05)-(1.58\pm 0.17)·logP\end{array}$$(8)

$$\begin{array}{c}\hfill K_{s0}=(17.29\pm 0.04)-(2.40\pm 0.15)·logP.\end{array}$$(9)

The coefficients of the current empirical PLK_s relation and their standard deviations are, within the errors, very similar to those obtained by Dall’Ora et al. (2004). The standard deviation of the PLJ relation is larger than in the PLK_s relation (0.05 vs. 0.04 mag), as suggested by theoretical predictions (0.06 mag, Marconi et al. 2015). Finally, we have estimated the true distance modulus (μ) of Reticulum using the new NIR mean magnitudes (J,K_s) and the theoretical Global PLZ relations provided by Marconi et al. (2015) and Marconi et al. (in prep.). We have adopted the spectroscopic iron abundance obtained by Suntzeff et al. (1992) from Reticulum red giants, transformed into the Carretta et al. (2009) metallicity scale ([Fe/H] = –1.70). We found μ_J = 18.47 ± 0.10 (rand.) ± 0.03 (syst.) and μ_Ks = 18.49 ± 0.09 ± 0.05 mag, where the first is the standard error of the mean and the second the standard deviation. The latter was computed as the squared sum of the average uncertainty on the mean magnitudes only, since the uncertainty on the extinction and the propagation of the uncertainties in the calibrating PLZ coefficients vanish when square-summed.

The true distance modulus obtained by Dall’Ora et al. (2004) from the same data, but using a different theoretical K_s-band PLZ (Bono et al. 2003) relation, was μ = 18.52 ± 0.05 mag. The three distance determinations agree wich each other quite well, and indeed the difference is within 1σ.

The distance to Reticulum was estimated by Kuehn et al. (2013) using the visual mean magnitude-metallicity relation relation provided by Catelan & Cortés (2008), a cluster metallicity of [Fe/H] = –1.66 (Mackey & Gilmore 2004, in the Zinn & West 1984 scale) and a cluster reddening of E(B − V) = 0.016 mag (Schlegel et al. 1998). They found a true distance modulus of 18.40 ± 0.20 mag. They also adopted the I-band PL relation provided by Catelan et al. (2004), the same cluster reddening and the Cardelli et al. (1989) reddening law and they found a true distance modulus of 18.47 ± 0.06 mag.

The Reticulum true distance modulus was more recently estimated by Muraveva et al. (2018a) using mid-infrared (MIR) mean magnitudes of 24 ([3.6]) and 23 ([4.5]) RRLs, respectively, collected with IRAC at Spitzer. They found true distance moduli of μ = 18.32 ± 0.06 mag ([3.6]) and 18.34 ± 0.08 mag ([4.5]) mag, adopting two empirical zero-points based on Gaia DR1 (Gaia Collaboration 2016a) and Gaia DR2 (Gaia Collaboration 2018; Clementini et al. 2019) trigonometric parallaxes and a reddening of E(B − V) = 0.03 mag (Walker 1992), according to a detailed review of the reddening estimates of Reticulum available in the literature. They also adopted a third independent zero-point based on HST (Benedict et al. 2011) trigonometric parallaxes for five field RRLs, and found that this calibration provides distances that are 0.10 mag larger than those based on the Gaia calibrations. We note that Muraveva et al. (2018a) adopted a different metal content ([Fe/H] = –1.66, Mackey & Gilmore 2004, in the Zinn & West 1984 scale), but the difference in cluster metallicity affects the distance only at the level of 0.01 mag.

The cluster distance found by Muraveva et al. (2018a) is smaller than the geometric distance to the LMC found by Pietrzyński et al. (2013) (μ_LMC = 18.493 ± 0.008 ± 0.047 mag) from late-type eclipsing binaries and by Inno et al. (2016, μ_LMC = 18.48 ± 0.10 mag) from CCs with optical/NIR (VIJHK_s; ∼4000) and MIR (w₁, WISE photometric system; ∼2600)

are young (t <  300 Myr), intermediate-mass stars and mainly trace the disk/bar of the galaxy. On the basis of their relative distances (Inno et al. 2016) found an LMC depth of the order of ∼±0.2 mag. This suggests that the intrinsic spread in distance along the line of sight is roughly the 10% of its distance (∼±5 kpc).

To discuss the position of Reticulum compared with LMC barycentre in greater detail, we provide independent and homogeneous distance moduli based on both optical and MIR measurements available in the literature. This approach is further strengthened by the recent findings by Muraveva et al. (2018b). Based on a large sample of Gaia DR2 trigonometric parallaxes (Arenou et al. 2018), Muraveva et al. suggest that the coefficients of the metallicity term predicted by pulsation models agree quite well with observations. We adopted the MIR mean magnitudes provided by Muraveva et al. (2018a) and the MIR theoretical PLZ relations provided by Marconi et al. (2015); Neeley et al. (2017, see Fig. 11). According to (Muraveva et al. 2018b), we adopted E(B − V) = 0.03 mag and we find the following empirical PL relations

$$\begin{array}{c}\hfill [3.6{]}_0=(17.24\pm 0.06)-(2.12\pm 0.21)·logP\end{array}$$(10)

$$\begin{array}{c}\hfill [4.5{]}_0=(17.13\pm 0.08)-(2.52\pm 0.29)·logP\end{array}$$(11)

and true distance moduli of μ_[3.6] = 18.30 ± 0.06 ± 0.05 and of μ_[4.5] = 18.31 ± 0.08 ± 0.08 mag. The new true distance moduli are in remarkable agreement with distances provided by Muraveva et al. (2018a). The distances based on MIR mean magnitudes are systematically smaller then those based on NIR mean magnitudes, but the difference is of the order of 1σ. To further investigate the possible systematics affecting the current distance determinations we also estimated the true distance modulus from the I-band mean magnitudes provided by Kuehn et al. (2013). We find the following empirical PL relation

$$\begin{array}{c}\hfill I_0=(18.30\pm 0.03)-(1.00\pm 0.10)·logP\end{array}$$(12)

and a true distance modulus of μ_I = 18.51 ± 0.07 ± 0.05 mag (see the purple line in Fig. 11). Finally, we also adopted the optical Period-Wesenheit (PW Madore 1982) relations for a threefold reason. i) These distance diagnostics are independent, by construction, of the reddening uncertainties. ii) Using some specific combinations of filters, they are minimally affected by metal content (Marconi et al. 2015). iii) They mimic a period-luminosity-colour relation (Madore 1982; Marconi et al. 2015; Neeley et al. 2017). However, they rely on the assumption that the adopted reddening law is universal. We adopted the PW(V,B − I) relation by Marconi et al. (2015) and the optical VBI mean magnitudes provided by Kuehn et al. (2013) and we found μ = 18.52 ± 0.03 ± 0.03 mag. The mean of the homogeneous NIR (PLZ: J,K_s), MIR (PLZ: [3.6], [4.5]) and optical (PLZ: I; PW(V,B − I)) distance determinations gives a mean cluster distance of μ = 18.47 ± 0.02 ± 0.06 mag¹.

Fig. 11. I-, J-, Ks- [3.6]- and [4.5]-band PL relations of Reticulum RRLs. Diamonds display RRab variables, triangles RRc variables and circles the RRd variables. Purple, blue, green, orange, and red symbols display the un-reddened mean magnitudes in the I, J, Ks, [3.6] and [4.5] bands, respectively. The J, Ks, [3.6] and [4.5] magnitudes were artificially shifted by –0.4, –0.8, –1.2 and –1.6 mag for more clarity. The solid lines display the empirical PL relations (Eqs. (8)–(12)). The dashed black lines display the theoretical PLZ relations by Marconi et al. (2015), Neeley et al. (2017), at [Fe/H] = –1.70 (Suntzeff et al. 1992, transformed into the Carretta et al. 2009 metallicity scale) and artificially shifted in magnitude. The PLZI relation was only shifted for the current value of the true distance modulus. The standard deviation of the relations are labelled on the top-left corner.

The current estimates support the evidence that Reticulum belongs to the LMC halo. In particular, the use of optical, NIR, and MIR data suggests that it is located ∼1 kpc closer than the LMC barycentre, although it must be kept in mind that the systematics are of the same order of magnitude of this shift. Distance determinations based on MIR data and on Gaia trigonometric parallaxes suggest that Reticulum might be even closer (∼3 kpc, Muraveva et al. 2018a). More accurate estimates require a novel approach, as suggested by Bono et al. (2019) to simultaneously estimate the cluster mean metallicity, reddening and distance.

7. Summary and final remarks

In this work, we have provided the NIR (JHK_s) light-curve templates of RRab and RRc variables. In the following, we summarize the most interesting results and discuss in more detail some relevant issues.

Homogeneous photometry. Here we publish JHK_s time series, in the 2MASS photometric system, of 254 RRLs in the GGCs ω Cen, M 4 and in the field of the Milky Way. The latter sample was obtained from heterogeneous literature data in four different photometric systems (CIT, SAAO, UKIRT, and ESO) which were homogenized. The overall sample includes both photoelectric and CCD data, collected at telescopes in a wide range of diameter classes (1.3 m–8 m). We provide NIR (JHK_s) characterization (mean magnitudes, light amplitudes, epochs of the mean magnitude on the rising branch) for 94 RRab and 51 RRc variables that were used to generate the light-curve templates.

Light-curve templates. We provide a total of 24 light-curve templates of RRLs: these are divided into Fourier and multi-Gaussian series (PEGASUS) fits of four period bins (one for the RRc and three for the RRab variables) and three photometric bands (J, H, and K_s). The Fourier and PEGASUS series range from the fourth to the seventh order and from the second to the sixth order, respectively. The Fourier templates show residuals with respect to the normalized cumulated light curves used to generate them that are smaller than those corresponding to the PEGASUS templates. However, the latter show fewer secondary, unphysical features (bumps and dips) and their residuals are still smaller than 0.005 normalized mag. We also provide the phases of minimum and maximum light for all the light-curve templates, in order to make it easier for future users to adopt the template even when lacking the epoch of the mean magnitude on the rising branch, which is reported less frequently than the epoch of maximum in large surveys.

Template validation. We have validated our templates and compared our K_s-band templates to those by J96. The tests were performed on both a sub-sample of four RRLs in ω Cen (one per template bin), that were not used to generate the templates, and on a set of 22 Galactic bulge RRLs for which we have VI time series from OGLE and K_s-band time series from the VVV survey. We have checked that, within the dispersion, the mean magnitudes derived by applying the template and the best estimate of the mean magnitude (i.e. the integral over the fit of the light curve, converted into fluxes) are the same. The largest offset is of 0.01 mag (with a standard deviation of 0.04 mag), for the H-band template of short-period RRab variables (RRab1 template bin), which are also the ones with the largest amplitudes, meaning that they are more prone to uncertainties. Compared to our K_s templates, the J96 templates provide results which are similar, showing offsets either comparable or–sometimes–larger than ours.

Reticulum. We have collected literature JK_s time series for 30 over 32 RRLs in the LMC globular cluster Reticulum (Dall’Ora et al. 2004). Using BV time series for the same RRLs (Kuehn et al. 2013), we derived the periods and t_ris to apply our templates and estimated NIR mean magnitudes. We derived new empirical PLJK_s relations, and in turn, new accurate and precise estimates of the distance to Reticulum. We found true distance moduli that agree quite well with each other (μ_J = 18.47 ± 0.10 ± 0.03 mag, μ_Ks = 18.49 ± 0.09 ± 0.05 mag) and with literature values. We adopted homogeneous calibrations for MIR ([3.6], [4.5]) and optical (I) PLZ relations and for the optical PW(V,B − I) relation together with mean magnitudes provided by Muraveva et al. (2018a) and by Kuehn et al. (2013). We find a mean cluster true distance modulus of μ = 18.47 ± 0.02 ± 0.06 mag. According to the most accurate and recent LMC distance determinations (Pietrzyński et al. 2013; Inno et al. 2016), the current μ estimate for Reticulum indicates that this cluster is ∼1 kpc closer to us than the LMC itself.

We briefly outline below some of the most relevant developments of the current project supporting the non-trivial effort for new NIR light-curve templates for RRLs.

Distance Scale. Future ground-based Extremely Large Telescopes (ELT, TMT, GMT) and space observing facilities (JWST, Euclid, WFIRST) have been designed to reach their peak performance in the NIR regime. This means that a few NIR measurements of variables already identified and characterized in the NIR will allow us to fully exploit the RR Lyrae distance scale in Local Group and in Local Volume galaxies. We note that this opportunity fits within a context in which Gaia will provide exquisite calibration for both the zero point and the slope of the diagnostics we are currently using to estimate individual RR Lyrae distances (Gaia Collaboration 2016a,b; Arenou et al. 2018). Moreover, LSST will provide an unprecedented wealth of optical time series, and in turn a complete census of evolved variables in the nearby Universe (Oluseyi et al. 2012). These are crucial prior conditions to reach a precision of the order of 1% on individual RRL distances and an accuracy better than 3% on the Hubble constant (Carnegie RR Lyrae Program, Beaton et al. 2016).

Light curve characterization. Light-curve templates also provide the opportunity to improve the accuracy of the fit of the light curve when either a single or a few measurements are available. We note that this opportunity becomes even more relevant for NIR photometric surveys such as VVV+VVV-X (Minniti et al. 2010), which collect time-series data in the K_s-band and just a few measurements in the J and H bands. Accurate NIR mean magnitudes are, together with optical mean magnitudes, a fundamental ingredient for constraining the distance, the reddening and the metal content of field and cluster RR Lyrae using the recent algorithm (REDIME) suggested by Bono et al. (2019).

Envelope tomography. Our knowledge of linear and non-linear phenomena taking place along the pulsation cycle of a variable star is still limited to a handful of objects. There is solid evidence that moving from the optical to the NIR regime luminosity changes are mainly dominated by variations of radius instead of temperature (Bono et al. 2001; Madore et al. 2013). However, we are still lacking accurate investigations of shock formation and propagation based on NIR spectroscopic diagnostics. The NIR light-curve templates provide the opportunity to trace the colour (V − K) variation along the pulsation cycle, and in turn, the temperature variation. This information is crucial for estimating atmospheric parameters of spectra including a limited number of ionized/neutral heavy element lines (Sollima et al. 2006; Magurno et al., in prep.).

It is a real pleasure to develop a new tool for use by the astronomical community, but it is even more appealing to use it on a broad range of stellar systems. This means the opportunity to provide accurate multi-band (J, H, K_s) individual distances for RRLs in the Magellanic Clouds (OGLE IV, Soszyński et al. 2016) and in nearby dwarf galaxies for which are available both optical and NIR photometry collected with both ACS and WFC3 at HST (Hatt et al. 2017; Monelli et al. 2017). The improved accuracy on individual distances open the paths not only for more quantitative constraints on their 3D shape, but also the opportunity to constrain simultaneously their metallicity and their reddening (Bono et al. 2019).

Acknowledgments

We are grateful to the anonymous referee for the useful suggestions that helped us to improve our paper. We gratefully acknowledge data from the ESO Public Survey program ID 179.B-2002 taken with the VISTA telescope, and products from the Cambridge Astronomical Survey Unit (CASU). Support is provided by the Ministry for the Economy, Development and Tourism, Programa Iniciativa Cientifica Milenio grant IC120009, awarded to the Millennium Institute of Astrophysics (MAS), and by the BASAL Center for Astrophysics and Associated Technologies (CATA) through grant AFB-170002. We would also like to acknowledge the financial support of INAF (Istituto Nazionale di Astrofisica), Osservatorio Astronomico di Roma, ASI (Agenzia Spaziale Italiana) under contract to INAF: ASI 2014-049-R.0 dedicated to SSDC. G.F. has been supported by the Futuro in Ricerca 2013 (grant RBFR13J716). D.M. acknowledges support from FONDECYT Regular grant No. 1170121. M.M. was partially supported by the National Science Foundation under Grant No. AST-1714534.

References

Appendix A: Estimate of the phase of the anchor point (t_ris)

To derive t_ris, we adopted the following approach. For each star we selected one filter for which either the optical (B, V) or the NIR (J) light curve is regular and well sampled. Then, we fit the light curve with either a PLOESS (all the MW RRLs and part of the ω Cen RRLs) or a spline fit (all the M 4 RRLs and the remaining part of the ω Cen RRLs) and derived the mean magnitude from the flux integral of the analytic fit.

We then interpolated the phase at which the rising branch of the fit intersects the mean magnitude (ϕ₀). Finally, t_ris could be obtained as t_i − ((ϕ_i − ϕ₀)⋅P)−P, where t_i and ϕ_i are the epoch and the phase of the ith phase point of the light curve, respectively. This specific phase point was called the “anchor point”. In principle, any of the phase points in the light curve could be an anchor point. However, we have used an interactive procedure, firstly to avoid selecting an anchor point which deviates from the others due to either period changes or phase shifts with time; and secondly to sort the phase points from the most recent to the oldest, to obtain a final estimate of t_ris that is as recent as possible. The latter might seem to be a non-necessary requirement, but it is crucial to have reference epochs which are as close as possible in time to the NIR measurements over which the templates will be applied.

All Tables

All Figures

Fig. 1. Bailey diagram, V-band amplitude versus period, for ω Cen RRLs. Blue circles mark RRc variables, while red squares mark RRab variables. The ranges in period for the RRab light-curve templates are indicated by vertical black dashed lines. The blue striped areas show the ranges in amplitude adopted for the light-curve templates by J96. Note that they provided thresholds in the B band, but here they have been rescaled by 1.25, i.e., the typical amplitude ratio (A(B)/A(V)) for RRab variables (Braga et al. 2016).

In the text

	Fig. 2. Left panels: from left to right the different panels display the Fourier fits of the normalized JHKs light curves. Top to bottom panels: four (RRc, RRab1, RRab2, and RRab3) different ranges in period. Right panels: same as the left, but for the residuals of the normalized light curves with the Fourier fits. The median and standard deviation of the median are labelled in red.
In the text

Fig. 3. Left panels: from left to right the different panels display the PEGASUS fits of the normalized JHKs light curves. From top to bottom the panels show the four (RRc, RRab1, RRab2, and RRab3) different ranges in period. Right panels: same as the left, but for the residuals of the normalized light curves with the PEGASUS fits. The median and standard deviation of the median are labelled in red.

In the text

Fig. 4. Panels a and b: Black crosses show the randomly extracted Ks-band phase points over the light curve of the RRc variable ωCen-V83. Grey dashed lines display the fit of the light-curve template (Fourier, a; PEGASUS, b) to the extracted phase points. The ID of the RRL is labelled. Panels c and d: Same as panels a and b, but for the RRab variable ωCen-V107. The RRab1 light-curve templates were adopted. Panels e and f: Same as panels a and b, but for the RRab variable ωCen-V125. The RRab2 light-curve templates were adopted. Panels g and h: Same as panels a and b, but for the RRab variable ωCen-V15. The RRab3 light-curve templates were adopted.

In the text

	Fig. 5. Same as Fig. 4, but for the H band light-curve templates.
In the text

	Fig. 6. Left (a,d,g,j) and middle (b,e,h,k) panels are the same as Fig. 4, but for the Ks band light-curve templates. The right panels (c,f,i,l) display the fit based on the Ks band the J96 light-curve template.
In the text

Fig. 7. Top-left: Fourier, PEGASUS and J96 template fits applied to RRc variable ωCen-V83. Grey crosses show the randomly extracted phase points. The red crosses display the three phase points of a single random extraction. The horizontal red line shows the mean magnitude of the variable based on the flux average fit of the empirical data. The black dashed curve and the horizontal black long-dashed line show the fit with the Fourier template and its mean magnitude. The magenta dashed curve and the horizontal magenta long-dashed line display the same, but for the PEGASUS fit. The blue dotted curve and the horizontal blue dotted line display the same, but for the J96 fit. Top-right: Same as the top-left, but for the RRab variable ωCen-V107 (RRab1). Bottom-left: Same as the top-left, but for the RRab variable ωCen-V125 (RRab2). Bottom-right: Same as the top-left, but for the RRab variable ωCen-V15 (RRab3). We note that, in these panels, the light curves are not phased using tris, but to an arbitrary epoch (HJD = 2 350 000) to underline that the three phase points method is independent of the reference epoch.

In the text

Fig. 8. Left: Same as Fig. 7, but for random extraction on the light curve of the RRab variable ωCen-V125 (RRab2). In this specific case the three randomly extracted phase points are close in phase. The J96 fit does not provide an accurate estimate of the mean magnitude. Right: Same as the left, but for an extraction where the three randomly extracted phase points have a smaller difference in phase. In this case the light-curve templates do not provide an accurate estimate of the Ks band mean magnitude.

In the text

Fig. 9. Panels a–c: black crosses represent the randomly extracted Ks-band phase points over the light curve of the OGLE RRc variable 11254. Grey dashed lines display the template fits to the individual phase points. The Fourier (panel a), PEGASUS (panel b), and J96 (panel c) RRc light-curve templates are also displayed. The ID of the RRL is labelled. Panels d–f: same as panels a–c, but for the OGLE RRab variable 14578. The RRab1 light-curve templates are displayed. Panels g–i: same as panels a–c, but for the OGLE RRab variable 14806. The RRab2 light-curve templates are displayed. Panels j–l: same as panels a–c, but for the OGLE RRab variable 14958. The RRab3 light-curve templates are displayed.

In the text

Fig. 10. Top: J (left panel) and Ks (right panel) band light curve for the RRc variable V10. The red diamonds display the binned phase points. The black dashed curves show the light-curve template applied to the individual binned phase points. The thick green line displays the mean magnitude listed in Table 11. Middle: same as the top panel, but for the RRab variable V19. Bottom: same as top panel, but for the RRab variable V5.

In the text

Fig. 11. I-, J-, Ks- [3.6]- and [4.5]-band PL relations of Reticulum RRLs. Diamonds display RRab variables, triangles RRc variables and circles the RRd variables. Purple, blue, green, orange, and red symbols display the un-reddened mean magnitudes in the I, J, Ks, [3.6] and [4.5] bands, respectively. The J, Ks, [3.6] and [4.5] magnitudes were artificially shifted by –0.4, –0.8, –1.2 and –1.6 mag for more clarity. The solid lines display the empirical PL relations (Eqs. (8)–(12)). The dashed black lines display the theoretical PLZ relations by Marconi et al. (2015), Neeley et al. (2017), at [Fe/H] = –1.70 (Suntzeff et al. 1992, transformed into the Carretta et al. 2009 metallicity scale) and artificially shifted in magnitude. The PLZI relation was only shifted for the current value of the true distance modulus. The standard deviation of the relations are labelled on the top-left corner.

In the text