NEW NEAR-INFRARED PERIOD–LUMINOSITY–METALLICITY RELATIONS FOR RR LYRAE STARS AND THE OUTLOOK FOR GAIA^*

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 catalog (with RRc stars "fundamentalized" by adding 0.127 to the logarithm of the period) we can now fit the ${\mathrm{PL}}_{{K}_{{\rm{s}}}}{\rm{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 that 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 that avoids the problems highlighted above by statistical modeling 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

$$\begin{eqnarray}&&{Ls}(P,[\mathrm{Fe}/{\rm{H}}])=A\;\mathrm{log}\;P+B\;[\mathrm{Fe}/{\rm{H}}]+C,\quad \end{eqnarray} \tag{ 2 }$$

where A is the slope in the $\mathrm{log}\;P$ 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 utilizes 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 data set, and returns the maximum a posteriori probability 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 log P, the slope in metallicity, the zero-point, and the intrinsic dispersion perpendicular to the magnitude axis.

By applying this method we found the following relation between period, metallicity and mean apparent K_s magnitude:

$$\begin{eqnarray}{K}_{{\rm{s}},0} & = & (-2.73\pm 0.25)\mathrm{log}\;P+(0.03\pm 0.07)[\mathrm{Fe}/{\rm{H}}]\\ & & +(17.43\pm 0.01).\end{eqnarray} \tag{ 3 }$$

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 2 presents the ${\mathrm{PL}}_{{K}_{{\rm{s}}}}{\rm{Z}}$ relation (Equation (3)) of the 70 LMC RR Lyrae stars in the period–luminosity–metallicity space, whereas the right panels show the projection of the ${\mathrm{PL}}_{{K}_{{\rm{s}}}}{\rm{Z}}$ on the $\mathrm{log}(P)-{K}_{{\rm{s}}}$ (top-right panel) and ${K}_{{\rm{s}}}-[\mathrm{Fe}/{\rm{H}}]$ (bottom right panel) planes. The gray lines in the figure are lines of equal metallicity (top-right) or equal period (bottom right). The method finds the relation (values of A, B, and C for the relation ${K}_{{\rm{s}}}=A\;\mathrm{log}\;P+B\;[\mathrm{Fe}/{\rm{H}}]+C)$ in the three dimensions (log P, K_s, and [Fe/H]). Each of the gray lines in the top-right plot are therefore ${K}_{{\rm{s}}}=A\;\mathrm{log}\;P+B\;[\mathrm{Fe}/{\rm{H}}]+C$ for the full range of periods, at the metallicity of each star (one line per star). Thus, by following the line up and down it is seen how K_s changes with period at some specific metallicity. The lines do not always cross the points on the diagram because the line is the result of the fit, and the points are affected by errors and intrinsic dispersion so they may be above or below the fit. In the bottom right plot the lines are ${K}_{{\rm{s}}}=A\;\mathrm{log}\;P+B\;[\mathrm{Fe}/{\rm{H}}]+C$ for the full range of metallicity with log P taken from each star.

Zoom In Zoom Out Reset image size

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 ${\mathrm{PL}}_{{K}_{{\rm{s}}}}{\rm{Z}}$ relation we derived this relation also for MW RR Lyrae stars. We discuss the results in Section 3.3.

To use the derived ${\mathrm{PL}}_{{K}_{{\rm{s}}}}{\rm{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.

3.2.1. Zero-point Based on the LMC Distance

The LMC is widely considered to be 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 stars, etc., allowing the galaxy distance to be determined using 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\pm 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) analyzed 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\pm 0.085$ mag. The most recent compilation of LMC distance moduli is that of de Grijs et al. (2014) who compiled 233 separate distance determinations, published from 1990 March until 2013 December, and concluded that the canonical distance modulus of ${(m-M)}_{0}=18.49\pm 0.09$ mag may be used for all practical purposes. The compilation of de Grijs et al. (2014) includes the distance modulus of ${(m-M)}_{0}=18.46\pm 0.03$ estimated by 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; Keller & Wood 2002, 2006; Marconi & Clementini 2005; McNamara 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}_{{\rm{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):

$$\begin{eqnarray}{M}_{{K}_{{\rm{s}}}} & = & (-2.73\pm 0.25)\mathrm{log}\;P+(0.03\pm 0.07)[\mathrm{Fe}/{\rm{H}}]\\ & & -(1.06\pm 0.01).\end{eqnarray} \tag{ 4 }$$

3.2.2. Zero-point Based on Trigonometric Parallaxes of Galactic RR Lyrae Stars

In order to obtain an estimate of the ${\mathrm{PL}}_{{K}_{{\rm{s}}}}{\rm{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 ${\mathrm{PL}}_{{K}_{{\rm{s}}}}{\rm{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:

$$\begin{eqnarray}{M}_{{K}_{{\rm{s}}}} & = & (-2.73\pm 0.25)\mathrm{log}\;P+(0.03\pm 0.07)[\mathrm{Fe}/{\rm{H}}]\\ & & -(1.27\pm 0.08).\end{eqnarray} \tag{ 5 }$$

A recent analysis (A. Monson 2015, private communication) shows that there is likely a typo in Benedict et al.'s (2011) parallax for the RR Lyrae star RZ Cep. Hence, we excluded this star from the sample and derived the zero-point based on parallaxes of the remaining four RR Lyrae stars (XZ Cyg, UV Oct, SU Dra, and RR Lyr):

$$\begin{eqnarray}{M}_{{K}_{{\rm{s}}}} & = & (-2.73\pm 0.25)\mathrm{log}\;P+(0.03\pm 0.07)[\mathrm{Fe}/{\rm{H}}]\\ & & -(1.25\pm 0.06).\end{eqnarray} \tag{ 6 }$$

The situation improves, but there is still a difference of ∼0.2 mag between the two zero-point obtained based on the distance to the LMC (Equation (4)) and the one based on the HST parallaxes of four RR Lyrae stars (Equation (6)). In fact, if we apply our ${\mathrm{PL}}_{{K}_{{\rm{s}}}}{\rm{Z}}$ relation with zero-point calibrated on Benedict et al. (2011) parallaxes (Equation (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\pm 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.'s (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\pm 0.065$ mag), published by Guinan et al. (1998), Fitzpatrick et al. (2002, 2003), and Ribas et al. (2002). Furthermore, in using the late-type EBs to calibrate the RR Lyrae ${\mathrm{PL}}_{{K}_{{\rm{s}}}}{\rm{Z}}$ relation we have implicitly assumed that RR Lyrae stars and EBs are at the 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 ${\mathrm{PL}}_{{K}_{{\rm{s}}}}{\rm{Z}}$ relation is the same in the MW and in the LMC, which may not be true (see Section 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 noting that, according to van Leeuwen (2007), the 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.'s (2011) 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 ${\mathrm{PL}}_{{K}_{{\rm{s}}}}{\rm{Z}}$ relation is expected from the ESA astrometric satellite Gaia. We discuss this topic in Section 4.

In spite of many studies in the literature, it still remains unsettled whether the RR Lyrae ${\mathrm{PL}}_{{K}_{{\rm{s}}}}{\rm{Z}}$ is a universal relation. To investigate this matter we have derived the ${\mathrm{PL}}_{{K}_{{\rm{s}}}}{\rm{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 absolute magnitudes known from B–W studies based on near-infrared data (Jones et al. 1988, 1992; Fernley et al. 1990; Liu & Janes 1990; Cacciari et al. 1992; Skillen et al. 1993; Fernley 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 field RR Lyrae stars analyzed 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 suspicion that the star is a dwarf Cepheid (see Feast et al. 2008 and discussion therein). Furthermore, following Fernley (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.

Table 2.  Properties of 23 Bright RR Lyrae Stars in the Milky Way

Star	R.A.	Decl.	P	[Fe/H]	σ[Fe/H]	$E(B-V)$	V	K	Ks	MV	MK	${\sigma }_{{M}_{V},{M}_{K}}$
	(deg)	(deg)	(day)	(dex)	(dex)	(mag)	(mag)	(mag)	(mag)	(mag)	(mag)	(mag)
UU Cet	1.02135	−16.99764	0.606081	−1.38	0.08	0.015	12.08	10.85	10.837	0.62	−0.55	0.15
SW And	5.92954	29.40101	0.442279	−0.06	0.08	0.045	9.71	8.54	8.505	0.94	−0.07	0.15
RR Cet	23.03405	1.34173	0.553025	−1.38	0.08	0.015	9.73	8.56	8.520	0.68	−0.42	0.15
X Ari	47.12869	10.44590	0.651139	−2.50	0.08	0.18	9.57	7.95	7.941	0.57	−0.575	0.15
AR Per	64.32165	47.40018	0.425549	−0.34	0.16	0.31	10.51	8.66	8.642	0.87	−0.06	0.25
RX Eri	72.43455	−15.74118	0.587246	−1.63	0.16	0.03	9.69	8.42	8.429	0.66	−0.51	0.25
RR Gem	110.38971	30.88318	0.397316	−0.32	0.16	0.075	11.38	10.26	10.275	0.89	−0.02	0.25
TT Lyn	135.78245	44.58559	0.597438	−1.64	0.16	0.015	9.86	8.65	8.611	0.65	−0.55	0.15
T Sex	148.36833	2.05732	0.324698	−1.20	0.16	0.015	10.04	9.18	9.200	0.66	−0.07	0.25
RR Leo	151.93108	23.99176	0.452387	−1.37	0.16	0.015	10.73	9.70	9.730	0.76	−0.16	0.15
WY Ant	154.02061	−29.72845	0.574312	−1.32	0.16	0.06	10.87	9.64	9.674	0.55	−0.55	0.15
W Crt	171.62351	−17.91435	0.412013	−0.89	0.16	0.03	11.54	10.56	10.539	0.96	+0.08	0.15
TU UMa	172.45205	30.06733	0.557659	−1.38	0.16	0.015	9.82	8.67	8.660	0.70	−0.41	0.15
UU Vir	182.14613	−0.45676	0.475606	−0.64	0.16	0.015	10.56	9.51	9.414	0.80	−0.195	0.15
SW Dra	184.44429	69.51062	0.569670	−0.91	0.16	0.015	10.48	9.33	9.319	0.68	−0.38	0.15
RV Oct	206.63230	−84.40177	0.571130	−1.92	0.16	0.09	10.98	9.51	9.526	0.68	−0.40	0.15
TV Boo	214.15242	42.35992	0.312559	−2.31	0.16	0.015	10.97	10.22	10.248	0.58	−0.20	0.25
RS Boo	218.38839	31.75462	0.377337	−0.37	0.16	0.015	10.37	9.45	9.507	0.85	+0.00	0.25
VY Ser	232.75803	1.68382	0.714094	−1.71	0.08	0.03	10.13	8.78	8.826	0.61	−0.665	0.25
V445 Oph	246.17171	−6.54165	0.397023	+0.17	0.08	0.195	11.05	9.24	9.262	1.09	+0.30	0.25
TW Her	268.63000	30.41048	0.399601	−0.58	0.16	0.06	11.28	10.22	10.239	0.80	−0.07	0.15
AV Peg	328.01164	22.57483	0.390380	−0.03	0.16	0.06	10.50	9.36	9.346	1.10	+0.14	0.15
RV Phe	352.13106	−47.45362	0.596413	−1.69	0.16	0.015	11.94	10.72	10.768	0.86	−0.29	0.25

Note. The columns report: (1) name of the star; (2) R.A. (J2000) from SIMBAD database; (3) decl. (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 (1994); (12) absolute magnitude in the K passband obtained from B–W analyses and corrected to p = 1.38 (see the text for details); (13) errors in the absolute V, K magnitudes.

Download table as:  ASCIITypeset image

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 ${\rm{\Delta }}S$ index (Preston 1959), from which metal abundances can be derived for the stars which lack abundance analysis. For the 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 Fe i and Fe ii abundances, adopting log (Fe i) = 7.56 and log (Fe ii) = 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 Fe ii abundance and adopting log (Fe i) = 7.51 for the Sun. In Appendix B we present the PLZ relation obtained using Lambert et al.'s (1996) metallicities and our approach.

By applying our fitting approach (Section 3.1) to the 23 MW RR Lyrae stars we derived the following ${\mathrm{PL}}_{K}{\rm{Z}}$ relation:

$$\begin{eqnarray}{M}_{K} & = & (-2.53\pm 0.36)\mathrm{log}\;P+(0.07\pm 0.04)[\mathrm{Fe}/{\rm{H}}]\\ & & -(0.95\pm 0.14).\end{eqnarray} \tag{ 7 }$$

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 because 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 log P 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, but 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 ${\mathrm{PL}}_{{K}_{{\rm{s}}}}{\rm{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 ${\mathrm{PL}}_{K}{\rm{Z}}$ relation (Equation (7)) on the $\mathrm{log}\;(P)-{M}_{K}$ and ${M}_{K}-[\mathrm{Fe}/{\rm{H}}]$ planes are shown in Figure 3. Gray lines are the same as in Figure 2 and are described in Section 3.1.

Zoom In Zoom Out Reset image size

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}_{{\rm{s}}}$ photometry in the VISTA system.⁹ To address this issue we have reported in column 10 of Table 2 the average ${K}_{{\rm{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}_{{\rm{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 significantly affect our comparison.

The near-infrared ${\mathrm{PL}}_{K}{\rm{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.'s (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 ${\mathrm{PL}}_{K}{\rm{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 ${\mathrm{PL}}_{K}{\rm{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:

$$\begin{eqnarray}&&{M}_{K}=-2.353\mathrm{log}P+0.175\mathrm{log}Z-0.597.\end{eqnarray} \tag{ 8 }$$

By using Equations (9) and (10) in Catelan et al. (2004) and assuming [α/Fe] ∼ 0.3 (e.g., Carney 1996) we transformed Equation (8) into the form, presented in the second row of Table 3.

Table 3.  ${\mathrm{PL}}_{{K}_{{\rm{s}}}}{\rm{Z}}$ Relations from the Literature

Relation	a	b	ZPa from the Original Relation	ZPb from Benedict et al. (2011)
Theoretical or semi-theoretical relations
Bono et al. (2003)	−2.101	0.231 ± 0.012	−0.770 ± 0.044	−0.58 ± 0.04
Catelan et al. (2004)	−2.353	0.175	−0.869	...
Empirical relations
Dall'Ora et al. (2004)	−2.16 ± 0.09	...	...	−0.56 ± 0.02
Del Principe (2006) c,d	−2.71 ± 0.12	0.12 ± 0.04	...	−0.57 ± 0.02
Sollima et al. (2006) e	−2.38 ± 0.04	0.08 ± 0.11	−1.05 ± 0.13	−0.57 ± 0.03f
Sollima et al. (2008) e	−2.38 ± 0.04	0.08 ± 0.11	−1.07 ± 0.11	−0.56 ± 0.02
Borissova et al. (2009) c,g	−2.11 ± 0.17	0.05 ± 0.07	−1.05	−0.56 ± 0.03
This paper (LMC)c,g	−2.73 ± 0.25	0.03 ± 0.07	−1.06 ± 0.01	−0.55 ± 0.06i
This paper (MW)h	−2.53 ± 0.36	0.07 ± 0.04	−0.95 ± 0.14	−0.56 ± 0.06i

Notes.

^aZero-point of the original relation from the literature in the form: $a\mathrm{log}\;P+b[\mathrm{Fe}/{\rm{H}}]+\mathrm{ZP}$. ^bZero-point of the relation in the form: $a(\mathrm{log}\;P+0.28)+b([\mathrm{Fe}/{\rm{H}}]+1.58)+\mathrm{ZP}$, as recalibrated by Benedict et al. (2011). ^cNear-infrared photometry in the K_s band. ^dMetallicity is on the Zinn & West (1984) metallicity scale. ^eMetallicity is on the Carretta & Gratton (1997) metallicity scale. ^fZero-point was calibrated by Benedict et al. (2011) neglecting the metallicity term. ^gMetallicity on the scale adopted by Gratton et al. (2004), which is , on average, 0.06 dex higher than Zinn & West (1984) scale. ^hMetallicities from Clementini et al. (1995), they are on the scale of the high dispersion spectra (i.e., the Carretta et al. 2009). ⁱZero-point was calibrated by us considering only four RR Lyrae stars (XZ Cyg, UV Oct, SU Dra, and RR Lyr) from Benedict et al. (2011) and excluding RZ Cep, since there are concerns about the parallax of this star.

Download table as:  ASCIITypeset image

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 (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 ${\mathrm{PL}}_{K}{\rm{Z}}$ relation by analyzing 538 RR Lyrae stars in 15 Galactic clusters and in the LMC globular cluster Reticulum. This relation spans the metallicity range $-2.15\lt [\mathrm{Fe}/{\rm{H}}]\lt -0.9$ dex. Mean K magnitudes were estimated by combining 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.'s (2006) ${\mathrm{PL}}_{K}{\rm{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 five 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.'s (2008) data set, and derived the ${\mathrm{PL}}_{{K}_{{\rm{s}}}}{\rm{Z}}$ relation based on the total sample of 107 LMC variables. The zero-point was calculated using Sollima et al.'s (2008) mean K magnitude, the reddening, and the trigonometric parallax of RR Lyrae.

Benedict et al. (2011) recalibrated all the literature relations listed in Table 3, but Catelan et al. (2004)'s one, by fitting to equations in the form: $a(\mathrm{log}\;P+0.28)+b([\mathrm{Fe}/{\rm{H}}]+1.58)+\mathrm{ZP}$, the Lutz–Kelker–Hanson-corrected absolute magnitudes of the five MW RR Lyrae stars for which HST parallaxes are available. Since there are concerns about Benedict et al.'s (2011) RR Lyrae star (namely RZ Cep) we transformed the ${\mathrm{PL}}_{K}{\rm{Z}}$ relations of the LMC and MW RR Lyrae stars derived in this paper (Equations (4), (6)) to the form adopted by Benedict et al. (2011) and determined their zero-points on the basis of Benedict et al.'s (2011) HST parallaxes but excluding RZ Cep. The zero-points based on Benedict et al. (2011) parallaxes are presented in column 5 of Table 3, whereas the zero-point of our LMC ${\mathrm{PL}}_{{K}_{s}}{\rm{Z}}$ relation calculated by assuming the distance to the LMC (Section 3.2) and the zero-point of our MW ${\mathrm{PL}}_{K}{\rm{Z}}$ relation based on the B–W studies are (Section 3.3) are presented in Column 4.

The slope in the period of the RR Lyrae ${\mathrm{PL}}_{K}{\rm{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 (2006), whereas the slope of the MW RR Lyrae ${\mathrm{PL}}_{K}{\rm{Z}}$ is in good agreement with that derived by Sollima et al. (2006, 2008).

The dependence on metallicity of the ${\mathrm{PL}}_{K}{\rm{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 (2006) are in the Zinn & West (1984) scale, whereas in Sollima et al. (2006, 2008) they 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 the Zinn & West scale. This difference is small and systematic, hence it 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 the Zinn & West scale using the relations provided in Carretta et al. (2009).

The slope in metallicity of the ${\mathrm{PL}}_{K}{\rm{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.'s (2009) 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 with on the RR Lyrae ${\mathrm{PL}}_{K}{\rm{Z}}$, independently of the specific environment.

Using parallax data for the 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). Nonlinear 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\mathrm{log}(\varpi )$, where ϖ is the parallax. Additionally, stars with a negative parallax measurement cannot 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 modeling the data and inferring the slope and zero-point of the relation via statistical methods. For a catalog of N stars we can define ${\boldsymbol{x}}=({{\boldsymbol{x}}}_{1},{{\boldsymbol{x}}}_{2},...,{{\boldsymbol{x}}}_{N})$ where the vector $({{\boldsymbol{x}}}_{i}=m,l,b,P,\varpi ,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 $({{\boldsymbol{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}=\rho {\rm{log}}{P}_{0}+\delta$, 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:

$$\begin{eqnarray}&&{\varphi }_{M}({{\boldsymbol{x}}}_{0}| \rho ,\delta ,{\sigma }_{\mathrm{PL}})=\displaystyle \frac{1}{{\sigma }_{\mathrm{PL}}\sqrt{2\pi }}{e}^{-0.5{\left(\frac{{M}_{0}-(\rho {\rm{log}}{P}_{0}+\delta )}{{\sigma }_{\mathrm{PL}}}\right)}^{2}}\end{eqnarray} \tag{ 9 }$$

where ${\sigma }_{\mathrm{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

$$\begin{eqnarray}&&{M}_{0}={m}_{0}+5\mathrm{log}({\varpi }_{0})+5-{A}_{0}.\end{eqnarray} \tag{ 10 }$$

The observations are normally distributed around the true values with a standard deviation given by the formal error on the measurement

$$\begin{eqnarray}{\mathcal{E}}({\boldsymbol{x}}| {{\boldsymbol{x}}}_{0}) & = & \displaystyle \frac{1}{{\epsilon }_{\varpi }{\epsilon }_{m}{\epsilon }_{A}{(2\pi )}^{3/2}}{e}^{-0.5{\left(\frac{\varpi -{\varpi }_{0}}{{\epsilon }_{\varpi }}\right)}^{2}}{e}^{-0.5{\left(\frac{m-{m}_{0}}{{\epsilon }_{m}}\right)}^{2}}\\ & & \times {e}^{-0.5{\left(\frac{A-{A}_{0}}{{\epsilon }_{A}}\right)}^{2}}\delta ({l}_{0},{b}_{0},{P}_{0}).\end{eqnarray} \tag{ 11 }$$

Assuming negligible errors on the position and period, the observations are described by a delta function. The terms ${\epsilon }_{\varpi }$, ${\epsilon }_{m}$, ${\epsilon }_{A}$ are the formal errors on the parallax, magnitude, and extinction.

With the above models defined, the joint PDF for the observations is

$$\begin{eqnarray}&&{\mathcal{P}}({{\boldsymbol{x}}}_{i}| \rho ,\delta ,{\sigma }_{\mathrm{PL}})\\ &&\quad ={\mathcal{S}}({\boldsymbol{x}}){\displaystyle \int }_{\forall {{\boldsymbol{x}}}_{0}}{\varphi }_{M}({{\boldsymbol{x}}}_{0}| \rho ,\delta ,{\sigma }_{\mathrm{PL}}){\mathcal{E}}({\boldsymbol{x}}{| {\boldsymbol{x}}}_{0})d{{\boldsymbol{x}}}_{{\bf{0}}};\end{eqnarray} \tag{ 12 }$$

the "true" parameters ${{\boldsymbol{x}}}_{0}$ are never known and so these values are marginalized through integration. The term ${\mathcal{S}}({\boldsymbol{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

$$\begin{eqnarray}{\mathcal{S}}({\boldsymbol{x}})=\left\{\begin{array}{ll}1, & \mathrm{if}\;\ G\lt 20,\\ 0, & {\rm{otherwise}}.\end{array}\right.\end{eqnarray} \tag{ 13 }$$

The Maximum Likelihood Estimation of the parameters is found by maximizing Equation (12) by varying the parameters (α, ρ, ${\sigma }_{\mathrm{PL}}$). This formulation avoids nonlinear transformations on error-affected data and includes a selection function which avoids the Malmquist bias.

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 catalog, 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 catalog. 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 color, 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

$$\begin{eqnarray}&&{M}_{K}=-2.53\;{\rm{log}}P-0.95,\end{eqnarray} \tag{ 14 }$$

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. Color information as $(V-I)$ is obtained from the Hipparcos catalog (Perryman & ESA 1997) where available. The apparent magnitude, position, color, and period data form the basis of a synthetic catalog of RR Lyrae stars, along with the distance obtained from the ${\mathrm{PM}}_{K}{\rm{Z}}$ relation, and is used as the input catalog of "true" parameters for GOG.

GOG then creates simulated Gaia observations for our sample. We take the ${\mathrm{PM}}_{K}{\rm{Z}}$ elation (Equation (14)) as true, as a study of the possible precision in ${\mathrm{PM}}_{K}{\rm{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 ${\mathrm{PM}}_{K}$ relation of

$$\begin{eqnarray}&&{M}_{K}=(-2.531\pm 0.038)\mathrm{log}P+(-0.95\pm 0.01)\end{eqnarray} \tag{ 15 }$$

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 become available for the much larger sample of RR Lyrae stars, we will apply the described method to fit the ${\mathrm{PL}}_{K}{\rm{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 the zero-point of the ${\mathrm{PL}}_{{K}_{{\rm{s}}}}{\rm{Z}}$ relation based on the LMC RR Lyrae stars. This provides a hint of what will be possible to achieve with Gaia parallaxes.

Because this case contains errors in more than one axis, they can be put together into a covariance tensor ${{\boldsymbol{S}}}_{i}$

$$\begin{eqnarray}{{\boldsymbol{S}}}_{i}\equiv \left[\begin{array}{ccc}{\sigma }_{{xi}}^{2} & {\sigma }_{{xyi}} & {\sigma }_{{xzi}}\\ {\sigma }_{{xyi}} & {\sigma }_{{yi}}^{2} & {\sigma }_{{yzi}}\\ {\sigma }_{{xzi}} & {\sigma }_{{yzi}} & {\sigma }_{{zi}}^{2}\end{array}\right].\end{eqnarray} \tag{ 23 }$$

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

$$\begin{eqnarray}&&p({x}_{i},{y}_{i},{z}_{i}{| {\boldsymbol{S}}}_{i},x,y,z)=\displaystyle \frac{1}{2\;\pi \;\sqrt{\mathrm{det}({{\boldsymbol{S}}}_{i})}}\\ &&\quad \times \mathrm{exp}\left(-\displaystyle \frac{1}{2}\;{\left[{{\boldsymbol{Z}}}_{i}-{\boldsymbol{Z}}\right]}^{\top }\;{{\boldsymbol{S}}}_{i}^{-1}\left[{{\boldsymbol{Z}}}_{i}-{\boldsymbol{Z}}\right]\right),\end{eqnarray} \tag{ 24 }$$

where we have implicitly made column vectors

$$\begin{eqnarray}&&{\boldsymbol{Z}}=\left[\begin{array}{c}x\\ y\\ z\end{array}\right];\quad {{\boldsymbol{Z}}}_{i}=\left[\begin{array}{c}{x}_{i}\\ {y}_{i}\\ {z}_{i}\end{array}\right].\end{eqnarray} \tag{ 25 }$$

In only two dimensions (e.g., y and z), the slope (e.g., B) can be described by a unit vector $\hat{{\boldsymbol{v}}}$ orthogonal to the line or linear relation (at any x)

$$\begin{eqnarray}&&\hat{{\boldsymbol{v}}}=\displaystyle \frac{1}{\sqrt{1+{B}^{2}}}\left[\begin{array}{c}-B\\ 1\end{array}\right]=\left[\begin{array}{c}-\mathrm{sin}\theta \\ \mathrm{cos}\theta \end{array}\right],\end{eqnarray} \tag{ 26 }$$

where the angle $\theta =\mathrm{arctan}B$ is made between the line and the y axis. The orthogonal displacement ${{\rm{\Delta }}}_{i}$ of each data point $({y}_{i},{z}_{i})$ from the line is given by

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{i}={\hat{{\boldsymbol{v}}}}^{\top }\left[\begin{array}{c}{y}_{i}\\ {z}_{i}\end{array}\right]-C\;\mathrm{cos}\theta .\end{eqnarray} \tag{ 27 }$$

Instead of extending fully into three dimensions, we will assume a negligible error in x (which will be the period, so that it has justifiably higher precision). The value of x can therefore be input directly into ${{\rm{\Delta }}}_{i}$ without worrying about the interplay between the other parameters:

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{i}={\hat{{\boldsymbol{v}}}}^{\top }\left[\begin{array}{c}{y}_{i}\\ {z}_{i}\end{array}\right]-(C\;\mathrm{cos}\theta +A\;{x}_{i}).\end{eqnarray} \tag{ 28 }$$

Assuming negligable errors in x also redefines the covariance matrix of the errors as

$$\begin{eqnarray}{{\boldsymbol{S}}}_{i}\equiv \left[\begin{array}{cc}{\sigma }_{{yi}}^{2} & {\sigma }_{{yzi}}\\ {\sigma }_{{yzi}} & {\sigma }_{{zi}}^{2}\end{array}\right].\end{eqnarray} \tag{ 29 }$$

Similarly, each data point's covariance matrix ${\boldsymbol{S}}$_i projects down to an orthogonal variance ${{\rm{\Sigma }}}_{i}^{2}$ given by

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{i}^{2}={\hat{{\boldsymbol{v}}}}^{\top }\;{{\boldsymbol{S}}}_{i}\;\hat{{\boldsymbol{v}}}\end{eqnarray} \tag{ 30 }$$

and then the log likelihood for $(A,B,C)$ or $(A,\theta ,C\;\mathrm{cos}\theta )$ can be written as

$$\begin{eqnarray}&&\mathrm{ln}{\mathcal{L}}=K-\underset{i=1}{\overset{N}{\displaystyle \sum }}\displaystyle \frac{{{\rm{\Delta }}}_{i}^{2}}{2\;{{\rm{\Sigma }}}_{i}^{2}},\end{eqnarray} \tag{ 31 }$$

where K is some constant. This likelihood can be maximized to find A, B, and C.

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 ${{\rm{\Sigma }}}_{i}^{2}$ defined above, with the Gaussian intrinsic scatter of variance V. Therefore the likelihood becomes

$$\begin{eqnarray}&&\mathrm{ln}{\mathcal{L}}=K-\underset{i=1}{\overset{N}{\displaystyle \sum }}\displaystyle \frac{1}{2}\;\mathrm{ln}({{\rm{\Sigma }}}_{i}^{2}+V)-\underset{i=1}{\overset{N}{\displaystyle \sum }}\displaystyle \frac{{{\rm{\Delta }}}_{i}^{2}}{2\;[{{\rm{\Sigma }}}_{i}^{2}+V]},\end{eqnarray} \tag{ 32 }$$

where again K is a constant and everything else is defined as above. We then solve for A, B, C, and V by maximizing the log likelihood. The optimization is performed using the adaptive MCMC sampler EMCEE developed by Foreman-Mackey et al. (2013). Any optimization 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.