New near-infrared period-luminosity-metallicity relations for RR Lyrae stars and the outlook for Gaia

We present results of the analysis of 70 RR Lyrae stars located in the bar of the Large Magellanic Cloud (LMC). Combining spectroscopically determined metallicity of these stars from the literature with precise periods from the OGLE III catalogue and multi-epoch $K_{\rm s}$ photometry from the VISTA survey of the Magellanic Clouds system (VMC), we derive a new near-infrared period-luminosity-metallicity (${\rm PL_{K_{\rm s}}Z}$) relation for RR Lyrae variables. In order to fit the relation we use a fitting method developed specifically for this study. The zero-point of the relation is estimated in two different ways: by assuming the value of the distance to the LMC and by using Hubble Space Telescope (HST) parallaxes of five RR Lyrae stars in the Milky Way (MW). The difference in distance moduli derived by applying these two approaches is $\sim0.2$ mag. To investigate this point further we derive the ${\rm PL_{K_{\rm s}}Z}$ relation based on 23 MW RR Lyrae stars which had been analysed in Baade-Wesselink studies. We compared the derived ${\rm PL_{K_{\rm s}}Z}$ relations for RR Lyrae stars in the MW and LMC. Slopes and zero-points are different, but still consistent within the errors. The shallow slope of the metallicity term is confirmed by both LMC and MW variables. The astrometric space mission Gaia is expected to provide a huge contribution to the determination of the RR Lyrae ${\rm PL_{K_{\rm s}}Z}$ relation, however, calculating an absolute magnitude from the trigonometric parallax of each star and fitting a ${\rm PL_{K_{\rm s}}Z}$ relation directly to period and absolute magnitude leads to biased results. We present a tool to achieve an unbiased solution by modelling the data and inferring the slope and zero-point of the relation via statistical methods.


Introduction
RR Lyrae stars are radially pulsating variables connected with low-mass helium-burning stars on the horizontal branch (HB) of the colour-magnitude diagram (CMD). These objects are Population II stars, which are abundant in globular clusters and in the halos of galaxies. RR Lyrae stars are a perfect tool for studying the age, formation and structure of their parent stellar system. Moreover, they are widely used for the determination of distances in the Milky Way (MW) and to Local Group galaxies. RR Lyrae stars are primary distance indicators because of the existence of a narrow luminosity-metallicity (M V − [Fe/H]) relation in the visual band and of period-luminosity-metallicity (PLZ) relations in the infrared passbands. The near-infrared PL K s Z relation of RR Lyrae stars was originally discovered by Longmore et al. (1986), and later was the subject of study by many different authors (e.g., Bono et al. 2003, Catelan et al. 2004, Del Principe et al. 2006, Sollima et al. 2006, Sollima et al. 2008, Borissova et al. 2009, Coppola et al. 2011, Ripepi et al. 2012a). The near-infrared PL K s Z relation has many advantages in comparison with the visual M V − [Fe/H] relation. First of all, the luminosity of RR Lyrae stars in the K s passband is less dependent on metallicity and interstellar extinction (A K s ∼ 0.1A V ). Furthermore, light curves of RR Lyrae stars in the K s band have smaller amplitudes and are more symmetrical than optical light curves, making the determination of the mean K s magnitudes easier and more precise.
In order to calibrate the PL K s Z relation a large sample of RR Lyrae stars is required, spanning a wide range of metallicities, for which accurate mean K s and [Fe/H] measurements are available. We have selected 70 RR Lyrae variables in the Large Magellanic Cloud (LMC) with spectroscopically determined metallicities in the range of −2.06 < [Fe/H] < −0.63 dex (Gratton et al. 2004). All of them have counterparts in the OGLE III catalogue (Soszyński et al. 2009), therefore very precise periods are available. In order to increase the accuracy of the determination of mean K s magnitudes, multi-epoch photometry is needed. For this reason we are using data from the near-infrared VISTA survey of the Magellanic Clouds System (VMC, Cioni et al. 2011), which is performing K s -band observations of the whole Magellanic System in 12 (or more) epochs, while in many previous studies only single-epoch photometry from the Two Micron All-Sky Survey (2MASS, Cutri et al. 2003) was used. To fit the PL K s Z relation we apply a fitting approach developed for the current study. This method takes into account errors in two dimensions, the intrinsic dispersion of the data and the possibility of inaccuracy in the formal error estimates.
One main issue in the determination of distances with the RR Lyrae PL K s Z relation is the calibration of the zero-point. Trigonometric parallaxes remain the only direct method of determining distances to astronomical sources, free of any assumptions (such as, for instance, the distance to the LMC, etc.) and hence calibrating the PL K s Z zero-point. However, reasonably well estimated parallaxes exist, so far, only for five RR Lyrae variables in the MW observed by Benedict et al. (2011) with the Hubble Space Telescope Fine Guidance Sensor (HST/FGS). In this study we use both a global estimate of the LMC distance and the HST parallaxes in order to calibrate the zero-point of our PL K s Z relation based on LMC RR Lyrae stars. Furthermore, to check whether the RR Lyrae PL K s Z relation is universal and could thus be applied to measure distances in the MW and to other galaxies, we analyse a sample of 23 MW RR Lyrae stars, for which absolute magnitudes in the K and V passbands are available from the Baade-Wesselink studies (e.g., Fernley et al. 1998b, and references therein). Based on these absolute magnitudes and applying our fitting approach we fit the RR Lyrae PL K s Z relation. Then we compare the PL K s Z relations derived for RR Lyrae stars in the MW and in the LMC.
Gaia, the European Space Agency (ESA) cornerstone mission launched in December 2013, is expected to provide a great contribution to the determination of the RR Lyrae PL K s Z relation and to the definition of its zero-point in particular. The satellite is designed to produce the most precise three-dimensional (3D) map of the MW to date (Perryman et al. 2001) by measuring parallaxes of over one billion stars during its five-year mission, among which are thousands of RR Lyrae variables. In the current study we present a method which avoids the problems of the non-linear transformation of trigonometric parallaxes (and negative parallaxes) to absolute magnitudes, and apply this method to fit the PL K s Z relation of the 23 MW RR Lyrae stars, based on simulated Gaia parallaxes.
In Section 2 we provide information about the 70 RR Lyrae stars in the LMC that form the basis of the present study. In Section 3 we present our method and results of fitting the RR Lyrae PL K s Z relation in the LMC and in the MW. In Section 4 we present the method to fit the PL K s Z relation with simulated Gaia parallaxes and apply this method to the 23 MW RR Lyrae stars analysed in Section 3. Section 5 provides a summary of the results. In the Appendix sections we present a detailed description of the fitting method which was developed for this study (Appendix A) and a compilation of metal abundances for the MW RR Lyrae stars (Appendix B).

