A Probabilistic Approach to Fitting Period-Luminosity Relations and Validating Gaia Parallaxes

Pulsating stars, such as Cepheids, Miras, and RR Lyrae stars, are important distance indicators and calibrators of the"cosmic distance ladder", and yet their period-luminosity-metallicity (PLZ) relations are still constrained using simple statistical methods that cannot take full advantage of available data. To enable optimal usage of data provided by the Gaia mission, we present a probabilistic approach that simultaneously constrains parameters of PLZ relations and uncertainties in Gaia parallax measurements. We demonstrate this approach by constraining PLZ relations of type $ab$ RR Lyrae stars in near-infrared W1 and W2 bands, using Tycho-Gaia Astrometric Solution (TGAS) parallax measurements for a sample of $\approx100$ type $ab$ RR Lyrae stars located within 2.5 kpc of the Sun. The fitted PLZ relations are consistent with previous studies, and in combination with other data, deliver distances precise to 6% (once various sources of uncertainty are taken into account). To a precision of 0.05 mas ($1\sigma$), we do not find a statistically significant offset in TGAS parallaxes for this sample of distant RR Lyrae stars (median parallax of 0.8 mas and distance of 1.4 kpc). With only minor modifications, our probabilistic approach can be used to constrain PLZ relations of other pulsating stars, and we intend to apply it to Cepheid and Mira stars in the near future.