Data
Optical photometry for the LMC RR Lyrae stars discussed in this paper was obtained by Clementini et al. (2003) and Di Fabrizio et al. (2005) using the Danish 1.54 meter telescope in La Silla, Chile. Two different sky positions, hereafter called fields A and B were observed. Both are located close to the bar of the LMC (Clementini et al. 2003, Di Fabrizio et al. 2005. As a result, accurate B, V and I light curves tied to the Johnson-Cousins standard system and pulsation characteristics (period, epoch of maximum light, amplitudes and mean magnitudes) for 125 RR Lyrae stars were obtained (Di Fabrizio et al. 2005). Low-resolution spectra for 98 of these variables were collected by Gratton et al. (2004) using the FOcal Reducer/low dispersion Spectrograph (FORS1) instrument mounted at the ESO VLT. They were used to derive metal abundances for individual stars by comparing the strength of the Ca II K line with that of the H lines (Preston 1959). For the calibration of the method, four clusters with metallicity in the range [−2.06; −1.26] dex were used. According to Gratton et al. (2004), the obtained metallicities are tied to a scale, which is, on average, 0.06 dex more metal-rich than the Zinn & West (1984) metallicity scale.
We cross-matched the sample of 98 RR Lyrae variables with known metallicities against the catalogue of RR Lyrae stars observed by the OGLE III survey (Soszyński et al. 2009). The OGLE III catalogue contains information about the position, photometric and pulsation properties of 24906 RR Lyrae stars in the LMC. We found that, respectively, 94, 2 and 2 objects are cross-identified with sources in the OGLE III catalogue within a pairing radius of 1 , 3 and 7 . The 2 stars with a counterpart at more than 5 are OGLE-LMC-RRLYR-10345 and OGLE-LMC-RRLYR-10509; for these two objects we checked both the OGLE III finding charts and Gratton et al. (2004) Figure 5 (field B1) in order to understand if they are affected by any problem. Star OGLE-LMC-RRLYR-10345 is an isolated slightly elongated star without any clear blending problem, while star OGLE-LMC-RRLYR-10509 is very close to another source possibly making more difficult to locate accurately the star center. Considering that Gratton et al. (2004) and OGLE III periods for these 2 stars agree within 0.5%, we kept these stars in our sample.
We compared the periods of the 98 RR Lyrae stars provided by Di Fabrizio et al. (2005) and those in the OGLE III catalogue (Soszyński et al. 2009). For 96 objects the periods agree within ∼ 2%, while for two objects periods differ significantly. For star A6332 the difference is of ∼ 25% and for star A5148 it is of ∼ 37% (star identifications are from Di Fabrizio et al. 2005). Moreover, star A5148 has been classified as a first-overtone RR Lyrae star (RRc) in the OGLE III catalogue, and as a fundamental-mode RR Lyrae (RRab) by Di Fabrizio et al. (2005). Since accurately estimated periods and classifications play a key role in the current study, we discarded these two objects from the following analysis.
Seven objects (B2811, B4008, B3625, B2517, A2623, A2119, A10360) from the sample are classified as RRc by Di Fabrizio et al. (2005) and as second-overtone RR Lyrae star (RRe) in the OGLE III catalogue. We removed them from our analysis because of the uncertain classification. Furthermore, since one of the main purposes of the current research is to study the PL K s Z relation of RR Lyrae stars of ab-and c-types we discarded seven objects, which were classified as double-mode RR Lyrae stars (RRd) by Di Fabrizio et al. (2005): A7137, A8654, A3155, A4420, B7467, B6470 and B3347. This left us with a final sample of 61 RRab and 21 RRc stars, which all have a counterpart in the OGLE III catalogue. The period search for the RR Lyrae stars in the OGLE III catalogue was performed using an algorithm based on the Fourier analysis of the light curves (Soszyński et al. 2009). The uncertainties in the OGLE III periods for the 82 RR Lyrae stars in our sample are declared to be less than 5 × 10 −6 days. Therefore we used the periods provided by the OGLE III catalogue in order to fit the PL K s Z relation for our sample, and did not consider errors in the periods since they are negligible in comparison to the other uncertainties.
In order to derive mean K s magnitudes for the RR Lyrae stars in our sample we used data from the VMC survey (Cioni et al. 2011). Started in 2009, the VMC survey covers a total area of 116 deg 2 in the LMC with 68 contiguous tiles. The survey is obtaining YJK s photometry. The K s -band photometry is taken in time-series mode over 12 (or more) separate epochs and each single epoch reaches a limiting K s magnitude ∼ 19.3 mag with a S/N ∼ 5 (see Figure 1 of Moretti et al. 2014). On the bright side, VMC is limited by saturation at K s ∼ 11.4 mag. The majority of RR Lyrae stars in our sample are located within the VMC tile LMC 5_5. Observations of the tile LMC 5_5 were performed in 15 epochs taken in the period from 2010, October 30, to 2012, January 11. For two epochs of observation the ellipticity was too high, so these data were not considered in the analysis. Among the remaining 13 epochs there are 11 deep and 2 shallow epochs. Since shallow observations were obtained in good seeing conditions their S/N was enough to detect the RR Lyrae stars. In the following analysis we used all 13 available epochs to fit the light curves of the RR Lyrae stars. PSF photometry of the time-series data for this tile was performed on the homogenised epoch-tile images (Rubele et al. 2012) using the IRAF Daophot packages (Stetson et al. 1990). On each epoch-tile image the PSF model was created using 2500 stars uniformly distributed, finally the Daophot ALLSTAR routine was used to perform the PSF photometry on all epoch images and time-series catalogues were correlated within a tolerance of one arcsec.
We have cross-matched our sample of 82 RR Lyrae stars against the PSF photometry catalogue of the VMC tile LMC 5_5. VMC counterparts for 71 objects were found within a pairing radius of 1 . Among them, 70 sources have 13 epochs in the K s -band, while for one object (B4749) we have observations only in 6 epochs. Six data points are not enough for a reliable fit of the light curve and, consequently, for the robust determination of the mean K s magnitude, hence, we discarded this source from the following analysis and proceeded with the 70 RR Lyrae stars, for which 13 epochs in the K s -band exist. We derived the mean K s magnitudes of these 70 RR Lyrae stars by Fourier fitting the light curves with the GRaphical Analyzer of TImes Series package (GRATIS, custom software developed at the Observatory of Bologna by P. Montegriffo, see e.g. Clementini et al. 2000). To fit the light curves we discarded obvious outliers. Nevertheless, after the σ-clipping procedure, each source still has 11 or more data points. Examples of the K s light curves are shown in Figure 1.
After deriving K s mean magnitudes we performed the dereddening procedure. Clementini et al. (2003) estimated reddening values of E(B − V) = 0.116 ± 0.017 and 0.086 ± 0.017 mag in LMC field A and B, respectively, using the method from Sturch (1966) and the colours of the edges of the instability strip defined by the RR Lyrae variables. Applying the coefficients from Cardelli et al. (1989) of A K /A V = 0.114 and assuming a ratio of total to selective absorption of R V = 3.1, we estimated the extinction in the K s -band as:  Fabrizio et al. (2005) and in the OGLE III catalogue, respectively. The table also shows coordinates and the classification of the stars from the OGLE III catalogue, metallicities with errors from Gratton et al. (2004) and dereddened mean K s magnitudes, determined with the GRATIS package, along with their errors.  Fabrizio et al. (2005), periods are from the OGLE III catalogue (Soszyński et al. 2009) and are given in days. Using the dereddened mean K s magnitudes of the 70 RR Lyrae stars derived as described in Section 2, spectroscopically determined metallicities from Gratton et al. (2004) and accurately estimated periods from the OGLE III catalogue (with RRc stars "fundamentalized" by adding 0.127 to the logarithm of the period) we can now fit the PL K s Z relation. The fit is performed using a fitting approach developed specifically for this work.
Fitting a line to data is a common exercise in science. Most common approaches use Minimum-Least-Squares methods, however these are often based on assumptions which do not always hold for real observational data. The most basic methods assume that data are drawn from a thin line with errors, which are Gaussian, perfectly known, and exist in one axis only. These assumptions do not hold in the present case, as we have an unknown but potentially significant intrinsic dispersion, non-negligible errors in two dimensions (K s and [Fe/H]), and the possibility of inaccuracy in the formal error estimates (e.g. in the determination of the precision metallicity estimates).
We therefore follow the prescription of Hogg et al. (2010), who develop a method for fitting a line to data which avoids the problems highlighted above by statistical modelling of the data. They present a method for use in two dimensions, which has been extended to three dimensions in this paper.
The method assumes that the data is drawn from a plane of the form Ls(P, [Fe/H] where A is the slope in the logP axis, B is the slope in the metallicity [Fe/H] axis, and C is the intercept. We assume a uniform Gaussian intrinsic dispersion around the luminosity axis, plus the scatter caused by the Gaussian observational errors. The exact mathematical definition is given in Appendix A. The method utilises adaptive Markov Chain Monte Carlo (MCMC) methods (Foreman-Mackey et al. 2013) to evaluate the posterior probability density function (PDF) of each parameter, given an input dataset, and returns the maximum a posteriori probability (MAP) estimate of each parameter, the formal error estimate, and the full posterior PDF. The formal error estimate is obtained from the 16% and 84% quartiles of the posterior PDF of the parameters, which give the 1σ formal error estimate assuming that the posterior PDF is approximately Normal. The free fit parameters are: the slope in logP, the slope in metallicity, the zero-point, and the intrinsic dispersion perpendicular to the magnitude axis.
The intrinsic dispersion of the relation is found to be 0.01 mag. The RMS deviation of the data around the relation, neglecting the intrinsic dispersion, is 0.1 mag. Since the reddening in the K s band is negligible we suggest that the effects of the LMC depth cause the intrinsic dispersion of the relation. The left panel of Figure  It is worth noting that we find a very small dependence of the K s magnitude on metallicity. However, the metal abundance range spanned by the adopted sample does not reach the highest values (up to solar and supersolar) observed in the MW bulge and disk RR Lyrae populations. In order to study the effect of the adopted range of metallicities on the slope of the PL K s Z relation we derived this relation also for MW RR Lyrae stars. We discuss the results in Section 3.3.

Zero-point of the PL K s Z relation in the LMC
To use the derived PL K s Z relation for determining distances it is necessary to calibrate its zero-point. This can be done in a number of different ways. In this paper we follow two different approaches: the first one is based on adopting a value for the distance of the LMC; in the second one we use the absolute magnitudes of Galactic RR Lyrae stars for which trigonometric parallaxes have been measured with the HST/FGS. Both approaches have their advantages and disadvantages, we discuss them in the following sections.

Zero-point based on the LMC distance
The LMC is widely considered the first rung of the cosmic distance ladder as it contains a large number of different distance indicators, such as Cepheids, RR Lyrae variables, eclipsing binaries (EBs), red giant branch (RGB) stars, etc., allowing the galaxy distance to be determined with several independent techniques. Figure 8 of Benedict et al. (2002) shows an impressive summary of LMC distance moduli published during the ten years from 1992 to 2001. Values from 18.1 to 18.8 mag were reported in the literature, with those smaller than 18.5 mag supporting the so-called "short" scale, and those larger than 18.5 mag, the "long" one. In more recent years the dramatic progress in the calibration of the different distance indicators has led the dispersion in LMC distance moduli to shrink significantly. Extreme values such as those listed in Benedict et al. (2002) are not very often seen in the recent literature (Clementini 2008). Still a general consensus on the LMC distance has not been fully reached yet. Moreover, there have been significant concerns about a possible "publication bias" affecting the distance to the LMC (Schaefer 2008, Rubele et al. 2012, Walker 2012. In particular, Schaefer (2008) claimed that from 2002 to 2007 June, 31 independent papers reported new measurements of the distance of the LMC, and the new values clustered tightly around the value (m − M) 0 = 18.5 ± 0.1 mag, adopted by the HST Key Project on the extragalactic distance scale (Freedman et al. 2001). Schaefer (2008) considered the effects of the "publication bias" to be the most likely cause of the clustering of LMC distance measurements.
A number of studies on the compilation of distances to the LMC as derived from different distance indicators can be found in the literature of the last 15 years (e.g., Gibson 2000;Benedict et al. 2002;Clementini et al. 2003;Schaefer 2008;de Grijs et al. 2014). Clementini et al. (2003) analysed the distance to the LMC measured using Population I and Population II standard candles and showed that all distance determinations converge within 1σ error on a distance modulus (m − M) 0 = 18.515 ± 0.085 mag. The most recent compilation of LMC distance moduli is that of de Grijs et al. (2014) Ripepi et al. (2012b) using LMC classical Cepheids observed by the VMC survey, and the recent determination of direct distances to eight long-period EBs in the LMC by Pietrzyński et al. (2013), which is claimed to be accurate to within ∼ 2%: D LMC = 49.97 ±0.19 (stat) ± 1.11 (syst) kpc, corresponding to the distance modulus (m − M) 0 = 18.493 ± 0.008 (stat) ± 0.047 (syst) mag. Furthermore, the model fitting of the light curves of different classes of pulsating stars in the LMC, also based on different samples and hydrodynamical codes, provides values consistent with 18.5 mag (see Bono et al. 2002;Marconi & Clementini 2005;Keller & Wood 2002McNamara et al. 2007).
The RR Lyrae stars in our sample are located in a relatively small area close to the center of the LMC bar. Neglecting depth/projection effects they can be considered as being all at the same distance from us and close to late-type EBs, which are all located relatively close to the barycenter of the LMC as derived by Pietrzyński et al. (2013). Therefore, in the following we adopt for the distance modulus of the LMC the value published by Pietrzyński et al. (2013). We subtracted this value from the dereddened mean K s apparent magnitudes of our 70 RR Lyrae stars and derived absolute magnitudes in the K s band (M K s ). Then by applying the technique described in Section 3.1 we derived the relation between K s -band absolute magnitudes, periods and metallicities, with the zero-point entirely based on the distance to the LMC by Pietrzyński et al. (2013):

Zero-point based on trigonometric parallaxes of Galactic RR Lyrae stars
In order to obtain an estimate of the PL K s Z relation zero-point which is independent of the distance to the LMC and, in turn, be able to measure the distance to this galaxy from the RR Lyrae PL K s Z relation, it is necessary to know the absolute magnitude of the RR Lyrae stars with good accuracy. Trigonometric parallaxes remain the only direct method to measure distances and hence derive absolute magnitudes. Benedict et al. (2011) derived absolute trigonometric parallaxes for five Galactic RR Lyrae stars (RZ Cep, XZ Cyg, SU Dra, RR Lyr and UV Oct) with the HST/FGS. With these parallaxes the authors estimated absolute magnitudes in the K s and V passbands, corrected for interstellar extinction and Lutz-Kelker-Hanson bias (Lutz & Kelker 1973, Hanson 1979. Absolute magnitudes in the K s -band, periods and metallicities from Benedict et al. (2011), and the slopes of the relation derived in Equation 3 were used in order to determine a zero-point from each of these five MW RR Lyrae stars. The metallicities in Benedict et al. (2011) are in the Zinn & West metallicity scale and were converted to the metallicity scale in Gratton et al. (2004) by adding 0.06 dex (see Section 2). The logarithm of the period of the RRc star RZ Cep was "fundamentalized" by adding 0.127. Then we calculated the weighted mean of the five zero-points, this corresponds to: −1.27 ± 0.08 mag. The relation between absolute magnitudes, periods and metallicities with the zero-point based on the five MW RR Lyrae stars from Benedict et al. (2011) is: Situation improves, however there is still a difference of ∼ 0.2 mag between the two zero-point obtained based on the distance to the LMC (Eq. 4) and the one based on the HST parallaxes of four RR Lyrae stars (Eq. 6). In fact, if we apply our PL K s Z relation with zero-point calibrated on Benedict et al. (2011) parallaxes (Eq. 6) to determine the absolute magnitudes of the 70 RR Lyrae stars in our sample, we obtain a distance modulus for the LMC of (m − M) 0 = 18.68 ± 0.10 mag. This distance modulus is about 0.2 mag longer than the widely adopted value of (m − M) 0 = 18.5 mag.
There are a number of possible explanations for the discrepancy between zero-points. First of all, we should remember that Pietrzyński et al. (2013) results have been called into question by Schaefer (2013) who, in addition to concerns regarding possible bandwagon effects, also pointed out that Pietrzyński et al. (2013) distance to the LMC differs significantly from the average distance inferred from four hot, early-type EBs, D=47.1±1.4 kpc ((m − M) 0 = 18.365 ± 0.065 mag), published by Guinan et al. (1998, and Ribas et al. (2002). Furthermore, in using the late-type EBs to calibrate the RR Lyrae PL K s Z relation we have implicitly assumed that RR Lyrae stars and EBs are at same distance from us. However, when pushing for distance comparisons at a few percent level the effects of sample size, spatial distribution, depth and geometric projection become important and properly accounting for the internal structure of the LMC may become necessary. The RR Lyrae stars in our sample could be distributed along the whole depth of the LMC. Moreover, RR Lyrae stars and EBs from Pietrzyński et al. (2013) could reside in different sub-structures of the LMC, which could be the reason for the systematic error in the determination of the zero-point (see e.g. Moretti et al. 2014 for different features of the LMC structure traced by classical Cepheids, RR Lyrae stars and hot EBs).
On the other hand, when calibrating the zero-point by applying parallaxes of the four MW RR Lyrae stars by Benedict et al. (2011) we implicitly assumed that the PL K s Z relation is the same in the MW and in the LMC, which may not be true (see Subsection 3.3). We may also wonder whether there might be unknown systematic errors affecting Benedict et al.'s parallaxes. These come from HST fields, which provide relative and not absolute trigonometric parallaxes. Absolute parallaxes of the reference stars in each field are estimated via a complex procedure of fitting the spectral type and luminosity class of each star. A general formal error of 0.5 mas is applied to the absolute parallax of the reference stars, equal for all stars in all fields, and without justification. This could result in miscalculated estimates of the precision of the final absolute parallax measurements of the four RR Lyrae stars. The Lutz-Kelker bias is corrected a posteriori. In this respect it is worth of notice that, according to van Leeuwen (2007), Hipparcos parallax of RR Lyrae itself, the only RR Lyrae variable for which the satellite measured the parallax with high precision (± 0.64 mas), is about 0.31 mas smaller than Benedict et al. (2011)'s parallax for the star, although consistent with it within the errors, hence, the corresponding distance modulus is about 0.17 mag longer. In any case, a great contribution to the determination of the zero-point of the RR Lyrae PL K s Z relation is expected from the ESA astrometric satellite Gaia. We discuss this topic in Section 4.

PL K s Z relation of RR Lyrae stars in the MW
In spite of many studies in the literature, it remains still unsettled whether the RR Lyrae PL K s Z is a universal relation. To investigate this matter we have derived the PL K s Z relation for RR Lyrae stars in the MW and compared it with the relation we have obtained in Section 3.2 for the LMC variables. To this end we selected 23 MW RR Lyrae stars which have their absolute magnitudes known from Baade-Wesselink (hereinafter B-W) studies based on near-infrared data (Jones et al. 1988(Jones et al. , 1992Fernley et al. 1990;Liu & Janes 1990;Cacciari et al. 1992;Skillen et al. 1993;Fernley et al. 1994, and references therein) and metallicities determined from abundance analysis of high resolution spectra. Information about these 23 RR Lyrae stars is presented in Table 2. Star's coordinates in the table are from the SIMBAD database; periods, apparent V and K magnitudes and reddening E(B-V) are from Fernley et al. (1998a). The sample contains two first-overtone mode RR Lyrae stars (namely, T Sex and TV Boo). As done for the LMC RRc stars, their periods were fundamentalized by adding +0.127 to the logarithm of the period. Absolute M V and M K magnitudes in Columns 10 and 12 were taken from the compilations of B-W results in Table 11 of Cacciari et al. (1992) and from Table 16 of Skillen et al. (1993) for the variable stars: WY Ant, W Crt and RV Oct. According to Cacciari et al. (1992) the K magnitudes of the stars analyzed with the B-W method are in the Johnson photometric system. Following the discussion in Cacciari et al. (1992) and Skillen et al. (1993) we retained only 23 of the original lists of 29 stars field RR Lyrae stars analysed with the B-W method, as stars which are likely to be evolved (DX Del, SU Dra, SS Leo, BB Pup and W Tuc) were discarded. We also discarded DH Peg as there is suspect the star is a dwarf Cepheid (see Feast et al. 2008 and discussion therein). Furthermore, following Fernley et al. (1994), original M V and M K values were revised (i) assuming for the p factor used to convert the observed pulsation velocity to true pulsation velocity in B-W analyses the value of p=1.38, and (ii) multiple determinations of individual stars were averaged.
Metal abundances with related errors are needed to apply our fitting approach. Several different spectroscopic studies have targeted the stars in Table 2. In Appendix B we provide a summary of their results. The largest and most homogeneous samples are those by Clementini et al. (1995) and Lambert et al. (1996). These Authors measured [Fe/H] abundances from high resolution spectra for several of the stars in Table 2 and provided recalibrations of the ∆S index (Preston 1959), from which metal abundances can be derived for the stars which lack abundance analysis. For sake of homogeneity and ease of use in this paper we adopt metallicities and metallicities errors for the MW RR Lyrae stars as they are listed, ready for use, in Table 21 of Clementini et al. (1995). These [Fe/H] values are the average of the FeI and FeII abundances, adopting log (FeI)= 7.56 and log (FeII)=7.50 for the sun. They are reported along with the related errors in columns 5 and 6 of Table 2. Metallicities in Lambert et al. (1996) were obtained from the FeII abundance and adopting log (FeI)= 7.51 for the sun. Corresponding [Fe/H] values and errors, are reported in columns 2 and 3 of Table 4 in Appendix B. In Appendix B we also present the PLZ relation obtained using Lambert et al. (1996)'s metallicities and our approach.
By applying our fitting approach (Section 3.1) to the 23 MW RR Lyrae stars we derived the following PL K Z relation: The intrinsic dispersion of the relation is found to be 0.007 mag. The RMS deviation of the data around the relation, neglecting the intrinsic dispersion, is 0.086 mag. It should be noted that the metallicities listed in Table 2 may differ slightly from the metal abundances used in the B-W analysis of these stars. However, this is not of great concern as the B-W based on near-infrared data is mildly affected by small changes in metallicity and reddening. We also point out that the rather large error of the logP term in Equation 7 is largely driven by the large errors (0.15-0.25 mag) in the K-band absolute magnitudes from the B-W analyses. This is confirmed by the exercise with Gaia simulated parallaxes we present in Section 4.
The slope in [Fe/H] in Equation 7 is higher than the slope obtained for the LMC RR Lyrae stars (Equations 4, 6), although they are still consistent within the respective errors. Equation 7 was derived over a wide range of metallicities [-2.5; 0.17] dex, nevertheless the slope of the metallicity term remains rather small. Thus, the relatively small metallicity range spanned by the LMC variables could be not responsible for the negligible dependence on metallicity of the RR Lyrae PL K s Z relation in the LMC.
The distribution of the 23 MW RR Lyrae stars in the period-luminosity-metallicity space and the projections of the PL K Z relation (Equation 7) on the log(P)− M K and M K −[Fe/H] planes are shown in Figure 3. Grey lines are the same as in Figure 2 and are described in Section 3.1.
Some concern may arise since the K magnitudes of the 23 MW RR Lyrae stars used to derive Equation 7 are in the Johnson photometric system (see, Cacciari et al. 1992), whereas for the LMC RR Lyrae variables we have K s photometry in the VISTA system 2 . To address this issue we have reported in column 10 of Table 2 the average K s magnitudes in the 2MASS system of the 23 MW RR Lyrae stars as provided by Feast et al. (2008). The difference with the Johnson K average magnitudes listed in column 9 is very small (of the order of about 0.03 mag, on average) and definitely much smaller than individual errors in the B-W K-band absolute magnitudes of the MW variables (0.15-0.25 mag), or errors in the K s average apparent magnitudes of the LMC RR Lyrae stars (see column 10 of Table 1). Hence, we are confident that the difference in photometric system does not affect significantly our comparison.

Comparison with the literature
The near-infrared PL K Z relation of the RR Lyrae stars has been studied by several authors both from a theoretical and an observational point of view. Longmore et al. (1986) pioneering work was followed by Liu & Janes (1990), Skillen et al. (1993) and Jones et al. (1996). A comprehensive analysis of the IR properties of RR Lyrae stars was performed by Nemec et al. (1994).
Some of the RR Lyrae PL K Z relations available in the literature are presented in Table 3. Bono et al. (2003) derived the semi-theoretical relation presented in the first row of Table 3. This theoretical relation has been derived from an extended set of RR Lyrae nonlinear hydrodynamical models spanning a wide range of chemical compositions (Z from 0.0001 to 0.02, which approximately corresponds to [Fe/H] from -2.45 to -0.15 dex). Catelan et al. (2004) presented a theoretical calibration of the RR Lyrae PL K Z relation based on synthetic HB models computed for several different metallicities, fully taking into account evolutionary effects besides the effect of chemical composition. They derived the relation: By using Eqs. 9 and 10 in Catelan et al. (2004) and assuming [α/Fe]∼ 0.3 (e.g., Carney 1996) we transformed Eq. 8 into the form, presented in the second row of Table 3.  Table 1.  Table 2.  Note. -The columns report: 1) Name of the star; 2) Right Ascension (J2000) from SIMBAD database; 3) Declination (J2000) from SIMBAD database; 4) Period from Fernley et al. (1998a); 5) Metallicity from Clementini et al. (1995); 6) Error in metallicity from Clementini et al. (1995); 7) Reddening from Fernley et al. (1998a); 8) V magnitude from Fernley et al. (1998a); 9) K magnitude in the Johnson system from Fernley et al. (1998a); 10) K s magnitude in the 2MASS system from Feast et al. (2008); 11) Absolute magnitude in the V passband from Fernley et al. (1994); 12) Absolute magnitude in the K passband obtained from B-W analyses and corrected to p=1.38 (see text for details); 13) Errors in the absolute V, K magnitudes. Dall'Ora et al. (2004) derived an empirical relation between apparent K magnitude and period for 21 RRab and 9 RRc stars in the LMC globular cluster Reticulum. Del Principe et al. (2006) obtained the relation between apparent K s magnitude, metallicity and period from the analysis of RR Lyrae stars in the Galactic globular cluster ω Cen. Sollima et al. (2006) derived a PL K Z relation by analysing 538 RR Lyrae stars in 15 Galactic clusters and in the LMC globular cluster Reticulum. This relation spans the metallicity range −2.15 < [Fe/H] < −0.9 dex. Mean K magnitudes were estimated by combining Two-Micron-All-Sky-Survey (2MASS, Cutri et al. 2003) photometry and literature data. The zero-point was calibrated on RR Lyrae itself, whose distance modulus was derived using the HST trigonometric parallax measured for this star by Benedict et al. (2002). Sollima et al. (2008) presented JKH time-series photometry of RR Lyrae and derived a new zero-point of Sollima et al. (2006)'s PL K Z relation. Borissova et al. (2009) presented near-infrared K s photometry and spectroscopically measured metallicity for a sample of 50 field RR Lyrae stars in inner regions of the LMC. These authors had 5 measurements in the K s passband for most of the stars in their sample and used templates from Jones et al. (1996) to fit the light curves and derive the mean K s magnitudes. To improve statistics they added to their sample LMC RR Lyrae stars from Szewczyk et al. (2008) Table 3, whereas the zero-point of our LMC PL K s Z relation calculated by assuming the distance to the LMC (Subsection 3.2) and the zero-point of our MW PL K Z relation based on the B-W studies are (Section 3.3) are presented in Column 4.
The slope in period of the RR Lyrae PL K Z relations (Column 2 of Table 3) differs significantly in the various studies. The slope we derived for the LMC RR Lyrae stars is in excellent agreement with that derived by Del Principe et al. (2006), whereas the slope of the MW RR Lyrae PL K Z is in good agreement with that derived by Sollima et al. (2006Sollima et al. ( , 2008. The dependence on metallicity of the PL K Z relations (Column 3 of Table 3) also varies among different studies and generally is larger in the theoretical and semi-theoretical relations. The comparison of the metallicity dependence in the different empirical relations is complicated by the inhomogeneity of the metallicity scales adopted in these studies. Metallicities in Del Principe et al. (2006) are in the Zinn & West (1984) scale, whereas in Sollima et al. (2006Sollima et al. ( , 2008 are in the Carretta & Gratton (1997) scale. In the current study for the LMC RR Lyrae stars we used the metallicity scale defined by Gratton et al. (2004) which is also the scale adopted by Borissova et al. (2009). As discussed in Gratton et al. (2004) this scale is systematically 0.06 dex higher than Zinn & West scale. This difference is small and systematic, hence should not affect the results of this comparison. Finally, for the MW RR Lyrae stars we used the metallicities measured by Clementini et al. (1995). Because the spectroscopic [Fe/H] values in Clementini et al. (1995) are derived from high dispersion spectra analyzed using standard reduction procedures, the derived metallicities are on the scale of the high dispersion spectra (i.e., the Carretta et al. 2009 scale) and could be transformed to Zinn & West scale using the relations provided in Carretta et al. (2009).
The slope in metallicity of the PL K Z relation based on the LMC RR Lyrae stars is the smallest among the various studies listed in Table 3 and it is close to Borissova et al. (2009)'s slope. This is consistent with the two studies both involving LMC variables and using exactly the same metallicity scale. The slope in metallicity we found for the MW RR Lyrae stars is larger than that of the LMC RR Lyrae stars and, in spite of the difference in metallicity scales, it is very close to the slope obtained by Sollima et al. (2006) for RR Lyrae stars in globular clusters. However, taken at face value, the metallicity slopes of the empirical relations in Table 3 appear to be all rather small and in agreement to each other within the relative uncertainties, thus generally suggesting a mild dependence the RR Lyrae PL K Z independently on the specific environment.

Gaia observation of RR Lyrae stars in the Milky Way
The Gaia astrometric satellite will revolutionise many fields of astronomy (Perryman et al. 2001). Of particular importance will be its catalogue of trigonometric parallaxes for more than one billion stars, with astrometric precision down to µas level. Due to Gaia's constant observation of the sky over the five year nominal mission, Gaia will repeatedly observe all stars brighter than its limiting magnitude, with an average of 70 observations per star. This will also make it possible for Gaia to discover and characterise many types of variables, including RR Lyrae stars and Cepheids.
Gaia is observing in the broad visual band G (Jordi et al. 2010) for its astrometric measurements, and is therefore not ideal for characterising the RR Lyrae PLZ relation, which exists only in the infrared passbands. However, since Gaia will provide accurate parallaxes for an expected tens of thousands of MW RR Lyrae stars, it could serve as a perfect tool for the determination of the zero-point of the PL K s Z relation through a combination with external datasets. As it was discussed in Subsection 3.2.2, the current largest limiting factor in zero-point calibration of PL K s Z and M V − [Fe/H] relations is the lack of a reliable and statistically significant sample of parallax measurements. The current state of the art is the sample of five RR Lyrae parallaxes from Benedict et al. (2011) using the HST. Gaia will improve this situation by several orders of magnitude in both precision and numbers of objects. Moreover, the distance to the LMC will be determined through the combination of Gaia parallaxes for a large sample of LMC bright stars, hence, a zero-point of the PL K s Z relation based on the distance to the LMC will be obtained with a high precision.

Method
Using parallax data for calibration of a PL relation is complicated by the presence of statistical biases (e.g. Lutz & Kelker 1973) and sample selection effects (e.g. Malmquist 1936). Non-linear transformations on the parallax cause a highly asymmetric uncertainty on the absolute magnitude when calculated using parallax and apparent magnitude information via the relation: m − M =−5 − 5 log ( ), where is the parallax. Additionally, stars with a negative parallax measurement can not be used to calculate an absolute magnitude, though they do contain information. For these reasons, calculating an absolute magnitude for each star and fitting a PL relation directly to period and absolute magnitude leads to a biased result.
An unbiased solution can be achieved through modelling the data and inferring the slope and zero-point of the relation via statistical methods. For a catalogue of N stars we can define x = (x 1 , x 2 , ..., x N ) where the vector (x i = m, l, b, P, , A) describes the observed data on each object. P is the period, m the apparent magnitude, l and b the position, the parallax, and A the extinction). We can additionally define that the vector (x 0 = m 0 , l 0 , b 0 , P 0 , r 0 , A 0 ) gives the 'true' underlying object properties.
We assume that the stars follow a PL relation of the form M 0 = ρlogP 0 +δ, although this can be changed to include other terms, such as metallicity, as needed. We can therefore model the true absolute magnitudes of the population as being Normally distributed around the PL relation, with the dispersion describing the intrinsic scatter on the relation: where σ PL is the intrinsic dispersion of the PL relation. The parameters ρ and δ are the slope and zero-point of the PL relation, which are to be found.
The true absolute magnitude is calculated through: The observations are Normally distributed around the true values with a standard deviation given by the formal error on the measurement: Assuming negligible errors on the position and period, the observations are described by a delta function. The terms , m , A are the formal errors on the parallax, magnitude and extinction.
With the above models defined, the joint probability density function for the observations is: the 'true' parameters x 0 are never known and so these values are marginalised through integration. The term S(x) is the selection function, which takes the probability of observing a star into account, given the properties of the star and the instrument's observational capabilities. To take the fact that Gaia is a magnitude limited sample into account, a step function is used with The Maximum Likelihood Estimation of the parameters are found by maximising equation 12 by varying the parameters (α, ρ, σ PL ). This formulation avoids non-linear transformations on error effected data, and includes a selection function which avoids the Malmquist bias.

Simulated Gaia data
In order to check the application of the method defined in Section 4.1 we have used the sample of 23 RR Lyrae stars in the MW discussed in Section 3.3 (see also Table 2). In order to investigate the performance of the Gaia satellite and the contents of the end-of-mission catalogue, Gaia's Data Processing and Analysis Consortium (DPAC) has a group working on the simulation of several aspects of the Gaia mission. One major product of this work is the Gaia Object Generator (GOG; Luri et al 2014), designed to simulate both individual Gaia observations and the full contents of the end-of-mission catalogue. GOG includes a full mathematical description of the nominal performance of the Gaia satellite, and is therefore capable of determining the expected precision in astrometric, photometric and spectroscopic observations. In general, the precision depends on the apparent magnitude of the star, its colour, and its sky position, which affects the number and type of observations made (due to the Gaia scanning law).
To obtain a distance for each RR Lyrae star from the sample, we use: as determined in Equation 7 to obtain an absolute magnitude (neglecting the metallicity term for simplicity). We then determine a distance by combining this absolute magnitude with the apparent magnitude and extinction as defined above. Colour information as (V-I) is obtained from the Hipparcos catalogue (Perryman and ESA 1997) where available. The apparent magnitude, position, colour, and period data form the basis of a synthetic catalogue of RR Lyrae stars, along with the distance obtained from the PM K Z relation, and is used as the input catalogue of 'true' parameters for GOG.
GOG then creates simulated Gaia observations for our sample. We take the PM K Z elation (Equation 14) as true, as a study of the possible precision in PM K Z calibration after the Gaia data will become available.
Using the fitting method described in Section 4.1 to the data including the simulated parallax observations and simulated errors applied to parallax and apparent magnitude, we find a PM K relation of: Comparison of these results to the input relation shows very good agreement. It proves that the fitting procedure given in Section 4.1 is accurate and unbiased. When Gaia parallaxes will become available for the much larger sample of RR Lyrae stars, we will apply the described method to fit the PL K Z relation of RR Lyrae variables in the MW. Moreover, precise distance to the LMC obtained from the combination of Gaia parallaxes for a large sample of the bright LMC stars, will allow us to calibrate zero-point of the PL K s Z relation based on the LMC RR Lyrae stars. This provides a flavor of what will be possible to achieve with Gaia parallaxes.

Summary
We studied a sample of 70 RR Lyrae stars in the LMC, for which multi-epoch K s photometry from the VMC survey, precise periods from the OGLE III catalogue and spectroscopically determined metallicities (Gratton et al. 2004), are available. There are 13 epoch data in the K s band for all stars in the sample, that allowed us to determine mean K s magnitudes with a great accuracy.
Specifically for this work we developed a fitting approach. This method has several advantages compared to the Minimum-Least-Squares fitting, such as taking into account potentially significant intrinsic dispersion of the data, non-negligible errors in two dimensions and the possibility of inaccuracies in the formal error estimates. We used this method to derive the PL K s Z relation of the 70 RR Lyrae stars in the LMC. Potentially the method could be used to fit any other sample of data.
The zero-point of the derived PL K s Z relation was estimated in two different ways: (i) by assuming the distance to the LMC determined by Pietrzyński et al. (2013); (ii) by applying HST parallaxes of four MW RR Lyrae stars by Benedict et al. (2011). The zero-point derived using the MW RR Lyrae stars is 0.2 mag larger and, consequently, gives a longer distance to the LMC: (m − M) = 18.68 ± 0.10 mag. In future studies we suggest to use the relation based on the precise distance to the LMC: We found a negligible dependence of the M K s on metallicity, which could be caused by the relatively small range in metallicity covered by the LMC RR Lyrae stars. Thus, we applied the fitting approach to 23 RR Lyrae stars in the MW, for which absolute M K and M V magnitudes are known from Baade-Wesselink studies. We derived the PL K s Z relation for MW RR Lyrae stars in the form: M K = (−2.53 ± 0.36)logP + (0.07 ± 0.04) [Fe/H] − (0.95 ± 0.14) Even though the metallicities of the MW RR Lyrae stars span a wide range [-2.5; 0.17] dex, the dependence on metallicity is relatively small and consistent, within the errors, with the slope in metallicity found for the LMC RR Lyrae variables. We concluded that the small range of metallicities doesn't cause the negligible dependence of the M K on metallicity for the LMC RR Lyrae stars.
To solve the problem of the PL K s Z zero-point, a large sample of RR Lyrae stars with precisely determined parallaxes is necessary. A great contribution to this field is expected by the Gaia satellite. By using GOG we simulated Gaia parallaxes of 23 MW RR Lyrae stars with observational errors. We present a method for the calibration of the PL relation which avoids several of the problems which arise when using parallax data. The method was tested by deriving the PL K s relation based on the simulated Gaia parallaxes. When combined with metallicity and photometry from other sources and a statistical tool such as the one developed in the present study, the extraordinary large sample of Gaia parallaxes for RR Lyrae stars will allow us to estimate these relations with unprecedented precision.
Therefore the likelihood is defined as: Taking the logarithm, which is effectively the least-squares solution. K is a normalisation coefficient. Returning to Bayes rule it is possible to define: , plus any other prior information which may be available. In our case we use uninformative (uniform) priors, making our inference method analogous to Maximum Likelihood Estimation.

A.1. Multiple errors, no dispersion
As in this case there exist errors in more than one axis, they can be put together into a covariance tensor With errors in several dimensions, our observed data point (x i ,y i ,z i ) could have been drawn from any true point along the plane (x,y,z). Making the probability of the data, given the model and the true position: where we have implicitly made column vectors In only two dimensions (e.g. y and z), the slope (e.g. B) can be described by a unit vectorû orthogonal to the line or linear relation (at any x):û where the angle θ = arctan B is made between the line and the y axis. The orthogonal displacement ∆ i of each data point (y i , z i ) from the line is given by Instead of extending fully into three dimensions, we will assume a negligible error in x (which will be the period, so has justifiably higher precision). The value of x can therefore be input directly into ∆ i without worrying about the interplay between the other parameters.
Assuming negligable errors in x also redefines the covariance matrix of the errors as: Similarly, each data point's covariance matrix S i projects down to an orthogonal variance Σ 2 i given by and then the log likelihood for (A, B, C) or (A, θ, C cos θ) can be written as where K is some constant. This likelihood can be maximized to find A, B and C.

A.2. Dispersion
The final step is to introduce an intrinsic variance in the line, V, orthogonal to the line.
According to Hogg et al. (2010), each data point can be treated as being drawn from a projected distribution function that is a convolution of the projected uncertainty Gaussian, of variance Σ 2 i defined above, with the Gaussian intrinsic scatter of variance V. Therefore the likelihood becomes: where again K is a constant, everything else is defined as above. We then solve for A, B, C, and V by maximising the log likelihood. The optimisation is performed using the adaptive MCMC sampler EMCEE developed by Foreman-Mackey et al. (2013). Any optimisation algorithm (e.g. Nelder-Mead, Powell, etc.) will find the maximum of the log likelihood. MCMC was chosen due to the evaluation of the full posterior PDF of the parameters, which is useful for the determination of formal errors.

B. Metallicities for the MW RR Lyrae stars analyzed with the B-W method
In this Appendix we provide a summary of spectroscopic metal abundances ([Fe/H]) derived for the field RR Lyrae stars analyzed with the B-W technique, by studies mainly based on high resolution spectroscopic material. Exceptions are the values with reference 1, 2 and 4 that come from compilations which include metallicities measured from low resolution spectroscopic data and photometric indices (see discussion in Section 4.2 of Cacciari et al. 1992).