INTRODUCTION
One of the main goals of observational astronomy is to ever more precisely and accurately measure the distance to astrophysical objects. Measuring distances to individual stars is critical to understanding a wide range of astronomical phenomena, from stellar structure to Galactic dynamics. An important tool in this endeavor are periodically-pulsating stars, such as Cepheids, Miras, and RR Lyrae stars, whose absolute magnitudes can be predicted using a period-luminosity (PL) relation 1 , and whose period-luminosity relations can be calibrated using trigonometric parallax measurements (e.g., Feast & Catchpole 1997;van Leeuwen et al. 1997;Benedict et al. 2007Benedict et al. , 2011. These stars and their PL relations are crucial, as they tie the extragalactic distance scale to the local one, and make the important second rung in the "cosmic distance ladder". The trigonometric parallaxes ( ) obtained by the Tycho-Gaia Astrometric Solution (TGAS; Michalik et al. 2015), and made public through the Gaia Data Release 1 (Gaia Collaboration et al. 2016;Lindegren et al. 2016), present us with an exciting opportunity to recalibrate PL relations and potentially improve the accuracy of the cosmic distance ladder. However, the majority of Cepheids, Miras, and RR Lyrae bsesar@mpia.de 1 The period-luminosity relation for Cepheids is also known as the Leavitt law (Leavitt 1908;Leavitt & Pickering 1912). stars in the TGAS sample are distant, and consequently, their parallaxes have a high fractional uncertainty (σ / > 0.2). When a parallax has a high fractional uncertainty, the posterior probability distribution over its true distance is complex and non-Gaussian, even when the uncertainties in parallax are Gaussian (Bailer-Jones 2015; Astraatmadja & Bailer-Jones 2016). This poses a serious problem for the traditional approach to PL relation fitting, which assumes Gaussian uncertainties in distance, and does a weighted least-squares fit in the distance (i.e., absolute magnitude) vs. period plane. We can always, of course, fit PL relations by using only stars with precise parallaxes (and thus precise distances with close to Gaussian errors), but then we would ignore a large amount of potentially useful data (that was also obtained at a significant cost).
To avoid ignoring valuable data, in this paper we present a probabilistic approach to inferring period-luminosity(metallicity) relations that makes a full use of parallax and other measurements, irrespective of their precision. To demonstrate the approach in practice, we use it to constrain period-luminosity-metallicity (PLZ) relations for type ab RR Lyrae stars (i.e., fundamental mode pulsators ;Smith 2004) in the near-infrared W1 and W2 bands used by the WISE mission (Wright et al. 2010). While the demonstration is done using type ab RR Lyrae stars (hereafter, RRab), we emphasize that the approach can be easily applied to other pulsating stars, by simply plugging different data and prior information into our framework.
In addition to constraining PLZ relations, our approach en-ables a straightforward validation of parallax measurements and their uncertainties. The validation of parallaxes is important as several studies have reported an offset in TGAS parallaxes (Stassun & Torres 2016;Jao et al. 2016) and an overestimation of TGAS parallax uncertainties (Casertano et al. 2016;Gould et al. 2016). Our paper is organized as follows. In Section 2, we describe the astrometric, photometric, and spectroscopic data used to calibrate PLZ relations for RRab stars in WISE W1 and W2 bands, and validate TGAS parallaxes. The construction of the likelihood function that is at the core of our probabilistic approach, is presented in Section 3. In Section 4, we present the results of applying the probabilistic method of Section 3, to data described in Section 2. We summarize and discuss our results in Section 5.
2. DATA To constrain the PLZ relations for RR Lyrae stars in WISE W1 and W2 bands, we use TGAS trigonometric parallaxes ( ), spectroscopic metallicities ([Fe/H]; Fernley et al. 1998), log-periods (log P , base 10), and apparent magnitudes 2 (m; Klein et al. 2014) for 102 RRab stars within ≈ 2.5 kpc from the Sun. The E(B − V ) reddening at a star's position is obtained from the Schlegel et al. (1998) dust map. We denote this data set as D = {d k }, where d k = { , [Fe/H], log P, m, EBV } is the data set associated with the kth star.
We note that the extinction correction changes by less than 0.01 mag for 96% of objects (0.01 to 0.02 mag for the remaining 4%), if instead of the E(B − V ) reddening from the Schlegel et al. (1998) dust map (which is essentially the upper limit on reddening at high galactic latitudes), we adopt effective reddening E(B − V )(1 − exp(−h/h dust )), where h is the distance from the Galactic plane, and h dust = 130 pc is the scale height of the dust layer (Gould & Popowski 1998;Klein et al. 2011). Such a small change is expected given the small values of extinction coefficients in near-IR W1 and W2 bands (0.13 and 0.17), the reddening distribution (0.18 mag in the 95th percentile), and the spatial distribution of RR Lyrae stars under consideration (most of them are located beyond a few scale heights of the galactic plane).
Briefly, the stars in our sample have metallicities ranging from −0.1 dex to −2.6 dex, and periods ranging from 0.36 to 0.73 days. The uncertainty in [Fe/H] is the same for all stars, σ [Fe/H] = 0.15 dex (see Note 8 in Table 1 of Fernley et al. 1998). Based on the analysis of Klein et al. (2011, see their Appendix A), we assume that the uncertainty in logperiod is σ log P = 0.02 log P . The uncertainty in E(B − V ) is assumed to be σ EBV = 0.1E(B − V ) (Schlegel et al. 1998). The average uncertainty in apparent magnitudes is σ m = 0.005 mag .
The TGAS parallax measurements have reported uncertainties ranging from 0.22 mas to 0.47 mas (5th to 95th percentile), with a median of 0.28 mas. The median fractional uncertainty of our sample is σ / = 0.17. As described by 2 Averaged in flux over a pulsation period. Lindegren et al. (2016), the reported parallax uncertainties were calculated as (their Equation 4) where ς is the formal parallax uncertainty, f = 1.4, and σ ,add = 0.2 mas. To facilitate validation of TGAS parallax uncertainties in Section 3, we calculate ς values as The Klein et al. (2014) sample also contains 16 type c RR Lyrae stars (i.e., first overtone pulsators; RRc) with measured W1 and W2 apparent magnitudes. However, we do not use these stars to constrain the PLZ relations of RRc stars because the sample is quite small, and because the median fractional parallax uncertainty of the sample is quite large (σ / ≈ 0.5). Instead, we use these stars to verify whether RRab and RRc stars follow the same PLZ relations. Before making such comparisons, the practice has been to first "fundamentalize" the periods of RRc stars (e.g., Dall'Ora et al. 2004). We follow the same practice, and calculate fundamentalized periods of RRc stars as where P f o is the original period. The data used in this work are provided in a machinereadable format in the electronic version of the Journal (filename "data.csv"). The is a standard normal random variable with zero mean and variance (aσ logP ) 2 + bσ [Fe/H] 2 + σ 2 M , where σ M accounts for the scatter in absolute magnitude M due to modeling uncertainties. When interpreting this scatter, however, it is important to keep in mind that, in reality, σ M also includes unaccounted measurement uncertainties (added in quadrature). Thus, σ M represents the so-called "intrinsic" scatter in a PLZ relation only if the measurement uncertainties are correctly estimated (which is difficult to do in practice).
As mentioned above, we also allow freedom in the error model used for the parallax measurements. We model the TGAS parallax measurements as being drawn from a Gaussian distribution centered on where r is the true heliocentric distance (in parsecs), and 0 represents the global offset of TGAS parallaxes with respect to the inverse distances (e.g., an offset due to the impact of imperfectly modeled basic-angle variations on the astrometric solution; Lindegren et al. 2016). The standard deviation of this Gaussian is equal to the uncertainty in the TGAS parallax, which we model using Equation 1, where f and σ ,add are also included as free parameters. The expressions defined by Equations 5 and 1 were motivated by conclusions of recent studies that have found TGAS parallaxes to be offset (Jao et al. 2016;Stassun & Torres 2016;De Ridder et al. 2016), and studies that have found TGAS parallax uncertainties to be overestimated (Gould et al. 2016;Casertano et al. 2016). 3 To constrain the PLZ relation in a probabilistic manner we need to calculate the joint posterior probability add }, the scale length parameter L (used in the distance prior, see Equation 18 below), and the set of nuisance parameters α k = {r, log P int , [Fe/H] int , EBV int } k that represent the true distance r, intrinsic log-period log P int , metallicity [Fe/H] int , and reddening EBV int for each star. For conciseness, we also define θ = (θ PLZ , θ , L). For this work, our main interest is in the marginal posterior probability of the TGAS validation and PLZ parameters p(θ | D), which is related to the marginal likelihood p(D | θ) through where p(θ) is the prior probability of the parameter value set θ, and is the marginal likelihood of the full data set given the combined PLZ, TGAS validation, and L parameters (θ), and assuming independent data points. The (un-marginalized) likelihood function for the kth star is given as follows: where: Fernley+98 Figure 1. The probabilistic graphical model that describes the dependencies between model parameters and data used in this work.
Double circles indicate likelihoods, single orange circles indicate nuisance parameters, while single green circles indicate model parameter sets. Fixed parameters, such as the extinction coefficient (A λ ; Schlafly et al. 2016), priors on nuisance parameters (e.g., log P int range), as well as data sources (e.g., Gaia), are not enclosed in circles. Parameters inside the square are specific to the kth star, while those on the outside are global. The arrows indicate conditional dependence. For example, the arrows from r and θ to indicate that the observed parallax depends on the heliocentric distance r and TGAS parallax validation parameters θ (i.e., p( | r, θ )).
is a normal distribution centered on µ with variance σ 2 , and Our statistical model is also visualized as a probabilistic graphical model (PGM) shown in Figure 1. One of the advantages of PGMs is that they enable a straightforward examination of dependencies between data and model parameters (Koller & Friedman 2009). For example, while the observed apparent magnitude m depends on the adopted extinction coefficient A λ , and parameters in α k and θ PLZ sets, the observed parallax depends only on distance and TGAS parallax validation parameters, θ , described above.
In Equations 12 to 14, we model the observed log-period, metallicity, and reddening as being drawn from Gaussian dis-tributions centered on some "intrinsic" values (represented by the superscript "int"), and with a standard deviation equal to the uncertainty in measurement or the model. The "int" parameters are "intrinsic" in the sense that they are not affected by uncertainties (model or observational).
Note that Equation 9 contains 7 global parameters (4 in θ PLZ and 3 in θ parameter sets), and 4 nuisance parameters for each star, α k . Since we are (at the moment) not interested in these nuisance parameters, we marginalize (i.e. integrate) Equation 9 over these parameters.
We use the following prior probability distributions for the nuisance parameters α k . For [Fe/H] int and log P int we choose priors that are uniform in ranges appropriate for RR Lyrae stars: 0 dex to -3 dex, and -1.0 log(day) to 0 log(day). For EBV int , we adopt a uniform prior in the 0 mag to 1 mag range (appropriate for low-extinction regions, as is the case here), and for distance we adopt an exponentially decreasing volume density prior with a scale length parameter L (Bailer-Jones 2015; Astraatmadja & Bailer-Jones 2016) The prior p(r | L) is positive in the 200 < r/pc < 2700 distance range, and zero elsewhere. The marginal likelihood for the kth star is which can be written in a simpler form by analytically performing the integrals over the Gaussian expressions in log P int , [Fe/H] int , and EBV int : where The likelihood for the entire data set D can now be calculated using Equation 7. Before we can calculate the (marginal) posterior distribution (Equation 6), we need to define prior probabilities for the global parameters θ PLZ , θ , and L. For σ M , σ ,add , and L parameters, we adopt Jeffreys log-uniform priors (p(x) ∝ 1/x; Jaynes 1968), and for the b parameter that scales the absolute magnitude with [Fe/H], we adopt a uniform prior that is positive for 0 < b/mag dex −1 < 3 (based on stellar evolution and pulsation models of RR Lyrae stars; Marconi et al. 2015). Since the parameter must be positive in the 200 pc to 2700 pc range (Equation 5), for 0 we adopt a uniform prior that is positive in the −5 < 0 /mas < 20 range. For the remaining model parameters we adopt wide uniform priors.
To efficiently explore the parameter space, we use the Goodman & Weare (2010)  We use 160 walkers and obtain convergence 5 after a burn-in phase of 1000 steps per walker. The chains are then evolved for another 1500 steps, and the first 1000 (burn-in) steps are discarded.
To describe the marginal posterior distributions of individual model parameters, we measure the median, the difference between the 84th percentile and the median, and the difference between the median and the 16th percentile of each marginal posterior distribution (for a Gaussian distribution, these differences are equal to ±1 standard deviation). We report these and maximum a posteriori (MAP) values in Table 1.
NOTE-In each band, the first line lists the median, the difference between the 84th percentile and the median, and the difference between the median and the 16th percentile of each marginal posterior distribution. The second line provides the maximum a posteriori (MAP) values.
4. RESULTS In this Section, we discuss PLZ relations and TGAS parallax validation parameters that we have constrained by applying the method described in Section 3, to data described in Section 2. To illustrate the correlations between various model parameters, in Figure 2 we show two-dimensional posterior distributions for parameter pairs. The parameters we constrain also allow us to calculate more precise distances to RR Lyrae stars in our sample (Table 2).

PLZ Parameters
We find that, within the uncertainties, the PLZ parameters obtained using W1 band data are consistent with those obtained using W2 band data. This result is not too surprising given i) fairly large uncertainties in parameters, ii) the proximity of W1 and W2 bands in wavelength, and iii) the fact that the parameters of PL(Z) relations change little with wavelength for bands redder than the H or K band (Catelan et al. 2004;Marconi et al. 2015; Figure 4 of Madore et al. 2013). Since the two PLZ relations are consistent, we use PLZ parameters for the W2 band when comparing relations to previous studies.
Overall, our PLZ relations are consistent with PL(Z) relations found by previous studies (within 1σ of uncertainties). For example, Madore et al. (2013), Klein et al. (2014), and Dambis et al. (2014) measure the period slope of the W2 band PLZ relation for RRab stars to be a = −2.6 ± 0.9, a = −2.4 ± 0.2, and −2.3 ± 0.1, respectively, while we find a = −2.4 ± 0.8. Our estimate of the absolute magnitude at P = 1 day of −1.2 ± 0.2 mag, is consistent with ≈ −1.1±0.1 mag estimated by the above studies. Regarding the metallicity slope of the PLZ relation, our measurement of b = 0.17 ± 0.10 is consistent at the 1σ level with the Dambis et al. (2014) value of 0.12±0.02, and (in)consistent at the 2σ level with the slope of zero reported by Madore et al. (2013) and Klein et al. (2014) (note: the latter two studies did not report the uncertainty of their estimate).
The σ M MAP values 6 of 0.04 mag and 0.05 mag for W1 and W2 bands, respectively, seem a bit high given the expectation of a small intrinsic scatter in near-IR PLZ relations (Table 6 of Marconi et al. 2015). However, recall that σ M includes (via Equation 21) unaccounted measurement uncertainties in apparent magnitude, metallicity, log-period, and reddening. One possible source of unaccounted uncertainties may be the scatter in the WISE magnitude zeropoint, which is ≈ 0.03 mag (  Fernley et al. (1998) are actually a compilation of metallicities measured by different studies.

TGAS Parallax Validation Parameters
We measure the 0 parameter (that models the global offset of TGAS parallax measurements) to be consistent with zero within 0.05 mas (i.e., 0 = 0.00 ± 0.05 mas), indicating that there is no statistically significant offset in parallax, at least for this sample of distant RR Lyrae stars (median distance ≈ 1 kpc).
Five RR Lyrae stars in our sample (SU Dra, RR Lyr, UV Oct, XZ Cyg, and RZ Cep) also have parallaxes measured by Benedict et al. (2011) using Hubble Space Telescope (HST) astrometric observations. We find the average of the difference between TGAS and HST parallax measurements to be 0.04 mas, in agreement with our 0 = 0.00 ± 0.05 mas measurement.
Regarding the renormalization of TGAS formal parallax uncertainties ς (see Equation 1), we find σ ,add = 0.18 mas and that f is consistent with 1 within uncertainties(i.e., f = 0.7 ± 0.3). In principle, it is possible for f to be less than 1, if the formal TGAS parallax uncertainties, ς , are overestimated. Since ς are calculated from the inverse 5 × 5 normal matrix of the astrometric parameters (Lindegren et al. 2016), this could happen if there is some serious problem with the software that performs this calculation. As we do not believe this to be likely, we assume that f ≥ 1, or given data, f = 1.

Do RRc Stars Follow the PLZ Relation for RRab Stars?
Since we did not use RRc stars to constrain the model described in Section 3 (only RRab stars were used), we can now use this model to verify whether RRc stars follow the same PLZ relation as RRab stars.
Given our data and the model, the most appropriate comparison between the data and the model is the one between the observed TGAS parallaxes, , and inverse distances, 1/r. To infer the inverse distance, we first need to calculate the marginal posterior distribution of the heliocentric distance for the kth star, p(r k | d k , θ), since the marginal posterior distribution of the inverse distance, p(1/r k | d k , θ), is equal to Using relations of conditional probability, the marginal posterior distribution of the heliocentric distance for the kth star can be written as where Equation 23 is being integrated (i.e., marginalized) over the scale length parameter L, as well as parameters con-tained in θ PLZ and θ sets 8 . We evaluate Equation 22 over the 1/(2700 pc) < (1/r)/arcsec < 1/(200 pc) inverse distance range and average it over 160 samples of L, θ PLZ , and θ parameters taken from the last step of the Markov chain (which has 160 walkers). The mean inverse distance and its uncertainty are obtained by fitting a Gaussian to p(1/r k | d k , θ) evaluated on the above grid. Similarly, we evaluate Equation 23 over the 200 < r/pc < 2700 distance range, and list the distances and their uncertainties in Table 2. The median fractional uncertainty in distance is 6%. As the top panel in Figure 3 shows, the inferred inverse distances of RRc stars seem to systematically deviate from observed TGAS parallaxes as the log P decreases (i.e., as the periods get shorter). This trend suggests that the PLZ relation of RRc stars may have a shallower slope than the PLZ relation of RRab stars. Indeed, when both RRab and RRc stars are used to constrain the PLZ relations in the W1 and W2 bands (periods of RRc stars are fundamentalized), the resulting period slopes are shallower, though still consistent (within 1σ) with the period slopes obtained using only RRab stars (a RRab+RRc = −1.5 ± 0.6 mag dex −1 vs. a RRab = −2.4 ± 0.8 mag dex −1 ).
The above result (i.e., the flattening of the period slope when RRc stars are added) is consistent with the finding of Klein et al. (2014) that PLZ relations of RRc stars have shallower period slopes in W1 and W2 bands than RRab stars (though, the RRab and RRc period slopes they measure are consistent with each other within 1σ). On the other hand, the above result is at odds with Braga et al. (2015), who find that the period slopes of RRab and RRc stars in the K−band 9 are quite consistent, a ≈ −2.4 ± 0.2 mag dex −1 .
Relative to the observed TGAS parallaxes, the inverse distances of RRc stars also show an interesting V-shaped trend as a function of [Fe/H] (bottom panel of Figure 3). If real, this trend would suggest that, for RRc stars, the metallicity slope of the PLZ relation changes sign at [Fe/H] ≈ −1.5 dex.
In conclusion, while the behavior of RRc stars in Figure 3 suggests that RRc stars may not follow the same PLZ relation as RRab stars, the evidence is not strong; the sample of RRc stars we use is quite small (only 16 stars) and their parallaxes are fairly uncertain. In the future, we intend to reexamine the PLZ relations of RRc stars using more precise Gaia DR2 parallaxes and a larger set of RRc stars (e.g., the one published by Gavrilchenko et al. 2014). 9 The period slope is not expected to change significantly between the K and WISE bands (see Figure 4    The errobars indicate the quadratic sum of the uncertainty in the observed parallax and the inferred inverse distance. The RRab and RRc stars are shown as red and blue symbols, respectively. The RR Lyrae stars whose light curves are affected by the Blažko (1907) effect (i.e., the amplitude or phase modulation of the light curve; as identified via http://www.univie.ac.at/tops/ blazhko/Blazhkolist.html) are denoted with open circles, and the stars not known to exhibit the Blažko effect are denoted with solid circles. For RRc stars, we are showing fundamentalized periods (Equation 3). Note how the observed parallaxes of RRc stars seem to systematically deviate from the inferred inverse distances, suggesting that RRc stars may not follow the same PLZ relations as RRab stars. Also, note that Blažko-affected stars do not seem to scatter more than stars not known to exhibit the Blažko effect, indicating that the inclusion of Blažko-affected stars does not significantly increase the scatter in the PLZ relations (i.e., the σM parameter).

SUMMARY AND CONCLUSIONS
In this work, we have presented a probabilistic method that simultaneously constrains a period-luminosity-metallicity (PLZ) relation and validates (TGAS) parallax measurements. Compared to the traditionally-used weighted least-squares fitting of a PLZ relation, our approach allows for a direct usage of parallax (and other) measurements, while accounting fully for their precision. In comparison, the traditional approach cannot use imprecise parallaxes as their transformation into distance, and the subsequent characterization of the uncertainty in distance, is not a trivial task (Bailer-Jones 2015; Astraatmadja & Bailer-Jones 2016).
The final product of our approach is the full posterior distribution of all model parameters (in the form of a Markov chain). The posterior distribution enables a more general description of model parameters and their correlations (see Figure 2). For example, the marginal posterior distribution (i.e., the histogram) of σ M in Figure 2 clearly shows that this pa-rameter cannot be described by a Gaussian, which is something the traditional approach would be forced to do. Furthermore, the non-linear correlation between f and σ ,add is beautifully described by the joint posterior distribution of these two parameters (see Figure 2). The least-squares approach would be hard-pressed to capture such complexities.
We have used this probabilistic approach to constrain the PLZ relations in the near-IR W1 and W2 bands used by the WISE mission. Overall, the PLZ parameters we recover are consistent (within uncertainties) with the parameters found by previous studies (Madore et al. 2013;Klein et al. 2014;Dambis et al. 2014). Due to the fairly high fractional parallax uncertainty of our sample (median σ / ≈ 0.17), we were not able to constrain the slopes in log-period and metallicity more precisely than Dambis et al. (2014). However, since Dambis et al. (2014) used iterative least-squares fitting to constrain PLZ relations and fairly ad hoc removal of outlying data points, it is possible that their uncertainties may be underestimated. The second release of Gaia data, scheduled for April 2018, will provide more accurate and precise parallax measurements, and based on our current experience, is expected to place significantly tighter constraints on PLZ relations of RR Lyrae and other pulsating stars.
When fitting PLZ relations, we did not reject any outlying data points. Judging by the size of error bars and the distribution of − 1/r values of RRab stars in Figure 3, there does not seem to be many potential outliers that could have biased our results. However, as the precision of Gaia parallaxes improves, outliers may appear and unless they are properly handled, they may bias the measurement of PLZ relation parameters. Within the probabilistic framework presented in this work, inliers and outliers can be modeled using a mixture model (e.g., see Section 3 and Equation 17 of Hogg et al. 2010). As an example of how a mixture model can be used to model inliers and outliers in the context of constraining PLZ relations, we refer the interested reader to Appendix B of Sesar et al. (2016).
Since obtaining precise and accurate measurements of [Fe/H] for RR Lyrae stars is not a trivial task (e.g., see Nemec et al. 2013), the optimal data set for constraining PLZ relations for RR Lyrae stars may need to contain field and globular cluster RR Lyrae stars. By being at the same distances and by having the same (and precisely measured) [Fe/H], globular cluster RR Lyrae stars could be used to constrain the period and metallicity dependence of the PLZ relation (i.e., parameters a and b), while a few field RR Lyrae stars with well-measured Gaia parallaxes and [Fe/H] (e.g., those observed by Nemec et al. 2013) would constrain the zero-point of the PLZ relation. Such datasets were used by Sollima et al. (2006) and Dambis et al. (2014), but these studies used leastsquares fitting to constrain PLZ relations. In the future, we may apply our method to the same datasets.
While constraining PLZ relations, we also simultaneously constrained parameters that model TGAS parallax measurements and their uncertainties. To a precision of 0.05 mas, we did not find a statistically significant offset in TGAS parallaxes (i.e., the global offset parameter 0 = 0.00±0.05 mas) using our sample of distant RR Lyrae stars (median parallax of 0.8 mas and distance of 1.4 kpc). This result is consistent with the conclusion of Casertano et al. (2016), who use photometric parallaxes of distant Cepheids (median distance of ≈ 2 kpc) and find no offset to a precision of 0.02 mas.
The fact that we do not detect an offset in TGAS parallaxes may even be consistent with the findings of Stassun & Torres (2016) and Jao et al. (2016), who measure a global offset of ≈ −0.25 mas in TGAS parallaxes, but suggest that it may become negligible for parallaxes smaller than a few mas (i.e., at large distances). However, given the manner in which the trigonometric parallax measurements are made (linear offsets on the detector), we do not understand what physical effect could cause a distance-dependent offset in TGAS parallaxes.
Regarding the uncertainty in TGAS parallaxes, we find no need to rescale formal parallax uncertainties for RR Lyrae stars (i.e., no need for f > 1), and recommend the following equation when calculating their uncertainty in parallax σ ,RRLyr = (1.0ς ) 2 + (0.18/1000) 2 , where the formal parallax uncertainty ς can be calculated using Equation 2. The f = 1.1 and σ ,add = 0.12 mas values obtained by Gould et al. (2016) are consistent at the 1σ level with our findings. We emphasize that our and Gould et al. (2016) results for f and σ ,add were obtained using RR Lyrae stars. Due to as yet uncalibrated systematic effects in Gaia measurements, stars with different properties (e.g., color, brightness) may have different f and σ ,add values. Various systematic effects may also explain the values of f = 1.4 and σ ,add = 0.20 mas that Lindegren et al. (2016) adopted for TGAS. Unlike us, Lindegren et al. (2016) used a much more diverse sample of stars when constraining these two parameters (see their Appendix C.1).
By using individual measurements (e.g., TGAS parallaxes), PLZ relations, and parameters that model uncertainties in TGAS parallaxes, we have constrained distances to ≈ 120 RR Lyrae stars within 2.5 kpc of the Sun, to a 6% precision. While this precision may seem quite low compared to precisions reported in some previous studies (e.g., a 0.8% precision reported by Klein et al. 2014), we note that our estimate includes uncertainties due to correlations (e.g., between the period and metallicity slopes in the PLZ relation), and accounts for underestimated or unaccounted uncertainties in measurements and the model (via the σ M parameter). Studies that did not take such uncertainties into account, and did not propagate them properly through the model, may have overestimated the precision of their distance measurements.
In conclusion, we are looking forward to applying the method and the experience developed in this work to the next Gaia data release, and doing some exciting Galactic science with Cepheid, Mira, and RR Lyrae stars in the near future.
B.S. and H.-W.R. acknowledge funding from the European Research Council under the European Unions Seventh Framework Programme (FP 7) ERC Grant Agreement n. [321035]. We thank Dr. Adam G. Riess for the thorough review, positive comments, and constructive remarks on this manuscript. This project was de-veloped in part at the 2016 NYC Gaia Sprint, hosted by the Center for Computational Astrophysics at the Simons Foundation in New York City. This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/ gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa.int/ web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This publication makes use of data products from the Widefield Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.