© The Authors 2026

1 Introduction

Pulsating stars serve as crucial standard candles for the (extragalactic) distance scale thanks to the period-luminosity (PL) relations – Leavitt laws (Leavitt & Pickering 1912, LL) – they obey. The highest distance accuracy to date has been demonstrated for classical Cepheids, which can be observed in distant galaxies and have been extremely valuable for determining Hubble’s constant (e.g., Riess et al. 2022b). However, several types of population-II pulsating stars are available for distance determination in rather more nearby, old, and metal-poor populations. These include RR Lyrae (RRL) stars and type-II Cepheids (T2Cep), both of which are found in large numbers throughout the Milky Way halo, globular clusters (GCs), and nearby dwarf galaxies. Indeed, RRL were initially referred to as “cluster-type variables” (Smith 2000). Additionally, several other stellar types exhibit PLR that can be used for determining distances (Catelan & Smith 2015).

RR Lyrae are horizontal branch stars that typically pulsate in either the fundamental mode (RRab) or the first overtone (RRc), with typical periods of ∼0.6 d and ∼0.3 d, respectively. They are abundant: over 270 000 have been identified in the third data release of the ESA Gaia mission (Gaia Collaboration 2016, 2023; Clementini et al. 2023, GDR3), including in the Large Magellanic Cloud (LMC), Small Magellanic Cloud (SMC), M31, and numerous other Local Group systems (Clement 2017; Martínez-Vázquez 2023; Cruz Reyes et al. 2024). RRL exhibit LL in the infrared and in reddening-free magnitudes (e.g., Bhardwaj et al. 2017; Neeley et al. 2019; Braga et al. 2020). T2Cep occupy the classical instability strip at luminosities exceeding the horizontal branch and are frequently grouped into three subclasses, by increasing period, including the BL Her, W Vir, and RV Tau types. They are rarer but trace more advanced post-horizontal-branch evolutionary phases, including stars evolving toward or from the asymptotic giant branch (for a review, cf. Bono et al. 2024). In total, 1551 T2Cep have been identified in the gaiadr3.vari_cepheid table (Ripepi et al. 2023), and all T2Cep subtypes found in GCs exhibit a common LL (Cruz Reyes et al. 2025). At infrared wavelengths and Wesenheit magnitudes, T2Cep appear to extend the RRab LL to longer periods, which is particularly apparent in the Magellanic Clouds (Bhardwaj et al. 2017; Soszyński et al. 2019).

Several studies have presented calibrations of RRL and T2Cep LL based on GDR3 parallaxes in conjunction with infrared photometry or using Gaia photometry to compute reddening-free Wesenheit magnitudes, $M_G^W$. For example, Garofalo et al. (2022) used a hierarchical Bayesian method to calibrate the Gaia Wesenheit LL of field RRL, and Bhardwaj et al. (2023) used GC distances from Bailer-Jones et al. (2021) to calibrate RRL and T2Cep (separately) using near-infrared photometry. Both studies focused extensively on metallicity effects of RRL stars (Jurcsik & Kovacs 1996), which have been studied in a variety of photometric systems and using different approaches (e.g., Mullen et al. 2021, 2023; Li et al. 2023; Narloch et al. 2024). While previous calibrations did include GDR3 parallax corrections provided by Lindegren et al. (2021), they did not correct for well-known residual parallax systematics demonstrated for stars brighter than G ∼ 12 mag by a variety of methods (e.g., Khan et al. 2023, and references therein).

Clusters hosting RRL or T2Cep offer compelling opportunities for improvements over the state-of-the-art. Firstly, their average parallaxes can be determined with high accuracy using stars in the magnitude range not affected by the residual bias, in analogy with classical Cepheids in open clusters (Riess et al. 2022a; Cruz Reyes & Anderson 2023). Secondly, GCs hosting multiple specimens of a given type of pulsator provide valuable constraints on LL slopes. Thirdly, GCs hosting multiple pulsator types provide crucial constraints on the relative offsets between the RRab, RRc, and T2Cep LL intercepts. Fourthly, the second and third points also apply in the case when GC parallaxes are not determined with sufficient S/N for inversion, thus creating an opportunity for considering a much larger sample of GCs. Finally, the multiplicative effect on parallax induced by an additive offset in magnitude space (cf. Riess et al. 2021) allows to validate the accuracy of the GC parallaxes zero-point.

In this work, we have exploited these benefits by performing the first joint calibration of the LL of RRab, RRc, and T2Cep stars across 93 GCs. Our methodology employs a global least-squares framework inspired by the SH0ES distance ladder (Riess et al. 2009, 2022b) that relies exclusively on GDR3 data, notably trigonometric parallaxes and photometry. With this approach, we aim to

• calibrate the LL of RRab, RRc, and T2Cep while accounting for their correlated zero points and treating metallicity dependencies consistently;

• derive accurate distances to all the GCs in the sample and to assess the validity of concerns raised about Gaia parallaxes in the low-S/N regime (Vasiliev & Baumgardt 2021);

• provide a robust foundation for an independent Population-II distance ladder anchored directly to Gaia.

The structure of the paper is as follows. Section 2 presents the data considered, including a description of the GCs and variable star samples (Sects. 2.1 and 2.2), photometry and reddening (Sect. 2.3), and metallicity measurements (Sect. 2.4). Section 3 describes the methodology, including the models (Sect. 3.1), fitting approach (Sect. 3.2), MCMC framework (Sect. 3.3), outlier rejection (Sect. 3.4), and model variants (Sect. 3.5). Section 4 presents the main results, with baseline and variant fits. Section 5 discusses the implications and comparison with the literature. Finally, Section 6 summarizes our findings.

2 Dataset

This work exploits the astrometric and photometric capabilities of the Gaia mission (Gaia Collaboration 2016), a space-based survey by the European Space Agency designed to map over a billion stars in the Milky Way with high-precision parallaxes, proper motions, and multi-band photometry. Specifically, we use data from Gaia Data Release 3 (Gaia Collaboration 2023), which provides the most precise all-sky measurements to date for large samples of variable stars.

Building on the results of Cruz Reyes et al. (2024,2025), we incorporate average parallaxes for 170 GCs, allowing us to include RRL and T2Cep within these systems. This approach ensures a homogeneous, self-consistent dataset, enabling accurate luminosity and distance estimates while minimizing systematics.

Fig. 1 Mean residual parallax difference as a function of G-band magnitude, computed with rolling bins of width 0.5 mag and step size 0.05 mag. The x-axis marks the bin midpoints, while the y-axis gives the corresponding mean residual. Points are color-coded by the number of stars per bin, and the shaded area shows the 1σ uncertainty of the mean. The inset zooms into the range populated by most RRL and T2Cep stars, where a typical residual trend of ∼3−4 μas is seen.

2.1 Galactic globular clusters

Our sample is taken from the 170 GCs listed in Table A.2 of Cruz Reyes et al. (2024), which provides cluster-level parameters. In particular, it includes weighted mean parallaxes computed with the method of Maíz Apellániz et al. (2021), incorporating the Lindegren et al. (2021) parallax bias corrections and the covariance on small angular scales. Among these, 115 host RRL, 39 contain T2Cep, and 37 include both types of variables, cf. Sect. 2.2. These GCs offer a unique opportunity for calibrating distance scales, as their parallaxes enable direct geometric distance measurements, while their member standard candles provide a photometric approach. This combined approach makes GCs strong benchmarks for refining the distance ladder. They also allow for consistency checks between individual-star parallaxes and GC parallaxes, which should agree if no systematic biases are present.

2.1.1 Magnitude-dependent parallax systematics

Even after applying the Lindegren et al. (2021) correction to Gaia DR3 parallaxes, residual biases can remain. To assess their impact, we compare the parallaxes of individual cluster member stars to the weighted mean GC parallaxes from Cruz Reyes et al. (2024). Previous studies report that cluster parallaxes generally show minimal or no residual offset, even in the low signal-to-noise (S/N) regime (Baumgardt & Vasiliev 2021; Maíz Apellániz 2022; Riess et al. 2022a; Cruz Reyes et al. 2024).

While a direct comparison between individual stars and cluster parallaxes is limited by the low S/N of single-star parallaxes, aggregating stars across similar magnitudes enables a more meaningful test. We therefore apply a rolling binning scheme in G-band magnitude, similar to a moving average, to boost S/N and reveal any systematic trends. This analysis includes all 427 034 stars used in the GC parallax determinations.

Figure 1 shows the residual parallax offset as a function of Gaia G magnitude. In the key range G ∈ [13, 16], where most RRL and T2Cep reside, we detect a mild systematic trend of ∼3−4 μas. A sharper deviation appears near G ≈ 16 mag, coinciding with a change in Gaia’s CCD window class configuration (Gaia Collaboration 2016). However, this faint-end feature does not affect our analysis, as stars in this regime have low S/N and thus contribute little weight to the cluster parallax estimates.

2.2 Pulsating star samples

All the stars used in this study are members of GCs as identified by Cruz Reyes et al. (2024,2025), and are drawn from the Gaia DR3 catalog. The RRL sample is listed in Table 1 of Cruz Reyes et al. (2024), and the T2Cep sample in Table A.1 of Cruz Reyes et al. (2025). We retain only stars flagged as part of the final samples (Final=True) and with valid (non-NaN) values in the Posterior column, ensuring the exclusion of stars near GC centers where astrometric and photometric measurements are typically less reliable.

We further restrict the sample to stars with Specific Object Studies (SOS) photometry (SOS=True), a dedicated pipeline optimized for pulsating variables (Ripepi et al. 2019; Gaia Collaboration 2023; Clementini et al. 2023). Additional photometric quality cuts, detailed in Table 1, are applied. We also retain only stars for which the Lindegren et al. (2021) parallax correction is defined, and exclude stars with missing (NaN) values in parameters relevant to the analysis.

2.3 Photometry and Wesenheit magnitudes

It is well established that RRL do not follow a PLR in visual bands such as V, consistent with their characterization as horizontal branch stars (Catelan et al. 2004; Bhardwaj 2022). However, they follow a PLR at longer wavelengths or in reddeningfree bands such as optical Wesenheit magnitudes, where the inclusion of color terms restores a tight PL correlation (Marconi et al. 2015; Braga et al. 2015; Neeley et al. 2017).

In this work, we use Wesenheit magnitudes as defined in the Gaia photometric system (Madore 1982; Ripepi et al. 2019). Specifically, we adopt the optical Wesenheit magnitude $M_G^W$, constructed from the G, BP, and RP bands: $$m_G^W=m_G-R_G^W(m_{BP}-m_{RP}).$$(1) Here, m_G, m_BP, and m_RP are the observed magnitudes in the respective bands, and the reddening coefficient $R_G^W$ is defined as $$R_G^W=\frac{A_G}{A_{BP}-A_{RP}}.$$(2) We adopt a fixed value of $R_G^W=1.867$ for both RRL and T2Cep, following Cruz Reyes et al. (2024).

2.3.1 Uncertainties reassessment

Photometric uncertainties increase with magnitude, reflecting lower S/N for fainter stars. However, a small subset of stars in the Gaia SOS catalog report unrealistically small uncertainties compared to others of similar brightness (Cruz Reyes et al. 2024). For example, the star with source_id = 4295843576778048384 has a BP-band uncertainty (int_average_bp_error) of just 7 × 10⁻¹⁵ mag, orders of magnitude below typical values for that magnitude range. Since the Wesenheit magnitude propagates errors from multiple bands, an underestimated uncertainty in a single band can artificially reduce the total error, giving the star disproportionate weight in the fit.

To address this, we examine the typical photometric uncertainties as a function of magnitude and identify outliers. We then reassess each star’s reported uncertainty by comparing it to the median uncertainty of stars in the same 1 mag-wide bin. If a star’s uncertainty in any band is more than three standard deviations below the bin median, we replace it with the bin’s median value. Stars with typical or larger-than-expected uncertainties remain unchanged. This correction affects approximately 10% of stars, depending on filter and type of variable, and ensures that stars of similar brightness contribute comparably to the model fit.

Leavitt laws exhibit intrinsic scatter arising from astrophysical effects. The finite extent of instability strips in temperature causes a spread in luminosity at fixed period (e.g., Madore & Freedman 2012; Neeley et al. 2017; Bhardwaj et al. 2023; Khan et al. 2025). For RRL, an additional contribution arises from horizontal-branch (HB) evolution. As RRL evolve away from the zero-age horizontal branch (ZAHB) due to core-helium burning, they gradually change in luminosity and effective temperature. These evolutionary shifts move RRL in the period-magnitude plane and introduce an intrinsic broadening of the PLR, even in the NIR where temperature sensitivity is reduced (Bono et al. 2001; Catelan et al. 2004; Marconi et al. 2015; Neeley et al. 2019). To account for these combined effects, we include an intrinsic scatter term following the SH0ES approach for Cepheids (Riess et al. 2022b). Specifically, we add 0.045 mag in quadrature to the uncertainties of $m_G^W$ for both RRL and T2Cep. This scatter is significantly smaller than the intrinsic scatter observed in classical Cepheids (0.07 mag, cf. Table 1 in Khan et al. 2025, and references therein) and is consistent with the dispersion seen in our baseline fit (see Sect. 5.1).

2.4 Metallicities

Metallicity affects stellar luminosity and is a known systematic of LL that requires calibration. Metallicity effects on the LL of classical Cepheids have been studied intensively in the recent literature, leading to convergence on the picture that more metalrich classical Cepheids are more luminous (Romaniello et al. 2022; Breuval et al. 2024; Riess et al. 2022b; Trentin et al. 2024; Khan et al. 2025). Analogous investigations have been conducted for Population II pulsators, although the same level of consensus has not yet been reached. Numerous studies have investigated metallicity effects of RRL LL in a various photometric passbands using both empirical and theoretical approaches (Marconi et al. 2015; Neeley et al. 2019; Mullen et al. 2023; Bhardwaj et al. 2023; Narloch et al. 2024). The situation has thus far remained least clear for T2Cep, although Wielgórski et al. (2022) reported the first evidence of their metallicity dependence, in contrast to most other studies (e.g., Matsunaga et al. 2006; Feast et al. 2008; Bhardwaj et al. 2017; Ngeow et al. 2022; Ripepi et al. 2023; Sicignano et al. 2024; Narloch et al. 2025).

For T2Cep, we adopt the metallicity of their host cluster, as listed in the Harris (2010) catalog. These values are based on spectroscopic measurements of red giant branch stars and serve as fiducial [Fe/H] estimates for Milky Way GCs. This approach assumes chemical homogeneity within clusters and does not account for intra-cluster variation.

The same catalog values can also be used for RRL. However, RRL light curves correlate with metallicity (Jurcsik & Kovacs 1996), enabling photometric [Fe/H] estimates derived from their pulsation morphology. We adopt such light-curve-based proxies as relative indicators of chemical composition. While the absolute accuracy of these calibrations may be uncertain, this has limited impact on our analysis: PL-based distances are primarily sensitive to relative metallicity differences. Since the calibration is applied homogeneously across our sample, it still allows robust measurement of the PLR’s metallicity dependence without introducing systematic distance offsets.

For Gaia G-band SOS light curves (Gaia Collaboration 2023; Clementini et al. 2023), we adopt the recent calibrations by Li et al. (2023), derived from a sample of over 2700 RRL with spectroscopic metallicities. The P−φ₃₁−R₂₁−[Fe/H] relation is used for RRab stars, and the P−R₂₁−[Fe/H] relation for RRc stars, as defined below: $$\begin{array}{rl}[\mathrm{F}\mathrm{e}/\mathrm{H}{]}_{RRab}=& -1.888-5.772\cdot (P-0.6)\\ & +1.065\cdot (R_{21}-0.45)+1.090\cdot ({\varphi }_{31}-2)\end{array}$$(3) $$\begin{array}{rl}[\mathrm{F}\mathrm{e}/\mathrm{H}{]}_{RRc}=& -1.737-9.968\cdot (P-0.3)\\ & -5.041\cdot (R_{21}-0.20).\end{array}$$(4)

To assess the impact of systematics, we also consider an alternative calibration by Muraveva et al. (2025). Although the functional form differs from that of Li et al. (2023), it is similarly based on the idea that the shape of RRL light curves encodes information about metallicity. Specifically, this calibration employs different combinations of light-curve parameters: $${[\mathrm{Fe}/\mathrm{H}]}_{RRab}=-0.37-5.55·P+0.94·{\varphi }_{31}$$(5) $${[\mathrm{Fe}/\mathrm{H}]}_{RRc}=0.80-8.54·P+0.23·{\varphi }_{31}-9.34·A_2.$$(6) These calibrations are particularly advantageous as they are derived from a large homogeneous dataset using Gaia photometry, allowing us to estimate metallicities homogeneously for the RRL in our sample, many of which lack direct spectroscopic measurements. This approach enables us to account for both inter- and intra-cluster metallicity variations while minimizing systematic biases.

3 Methodology

Our method simultaneously calibrates the PLR for RRL and T2Cep by leveraging the precise average parallaxes of GCs. The objective is to jointly fit the various PLR (RRab, RRc, T2Cep) across different clusters. Any cluster containing multiple stars of a given type (e.g., RRab) contributes information on the slope of the corresponding PLR, regardless of its distance precision. Moreover, GCs hosting multiple types of pulsating stars provide valuable constraints on the relative luminosity scales between the different candles. To ensure consistency and precision in the calibration, we adopt the methodology employed in the Cepheid distance ladder, following the approach developed by the SHOES team (Riess et al. 2016, 2022b).

This method requires working in magnitude space, which involves inverting parallaxes, a process that can introduce biases when parallax uncertainties are large. To address this, we classified clusters into two categories based on the quality of their parallaxes. Clusters with reliable parallaxes (typically nearer, with low extinction and high S/N) were designated as anchors, as they provide the geometric foundation for anchoring the LL. More distant clusters with less reliable parallaxes or significant extinction are designated as hosts. Their distances are not directly constrained by their parallaxes but are instead inferred from the fitted PLR, allowing us to bypass the use of their less reliable parallax measurements.

This two-class classification enabled us to incorporate variable stars from a wide range of GCs. Even when their geometric distances are uncertain, these clusters still contribute essential information, particularly in constraining the slopes of the PLR. The criteria used to classify clusters as anchors or hosts are summarized in Table 1.

3.1 Models

For each type of variable, a separate PLR is fitted. The model differs slightly depending on whether the star belongs to a host or an anchor cluster. For stars in hosts, where the distance modulus μ_{0, i} is a free parameter, the PLR for RRab stars is $$\begin{array}{rl}{m}_{i,j}^W={\mu }_{0,i}& +{\alpha }_{0,\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}\\ & +{\beta }_{0,\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}(\log P_{i,j}-\log P_{0,\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}})\\ & +{\alpha }_{M,\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}([\mathrm{F}\mathrm{e}/\mathrm{H}{]}_{i,j}-[\mathrm{F}\mathrm{e}/\mathrm{H}{]}_{0,\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}).\end{array}$$(7) Here, indices i and j refer to the cluster and the individual star within that cluster, respectively. The quantity $m_{i,j}^W$ denotes the observed Wesenheit magnitude of star j in cluster i, and μ_0,i is the extinction-corrected distance modulus of that cluster. The term α_0,RRab represents the absolute Wesenheit magnitude at the pivot period P_0,RRab and metallicity [Fe/H]_0,RRab, while β_0,RRab and α_{M, RRab} are the slopes with respect to log P and [Fe/H].

For RRc and T2Cep, we adopt the same model structure. However, we calibrated their LL intercepts relative to the fiducial magnitude of the RRab stars: $$\begin{array}{rl}{m}_{i,j}^W={\mu }_{0,i}& +{\alpha }_{0,\text{\hspace{0.17em}RRab\hspace{0.17em}}}+\mathrm{\Delta }{\alpha }_{0,\text{\hspace{0.17em}RRc\hspace{0.17em}}/\text{\hspace{0.17em}T2Cep\hspace{0.17em}}}\\ & +{\beta }_{0,\text{\hspace{0.17em}RRc\hspace{0.17em}}/\text{\hspace{0.17em}T2Cep\hspace{0.17em}}}(\log P_{i,j}-\log P_{0,\text{\hspace{0.17em}RRc\hspace{0.17em}}/\text{\hspace{0.17em}T2Cep\hspace{0.17em}}})\\ & +{\alpha }_{M,\text{\hspace{0.17em}RRc\hspace{0.17em}}/\text{\hspace{0.17em}T2Cep\hspace{0.17em}}}([\mathrm{F}\mathrm{e}/\mathrm{H}{]}_{i,j}-[\mathrm{F}\mathrm{e}/\mathrm{H}{]}_{0,\text{\hspace{0.17em}RRc\hspace{0.17em}}/\text{\hspace{0.17em}T2Cep\hspace{0.17em}}}).\end{array}$$(8) Hence, the offsets Δ α_0,RRc or Δ α_{0, T2Cep} represent the differences between the fiducial RRab intercept and the fiducial LL intercepts of RRc and T2Cep, respectively, each defined at their own reference period and metallicity. This formulation enables a consistent calibration across all pulsator classes while anchoring the absolute scale to RRab stars.

With this structure, the global fit is anchored by a single absolute zero point (α_0,RRab), while relative offsets (Δ α_0,RRc/T2Cep) capture the luminosity differences between classes. If each class were instead fitted with its own independent zero point, all three calibrations would shift coherently with any change in the anchor distances, leading to strong correlations between the zero points and the GC distance moduli. Using one common α₀ absorbs this global correlation, while the relative offsets are sensitive only to differences between classes and remain largely decoupled from the distance scale. In practice, this parametrization improves consistency across classes and confines the dominant uncertainty to α_0,RRab, while Δ α₀ is constrained primarily by clusters that host multiple types of variables.

We adopt pivot values P₀ and [Fe/H]₀ near the median period and metallicity of each class to minimize covariances between fit parameters and to improve numerical stability. These values are chosen close to the sample averages and rounded for convenience, facilitating the interpretation of the fitted slopes and intercepts.

For anchors, the equations take a slightly different form: the geometric distance modulus appears on the left-hand side (and is therefore not a fitted parameter), while a cluster-specific correction term Δ μ_{0, i} is introduced on the right-hand side to account for residuals relative to the geometric distance: $$\begin{array}{rl}{m}_{i,j}^W-{\mu }_{0,i}=\mathrm{\Delta }{\mu }_{0,i}& +{\alpha }_{0,\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}\\ & +{\beta }_{0,\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}(\log P_{i,j}-\log P_{0,\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}})\\ & +{\alpha }_{M,\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}([\mathrm{F}\mathrm{e}/\mathrm{H}{]}_{i,j}-[\mathrm{F}\mathrm{e}/\mathrm{H}{]}_{0,\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}})\end{array}$$(9) $$\begin{array}{rl}{m}_{i,j}^W-{\mu }_{0,i}=\mathrm{\Delta }{\mu }_{0,i}& +{\alpha }_{0,\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}+\mathrm{\Delta }{\alpha }_{0,\mathrm{R}\mathrm{R}\mathrm{c}/\mathrm{T}2\mathrm{C}\mathrm{e}\mathrm{p}}\\ & +{\beta }_{0,\mathrm{R}\mathrm{R}\mathrm{c}/\mathrm{T}2\mathrm{C}\mathrm{e}\mathrm{p}}(\log P_{i,j}-\log P_{0,\mathrm{R}\mathrm{R}\mathrm{c}/\mathrm{T}2\mathrm{C}\mathrm{e}\mathrm{p}})\\ & +{\alpha }_{M,\mathrm{R}\mathrm{R}\mathrm{c}/\mathrm{T}2\mathrm{C}\mathrm{e}\mathrm{p}}([\mathrm{F}\mathrm{e}/\mathrm{H}{]}_{i,j}-[\mathrm{F}\mathrm{e}/\mathrm{H}{]}_{0,\mathrm{R}\mathrm{R}\mathrm{c}/\mathrm{T}2\mathrm{C}\mathrm{e}\mathrm{p}}).\end{array}$$(10) The additional term Δ μ_{0, i} captures small residual offsets between the model-predicted and geometric distance moduli. This term is constrained by $$0=\Delta {\mu }_{0,i}\pm \sigma \left({\mu }_{0,i}\right),$$(11) where σ(μ_{0, i}) is the uncertainty on the parallax-based distance modulus. This treatment, following Riess et al. (2022b), allows the use of multiple GC anchors to geometrically calibrate the global fit.

The equations are formulated in magnitude space, whereas Gaia provides parallaxes. This requires careful treatment when converting between distances and magnitudes. Estimating distances via d=1/ϖ becomes increasingly biased when the relative parallax uncertainty exceeds 0.1 (Luri et al. 2018). To mitigate this, we apply an S/N cut (see Table 1) when selecting anchors, ensuring that only parallaxes with sufficient precision are used. This allows us to safely invert the parallaxes without introducing significant bias. The corresponding transformation to distance modulus is $$\mu =5log\left(d\right)-5=5log\left(\frac{1}{\varpi }\right)-5,$$(12) where μ is the distance modulus in [mag], d the distance in [pc], and ϖ the parallax in [arcsec].

3.2 Joint fit

The full definitions of the matrices and the explicit derivation of the least-squares solution are provided in Appendix A. In short, the fit reduces to solving a system of linear equations, where each star and each geometric constraint contributes one relation. These can be expressed in matrix form as $$\mathit{y}=\mathit{L}·\mathit{q},$$(13) with y denoting the left-hand side of Eqs. (7)–(11), L the design matrix containing the known coefficients from the right-hand side, and q the vector of parameters to be estimated. In the baseline case, L is a large matrix (1207 × 102), but its block structure allows for an intuitive representation of the problem.

3.3 Markov chain Monte Carlo

The Markov chain Monte Carlo (MCMC) approach extends the least-squares solution by exploring the full posterior distribution of the parameters. This provides not only central values but also a comprehensive assessment of the associated uncertainties. Specifically, we infer the posterior distribution p(q ∣ L, C, y) of the model parameters q given the observational data y, the design matrix L, and the associated uncertainties in C: $$p(\mathbf{q}\mid \mathbf{L},\mathbf{C},\mathbf{y})\propto p\left(\mathbf{q}\right)p(\mathbf{y}\mid \mathbf{L},\mathbf{C},\mathbf{q}).$$(14) Here p(q) denotes the prior distribution, and p(y ∣ L, C, q) is the likelihood. Assuming a Gaussian likelihood with diagonal covariance $\mathbf{C}=diag({\sigma }_1^2,...,{\sigma }_n^2)$, we write $$p(\mathbf{y}\mid \mathbf{L},\mathbf{C},\mathbf{q})={\prod }_i\frac{1}{\sqrt{2\pi }{\sigma }_i}exp(-\frac{{(y_i-y_{\mathrm{fit},i})}^2}{2{\sigma }_i^2}),$$(15) which leads to a log-likelihood of the form $$lnp(\mathbf{y}\mid \mathbf{L},\mathbf{C},\mathbf{q})\propto -\frac{12}{}{\sum }_i{\chi }_i^2,$$(16) with $${\chi }^2={(\mathbf{y}-\mathbf{Lq})}^T{\tilde{\mathbf{C}}}_L^{-1}(\mathbf{y}-\mathbf{Lq}),$$(17) where C̃_L is the Cholesky decomposition of the diagonal covariance matrix C. Although C⁻¹ is analytically straightforward, the Cholesky decomposition improves numerical stability, particularly for large or ill-conditioned systems.

For the prior p(q), we adopt a uniform distribution centered on the analytic least-squares solution μ from Eq. (A.10): $$q_i\sim Uniform({\mu }_i-10{\sigma }_i,{\mu }_i+10{\sigma }_i).$$(18) A broad prior (20σ width) that allows the chains to explore the parameter space freely while remaining informed by the analytic solution.

Sampling was performed with the emcee Python package (Foreman-Mackey et al. 2013), using 500 walkers initialized uniformly within the prior bounds. Convergence was monitored via the autocorrelation time τ, estimated following Goodman & Weare (2010). A burn-in phase of 5 τ is discarded, after which the chains are extended to a total length of ∼100 τ, ensuring adequate sampling of the posterior. The final posteriors lie well within the central region of the prior range, confirming that the prior bounds do not bias the inferred parameters.

3.4 Outlier rejection

To mitigate the impact of potential outliers that pass through the selection criteria, an outlier-rejection algorithm was implemented within the fitting routine. The procedure follows the standard σ-clipping approach, in which the most deviant point is removed at each iteration until no star deviates by more than κ σ from the fit (e.g., Kodric et al. 2015).

Because the dataset contains multiple types of variable stars with distinct intrinsic scatter, the classical single-population σ-clipping method is not directly applicable. Both the dispersion σ and the rejection threshold κ differ between stellar populations, making a uniform rejection criterion unsuitable. To account for this, a modified algorithm (see Appendix B, Alg. 1) is adopted that incorporates population-dependent scatter. At each iteration, the most deviant star relative to the scatter of its own population is removed, and the process continues until the stopping criterion for each population is satisfied.

The values of κ_i for different populations can be freely chosen. Unless otherwise stated, they are set using Chauvenet’s criterion for each population.

3.5 Baseline fit and variants

To assess the robustness of the results, a baseline model is defined that incorporates the full dataset and reflects the most comprehensive set of assumptions. This model serves as the primary reference throughout the analysis. Its configuration is summarized in Table 2.

To evaluate the sensitivity of the conclusions to specific modeling choices and potential sources of systematic error, a series of targeted variants is constructed. Each variant isolates the impact of a particular assumption, input, or methodological choice. The variants are grouped into categories, each testing different sources of systematics in the analysis.

1. Baseline: the primary reference model.

2. Outlier rejection: sensitivity of the baseline model to different outlier-rejection criteria.

   • No OR: no outlier rejection.

   • κ=5.0: weak rejection with κ=5.0.

   • κ=3.5: same as above, with κ=3.5.

   • κ=3.0: same as above, with κ=3.0.

   • κ=2.5: strict rejection with κ=2.5.

3. Metallicities: impact of different metallicity calibrations.

   • Muraveva: calibration from Muraveva et al. (2025) for RRab and RRc.

   • Cluster: GC metallicities from Harris (2010) for all candles.

4. Parallaxes: alternative determinations of GC parallaxes.

   • Median: median parallaxes instead of weighted mean.

   • G<19: weighted mean parallaxes using only stars with G<19 mag.

   • 14.5<G<19: weighted mean parallaxes using only stars with 14.5<G<19.

   • VB21: parallaxes from Table A.1 of Vasiliev & Baumgardt (2021) (VB21).

5. Models: alternative assumptions to Eqs. (7)–(10).

   • No RRab: excluding RRab.

   • No RRc: excluding RRc.

   • No T2Cep: excluding T2Cep.

   • β_0,RRc=β_0,RRab: common slope β for RRab and RRc.

   • β_0,T2Cep=β_0,RRab: common slope β for RRab and T2Cep.

   • β_0,T2Cep=β_0,RRc: common slope β for RRc and T2Cep.

   • β_{0, T2 Cep}=β_0,RRc=β_0,RRab: all types share the same slope β.

   • α_{M, RRc}=α_{M, RRab}: common metallicity dependence α_M for RRab and RRc.

   • No α_M: metallicity dependence ignored (α_M=0).

   • β₀ → β₀+β_M[Fe/H]: period slope β depends on metallicity.

   • β₀ → β₀+β_M[Fe/H]+ No α_M: same as above, but no metallicity effect on the intercept.

   • β → β₀+β_M[Fe/H]_Cl: same as above, with cluster metallicities from Harris (2010).

   • β → β₀+β_M[Fe/H]_Cl+ No α_M: same as above, but no metallicity effect on the intercept, using cluster metallicities.

6. Pre-processing: sensitivity to pre-processing parameters.

   • σ_extra=0: no additional scatter from the instability strip.

   • σ_extra=0.025: extra scatter decreased by 50% relative to baseline.

   • σ_extra=0.075: extra scatter increased by 50% relative to baseline.

7. Anchor criteria: alternative thresholds for identifying GC anchors.

   • S/N=5: threshold set to half the baseline value.

   • S/N=20: threshold set to twice the baseline value.

   • E(B−V)=∞: no extinction-based criterion.

   • E(B−V)=0.5: threshold set to half the baseline value.

4 Results

This section presents the LL calibrations from the baseline analysis (Sect. 4.1), with uncertainties characterized using the MCMC approach (Sect. 4.1.1). The corresponding distance estimates to GCs are given in Sect. 4.1.2, while Sect. 4.2 examines the impact of model variants. The complete set of best-fit parameters, including both the baseline and all variant models, is listed in Table E.1. The distances to anchor and host GCs for the baseline configuration are reported in Tables E.2 and E.3.

4.1 Baseline results

Figure 2 presents the data and fit from the baseline analysis and illustrates the residual scatter for RRab, RRc, and T2Cep. The absolute magnitudes of all stars are standardized for metallicity differences to properly reflect the scatter in the residuals. The three types of variables occupy distinct period ranges with minimal overlap. The bottom panel shows the residuals, where RRab and T2Cep exhibit larger scatter (σ=0.085 mag and σ= 0.141 mag, respectively) compared to RRc (σ=0.050 mag). The fit yields χ² ≈ 1089 for 1101 degrees of freedom, indicating good agreement between model and data. During the fitting process, Chauvenet’s criterion results in threshold values of κ=3.44 for RRab, 3.20 for RRc, and 2.26 for T2Cep. Based on these thresholds, the outlier-rejection algorithm excludes about 5% of the total sample: 43 RRab, 21 RRc, and no T2Cep.

Equations (19)–(21) provide the Gaia-based Wesenheit absolute-magnitude calibrations derived from the baseline fit for RRab, RRc, and T2Cep. The fiducial magnitudes for RRc and T2Cep are computed from the fitted Δ α₀ values (cf. Eq. (9)), with uncertainties added in quadrature (cf. Table E.1): $$\begin{array}{rl}{M}_{\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}=& (-0.300\pm 0.020)\\ & +(-2.582\pm 0.047)\cdot (\log P+0.25)\\ & +(-0.028\pm 0.008)\cdot ([\mathrm{F}\mathrm{e}/\mathrm{H}]+1.5),\end{array}$$(19) $$\begin{array}{rl}{M}_{\mathrm{R}\mathrm{R}\mathrm{c}}=& (0.004\pm 0.021)\\ & +(-2.497\pm 0.095)\cdot (\log P+0.5)\\ & +(0.013\pm 0.010)\cdot ([\mathrm{F}\mathrm{e}/\mathrm{H}]+1.5),\text{and}\end{array}$$(20) $$\begin{array}{rl}{M}_{\mathrm{T}2\mathrm{C}\mathrm{e}\mathrm{p}}=& (-3.407\pm 0.027)\\ & +(-2.506\pm 0.039)\cdot (\log P-1)\\ & +(0.107\pm 0.043)\cdot ([\mathrm{F}\mathrm{e}/\mathrm{H}]+1.5).\end{array}$$(21)

The absolute magnitude of the fiducial candle (α_0,RRab) is calibrated to a precision of 0.020 mag (0.92% in distance). RRc and T2Cep reach comparable accuracy, with total uncertainties of 0.021 mag(0.97%) and 0.027 mag(1.3%), respectively. These values match the precision of Gaia DR3-based calibrations for classical Cepheids (0.9%; Cruz Reyes & Anderson 2023). Among the three populations, RRc show the lowest intrinsic scatter, indicating that they may offer the highest precision under ideal conditions. However, the larger amplitudes and distinctive light curves of RRab make them more easily identifiable and therefore more practical for distance measurements.

The three PL slopes are broadly consistent across types, though not identical: RRab stars have a slope of −2.582 ± 0.047 mag, while RRc and T2Cep yield −2.497 ± 0.095 mag and −2.506 ± 0.039 mag, respectively. A small but statistically significant metallicity dependence is detected for RRab (α_M= −0.028 ± 0.008 mag/dex), whereas RRc show a negligible effect (α_M=0.013 ± 0.010 mag/dex). In contrast, T2Cep display a more pronounced dependence of 0.107 ± 0.043 mag/dex, albeit at lower significance.

These baseline results represent the most accurate absolute calibrations to date for RRab, RRc, and T2Cep stars in the Gaia Wesenheit system. A direct comparison with the RRab calibration from Garofalo et al. (2022) is presented in Sect. 5.2.

Fig. 2 Baseline joint fit to fundamental-mode, first-overtone RR Lyrae, and type-II Cepheids in Galactic GC. Colors indicate different stellar types: blue for RRab, green for RRc, yellow for T2Cep, and red for stars excluded by the outlier rejection. Circles mark stars in hosts, while downward-pointing triangles denote stars in anchors. Top: absolute magnitudes, corrected for metallicity, plotted against log P. For hosts, $M_{i,j}^W=m_{i,j}^W-{\mu }_{0,i}-{\alpha }_M({[\mathrm{Fe}/\mathrm{H}]}_{i,j}-{[\mathrm{Fe}/\mathrm{H}]}_0)$; for anchors, an additional −Δ μ0, i term applies. Colored dashed lines show the fitted PLR. Bottom: residuals between the corrected magnitudes and the fitted PLR, with dashed lines marking the residual standard deviation. The inset (bottom right) shows the residual distributions.

4.1.1 Uncertainties and Markov Chain Monte Carlo

Because the problem is linear, the least-squares matrix formalism provides an adequate solution. Its main limitation is the requirement of a full covariance matrix, which in this case is simplified to a diagonal form C (Sect. 3). This assumption neglects potential correlations between data points (off-diagonal terms) but ensures analytical consistency under linearity and Gaussian-distributed errors. To test whether this simplification biases the results or underestimates the uncertainties, an MCMC analysis is carried out as an independent cross-check. The posterior distributions obtained with MCMC show excellent agreement with the baseline least-squares solution. Their near-Gaussian shape and the close match in parameter uncertainties confirm that the model is well calibrated and that the error budget is robust. Together, these results demonstrate that the least-squares estimates provide a reliable and efficient characterization of the fit.

Posterior distributions for selected parameters are shown in Fig. C.1, together with direct comparisons to the analytical (least-squares) estimates. The figure displays the primary fit parameters for RRab and RRc stars (results for T2Cep are similar), as well as representative examples of distance parameters for one host and one anchor cluster. The agreement is excellent; for example, the absolute calibration of RRab in the baseline fit is α_0,RRab,LS=−0.300 ± 0.020 mag from the least-squares method and α_0,RRab,MCMC=−0.299 ± 0.020 mag from the MCMC analysis. Across all parameters, differences between the two methods remain within 0.1σ, confirming that the least-squares approach neither introduces significant bias nor underestimates uncertainties.

Most parameters show minimal mutual correlation. A strong anti-correlation is evident between the absolute calibration parameter α_0,RRab and the distance-related terms μ and Δ μ, as expected: increasing the inferred distances implies intrinsically brighter stars, requiring lower α_0,RRab values to remain consistent with the observed magnitudes. Mild correlations also appear between the period slopes (e.g., β_RRab) and their associated metallicity terms (e.g., α_{M, RRab}), suggesting that these parameter pairs encode partially overlapping physical information. This motivates the inclusion of cross-terms such as β_M ·[Fe/H] · log P in some model variants (see Sect. 3.5).

Overall, the MCMC analysis confirms the robustness of the least-squares solution. The close agreement in both central values and uncertainties demonstrates that simplifying the covariance matrix to a diagonal form does not bias the results. The near-Gaussian posteriors indicate that the parameters are well constrained and free from strong degeneracies, while the consistency of the covariance estimates at the few-percent level provides confidence that the error budget is reliable. Together, these results confirm that the least-squares solution is unbiased, efficient, and statistically consistent with the full posterior distribution explored via MCMC.

Fig. 3 Visualization of the baseline mini-distance ladder. The lower-left panel shows anchor GC parallaxes on the x-axis against the joint-fit distances on the y-axis, including residuals. The mean difference between Gaia parallax and fit parallax for anchor GCs is −0.87 ± 1.52 μ as. The upper-right panel shows the host GC fit distances on the x-axis against the host GC parallaxes on the y-axis, with a mean difference of −5.9 ± 4.3 μ as affected by outliers. The two outlier hosts Terzan 1 and BH 229 are due to sparse membership and do not affect the fit. The comparison reveals no evidence for a distance dependent bias of GDR3 parallaxes.

4.1.2 Homogeneous distances to 93 GCs

A key result of the baseline analysis is a homogeneous determination of GC distances, calibrated to a zero point set by Gaia DR3 parallaxes. For anchors, the distances correspond to the geometric parallax-based estimates, modified by the fitted correction terms (Δ μ) from Eq. (11). For hosts, the distance moduli (μ) are inferred directly from the fit. Overall, Gaia parallaxes are consistent with the best-fit distances, aside from two identifiable outliers discussed below. Figure 3 illustrates the agreement between the Gaia-based and jointly fitted distances, highlighting their overall consistency.

For anchors, the methodology allows a direct comparison between fitted distances and Gaia-based values through the correction factor Δ μ from Eq. (11). This comparison assesses how well the fitted Δ μ values remain consistent with zero, within the uncertainties of the Gaia parallax measurements for GC.

This comparison is illustrated in Fig. 4, where both the initial and fitted values of Δ μ are shown for each anchor, ordered by distance. As expected, the uncertainties increase with distance for both the Gaia parallaxes and the fitted correction parameters. Statistically, about 70% of GC have correction factors within 1σ of their Gaia values, and all anchors fall within 2σ of their starting value.

A significant outlier appears in the case of NGC 6553, with a fitted Δ μ offset nearly 6σ from zero. Closer inspection shows that this discrepancy arises from a single RRc candidate (source_id = 4064796635702285056). Although the star is confirmed to be an RRc, its membership probabilities, both prior and posterior, indicate that it is unlikely to belong to NGC 6553. This star is therefore excluded from the final dataset, and since no other valid variables remain in that cluster, NGC 6553 is omitted from the final results.

Overall, the distances derived from the fit closely match those obtained from Gaia parallaxes for the anchors. This agreement is corroborated in the bottom-left panel of Fig. 3, where the mean parallax difference between Gaia measurements and the joint-fit values is −0.87 μ as, with a standard deviation of 9.26 μ as. A similar consistency is observed for the hosts in the upper-right panel of the same figure, with a mean difference of −5.89 μ as and a standard deviation of 31.84 μ as.

Caution is warranted when interpreting the comparison for the hosts, as the reliability of their Gaia parallaxes decreases at larger distances. Terzan 1 (HP2) and BH 229 (HP1), which appear as the two upper-left points in the plot, contain only 40 and 54 identified member stars, respectively. In such sparsely populated cases, membership assignments are less robust and may lead to biased cluster parallaxes. These GC are hosts rather than anchors, so their Gaia parallaxes are not used in the fit. Their uncertain parallaxes therefore do not affect the derived distances, which are determined from the fit itself. For example, the Gaia parallax for Terzan 1 is about 25 μas, whereas both the fit and independent literature estimates give values in the 175-200 μas range (Baumgardt & Vasiliev 2021).

Fig. 4 Correction parameters Δ μ for each anchor GC in the baseline. Blue symbols show the initial values (set to zero; Eq. (11)), and red symbols the fitted values (Eq. (A.10)). Clusters are ordered by increasing distance, so the typical uncertainty in Δ μ grows along the x-axis, although the number of member stars also contributes. Overall, 70% of anchors remain within 1σ of their Gaia parallax values, and all lie within 2σ.

4.2 Variants

A comprehensive visual summary of the variant fits is provided in Appendix D (Fig. D.1), while the numerical results are compiled in Table E.1. In total, 30 model variants, defined and categorized in Sect. 3.5, are fitted to assess the impact of specific assumptions, data subsets, and potential sources of systematics.

In brief, all variants confirm the robustness of the baseline solution. The fitted parameters remain consistent, with nearly all deviations lying within 1σ of the baseline values. This demonstrates that the absolute calibration is stable against changes in modeling assumptions or sample definitions. The whisker plots in the appendix further illustrate this point by showing both the statistical uncertainties of each variant and the narrow clustering of results around the baseline calibration.

4.2.1 Outlier rejection

The fit is remarkably robust once some level of outlier rejection is applied. The calibration stabilizes quickly for moderate thresholds and becomes effectively unchanged (<1σ deviation) for stricter criteria. The specific value of the threshold κ has only a minor effect on the final results, but the presence of an outlier-rejection step is essential. Omitting this step entirely introduces significant biases in several parameters; for instance, β_RRc shifts by nearly 5σ relative to the baseline model. In contrast, any variant that applies even a weak rejection (e.g., κ=5.0) yields results that remain within 1σ of the baseline.

As expected, the choice of κ has a pronounced impact on the goodness-of-fit metric. Without outlier rejection, the reduced chi-squared increases sharply to ${\chi }_{\text{dof}}^2=6.297$. A weak rejection threshold (κ=5.0) reduces it to ${\chi }_{\text{dof}}^2=1.410$, while a strict cutoff (κ=2.5) results in ${\chi }_{\text{dof}}^2=0.497$. These results confirm that the rejection step is critical to ensure a robust calibration, while the precise value of κ has limited influence. Chauvenet’s criterion provides a principled, data-driven threshold that balances robustness with minimal distortion of the underlying distribution.

Figure D.1 reveals a subtle but systematic trend: the metallicity term α_{M, RRab} decreases as the rejection threshold becomes stricter, which in turn yields slightly brighter absolute calibrations α_0,RRab. This behavior is expected, since both metallicity and outliers contribute to scatter in the PLR fit. Stricter rejection may inadvertently exclude stars with genuine residuals driven by metallicity effects. Overly aggressive rejection can therefore suppress intrinsic scatter associated with metallicity, artificially flattening α_M and biasing the calibration. These trends underscore the need for an objective, statistically motivated threshold. In this analysis, Chauvenet’s criterion provides such a compromise: it removes true outliers while minimizing the risk of discarding stars whose deviations arise from astrophysical causes such as metallicity variation. Moreover, a 5σ threshold yields results nearly identical to baseline, whereas omitting rejection shifts the absolute calibration by 1σ, indicating that the outliers affecting α_M lie beyond 5σ from the fit.

4.2.2 Metallicities – different sources for [Fe/H]

Replacing the photometric metallicity calibration with that of Muraveva et al. (2025) reduces the RRc sample by roughly half because the light-curve parameters φ₃₁ and A₂ are frequently unavailable, especially in the hosts. This variant increases the uncertainty on α_M for RRc, although the results remain consistent within 1σ.

Photometric metallicity estimates do not necessarily provide absolute values of [Fe/H], but rather proxies derived from light-curve morphology (via Fourier parameters). In practice, this means that the fit captures a correlation with light-curve shape rather than metallicity on an absolute scale. The absolute calibration of [Fe/H] is therefore not critical for determining GC distances. However, because these calibrations are light-curve dependent, any imperfections may induce a residual correlation between period and metallicity effects.

In the variant using GC metallicities from Harris (2010), the results again remain consistent with the baseline. The uncertainties on the metallicity coefficients α_M increase significantly, by more than a factor of five for both RRab and RRc. This outcome is expected: GC metallicities are determined from red giant branch stars and may not capture variations among RRL. Assigning a single metallicity per cluster removes any intracluster variation that would otherwise help constrain α_M, thereby reducing the leverage to measure metallicity dependence. In both cases, the absolute calibrations of the standard candles remain largely unchanged.

Interestingly, the α_M term for T2Cep shifts by ∼2σ in this variant, even though the metallicity calibration was only altered for RRab and RRc. This behavior suggests that the metallicity dependence inferred for T2Cep is largely driven by the more numerous RRab, through their weight in the global fit. The change also reflects the rejection of a few additional T2Cep once cluster metallicities were adopted, further amplifying the effect. While the shift remains modest, it highlights that the metallicity dependence of T2Cep should be treated with caution.

4.2.3 Parallaxes

The two variants that restrict the sample to narrower apparent-magnitude ranges (G<19 and 14.5<G<18) yield results nearly identical to those of the baseline, despite visible inconsistencies between faint-star parallaxes and GC-level parallaxes (see Fig. 1, where 〈ϖ_i−〈ϖ〉_GC〉 reaches up to −40 μ as). This minimal impact is expected, as faint stars typically carry large uncertainties and therefore contribute little statistical weight to the GC parallax estimates.

By contrast, the variant that adopts parallaxes from VB21 produces a 1σ shift in the absolute calibration α_0,RRab. While both VB21 and the baseline rely on Gaia parallaxes, they differ in the treatment of angular covariance and uncertainty propagation. As a result, the mean GC parallaxes used in VB21 are slightly offset from ours, with an average difference of 2 μ as across the anchor clusters.

This offset has a predictable effect: a change of 2 μ as in the input parallaxes translates into a ∼0.020 mag shift in the absolute-magnitude calibration. This matches the observed difference in α_0,RRab between the baseline and the VB21 variant. A similar shift appears in the variant that uses median parallaxes instead of the mean.

Fig. 5 Residuals of the baseline fit for RRab (bottom) and RRc (top). Blue histograms show host stars, orange histograms anchor stars. Overplotted are Gaussian distributions with identical mean and dispersion (blue for hosts, red for anchors). For anchor RRab, a Kolmogorov-Smirnov test reveals significant deviations from Gaussianity, with a sharper central peak and extended tails. RRc residuals are both more Gaussian-like and systematically tighter than those of RRab.

4.2.4 Models – exploring different assumptions

Removing RRab reduces the ${\chi }_{\text{dof}}^2$, while removing RRc increases it. This behavior arises from the use of a single extra-scatter term applied uniformly to all candles. As shown in Fig. 5, the residual spread is larger for RRab than for RRc. The shared extra-scatter value represents a compromise between the two (see Sect. 5.1), which explains the observed changes in the reduced chi-square when either type is removed. While the extra-scatter term could be fine-tuned for each candle type to achieve ${\chi }_{\text{dof}}^2\approx 1$ individually, it is kept fixed for simplicity. Overall, the model parameters demonstrate robustness, showing no significant changes when particular candles are excluded.

Model variants in which certain parameters are shared between standard candle types are also explored. Previous studies suggest that T2Cep may share their period slope with RRab (Braga et al. 2020), whereas the present results indicate that the RRc slope is somewhat closer to that of T2Cep. Nevertheless, all three slopes remain consistent within 2σ. To test these hypotheses, a variant enforcing β_T2Cep=β_RRab is implemented. As expected, the fit converges to an intermediate value between the independently determined β values from the baseline, while the absolute calibrations of both types remain unchanged. Similar behavior is observed in variants with shared β and α_M parameters, confirming that coupling these terms has minimal impact on the resulting absolute calibrations. The leverage on β is, however, limited for RRab and RRc due to their narrow period ranges (log P ≈[−0.4, 0] for RRab and log P ≈[−0.6,−0.4] for RRc).

Several variants are used to assess the impact of metallicity on the fit results. The most direct approach involves removing the metallicity dependence entirely by setting α_M=0 for all standard candles. Under this constraint, the other parameters for RRab and RRc remain nearly unchanged, as expected, given that their fitted metallicity terms on the intercept (α_M) are below 0.030 mag/dex. For T2Cep, however, α_M is more significant, and its removal results in a shift of β_T2Cep by 1.2σ. Despite these changes, the absolute calibrations of all candles remain stable, even when the metallicity dependence on the intercept is excluded.

Another set of variants introduces a metallicity dependence in the period slope by replacing β₀ with β₀+β_M[Fe/H] for both RRab and RRc, as has been suggested for classical Cepheid models (Anderson et al. 2016; Khan et al. 2025). This modification also allows probing the potential correlation observed in the MCMC results (see Fig. C.1) between β and α_M. For RRab, a significant metallicity dependence in the slope is detected, with β_0,2,RRab=−0.393 ± 0.078, corresponding to a 5σ detection. In contrast, no significant effect is found for RRc, with β_0,2,RRc=0.048 ± 0.110. In this variant, the metallicity effect on the intercepts (α_M) drops to zero, suggesting that the dependence initially attributed to α_M is effectively absorbed by the metallicity-dependent period slope term β_M. A separate variant where α_M is explicitly fixed to zero yields similar results. Repeating the analysis using GC metallicities from Harris (2010) shows that both β_M terms are consistent with zero within 1σ, indicating that the detection of metallicity dependence is sensitive to the adopted metallicity scale. This behavior most likely reflects limitations in the RRab photometric metallicity calibration of Li et al. (2023), underscoring the need for spectroscopic metallicities to reliably probe slope-metallicity dependencies. Nevertheless, in all of these variants, the absolute calibrations of the candles remain nearly unchanged.

4.2.5 Pre-processing

In the variants where the extra-scatter term is either suppressed, halved, or doubled, the expected changes appear in the overall χ², reflecting the adjusted uncertainties. These modifications, however, have only a marginal impact (<1σ) on the fitted parameters. Increasing the extra scatter primarily affects the uncertainties on the fitted parameters, particularly β and α_M, owing to the corresponding increase in χ². These results support the choice of adopting a single value of σ_extra=0.045 mag.

4.2.6 Anchors criteria

Changing the anchor selection criteria alters the distribution of clusters classified as anchors versus hosts. Stricter thresholds may reclassify marginal cases, with clusters that previously met the anchor requirements instead treated as hosts, and vice versa, potentially affecting the global fit.

The extinction-based criterion has limited impact on the overall results. Removing it produces only minor changes to the anchor/host classification, and tightening the threshold to E(B–V)=0.5 exerts negligible influence on the fitted parameters, despite a modest redistribution of clusters.

By contrast, the criterion based on signal-to-noise ratio (S/N) has a more significant effect. Doubling the S/N threshold from 10 to 20 induces a shift of about 1.5σ in the absolute calibration, along with slightly increased parameter uncertainties. While these changes remain statistically acceptable, the stricter threshold substantially reduces the number of anchor clusters, thereby increasing the sensitivity of the fit to individual clusters. This reduced anchor sample can introduce systematics that would otherwise be averaged out in a broader distribution.

To balance sample size with calibration stability, the baseline adopts an intermediate threshold of S/N=10, which provides a compromise between selectivity and robustness. This value is also consistent with the typical limit above which Gaia parallaxes can be reliably inverted (Luri et al. 2018).

5 Discussion

This section examines the residuals from the baseline analysis (Sect. 5.1). The RRab calibration is compared to similar results from the literature (Sect. 5.2), followed by a comparison of the derived GC distances with published values (Sect. 5.3). Finally, the impact of metallicity and the possibility of metallicity inhomogeneities within clusters are discussed (Sect. 5.4).

5.1 Residuals

Period–luminosity relations exhibit intrinsic scatter due to genuine astrophysical effects, as discussed in Sect. 2.3.1, arising from both the width of the instability strip and the HB evolutionary effect. It is therefore customary to include an extra-scatter term when fitting PLR, to account for these intrinsic contributions as well as observational uncertainties.

Figure 5 shows the distributions of fit residuals for RRab and RRc in anchor and host GCs. Most distributions are approximately Gaussian, except for RRab in anchors, which display a sharply peaked core and heavy tails, as confirmed by a Kolmogorov-Smirnov test. This deviation challenges the Gaussian error assumption. However, the MCMC analysis indicates that this did not bias the fit parameters.

We quantified the intrinsic scatter to apply using the interpercentile range (16th-84th percentiles), which is more robust to non-Gaussianity than the standard deviation. For anchor RRab, this yielded ∼0.045 mag, marginally larger than the scatter of RRc in anchors (σ=0.042 mag). For simplicity, we adopted a very small additional scatter of 0.045 mag for all variables, as the number of T2Cep was insufficient for a meaningful estimate. This is consistent with the 0.03−0.05 mag scatter predicted by theoretical pulsation models that account for HB evolutionary effects and the finite width of the instability strip (e.g., Bono et al. 2001; Marconi et al. 2015).

Several physical and observational explanations for the extended tails of the RRab anchor residuals were considered. Depth effects within GCs are ruled out, since their physical sizes (a few parsecs) are negligible compared to their several kpc distances. The Blažko (1907) effect, which induces amplitude and phase modulations in RRL light curves, was ruled out based on measured changes in mean magnitude over time (Skarka et al. 2020). Crowding and blending, which can bias photometric measurements, would likely produce skewed residuals rather than the symmetric heavy-tailed structure observed, and would be more applicable in host than in anchor GCs. Correlations between pulsation amplitudes across the three GDR3 photometric bands were also examined but revealed no linkages with the shapes of the fit residuals. Thus far, none of the mechanisms explored adequately explain the observed deviations from Gaussianity.

Variants using different values of the extra-scatter term showed expected changes in overall χ², while resulting in LL parameters well within 1σ of the baseline. This robustness indicates that while the choice of scatter influences statistical measures of goodness-of-fit, it does not bias the core results.

5.2 Comparison with the literature

The baseline results can be compared with the RRL calibration reported by Garofalo et al. (2022), which was based on Gaia EDR3 parallaxes and Gaia DR2 photometry of field RRab: $$M_G^W=(-0.{88}_{-0.09}^{+0.08})+(-2.{49}_{-0.20}^{+0.21})logP+(0.{14}_{-0.03}^{+0.03})[\mathrm{Fe}/\mathrm{H}].$$(22) A key difference with our results lies in the adopted pivot points: the present calibration is anchored at log P=−0.25 and [Fe/H]= −1.5, whereas theirs is defined at log P=0 and solar metallicity. Shifting our baseline fit to the same pivot values (Eq. (7), Table E.1) yields: $$\begin{array}{rl}{M}_G^W=& (-{0.986}_{-0.034}^{+0.034})+(-{2.582}_{-0.047}^{+0.047})\log P\\ & +(-{0.028}_{-0.008}^{+0.008})[\mathrm{F}\mathrm{e}/\mathrm{H}].\end{array}$$(23)

Thus we found the two calibrations to be broadly consistent in both the absolute magnitude zero point (${\alpha }_{0,\text{RRab}}^{\text{'}}$; agreement within 1.2σ) and the period slope (β_RRab; agreement within 1σ). However, the uncertainties from the present analysis are substantially smaller.

The most significant difference lies in the metallicity term. The calibration reported by Garofalo et al. (2022) implies a strong positive dependence of α_{M, RRab}=0.14 ± 0.03 mag/dex, corresponding to a 4.5σ detection. In contrast, the baseline solution yields a negative dependence of α_{M, RRab}=−0.028 ± 0.008 mag/dex, significant at the 3.5σ level and differing from their value by 5.4σ. However, as discussed in Sect. 4.2.2, the present analysis indicates that metallicity plays only a minor role in the absolute calibration of the candles.

Several methodological differences may underlie this discrepancy. Firstly, the metallicities adopted in Garofalo et al. (2022) were taken from Muraveva et al. (2018), who, in turn, had adopted values from Dambis et al. (2013). This latter catalog combined metallicity estimates from a variety of heterogeneous sources and methods and homogenized them onto a common scale. By contrast, our baseline analysis uses internally consistent photometric metallicities calibrated based on GDR3 light-curve parameters and homogeneous spectroscopic observations (Li et al. 2023). Secondly, the calibration based on bright (G ≲ 13 mag) field RRab in Garofalo et al. (2022) was not corrected for residual parallax systematics, which have been investigated in detail but remain complex in several respects (Khan et al. 2023, and references therein). Lastly, adopting mean-centered photometric metallicities minimizes parameter correlations and ensures that α₀ corresponds most closely to the observed stellar sample.

A recent study based on infrared photometry of field and GC RRL reported significantly positive values of α_M=0.20 ± 0.02 mag/dex both in individual passbands and infrared Wesenheit magnitudes (Bhardwaj et al. 2023). Several methodological differences between the Bhardwaj et al. analysis and the present one preclude a direct comparison, notably the distinct photometric bands employed and the adopted pivot points. Conceptually, our cluster variant, which yields null detections of α_M, is most similar in its approach. We note that the heavy element abundance scale adopted by Harris (2010) matches that by Carretta et al. (2009), which was used in Bhardwaj et al. (2023). A comparison of [Fe/H] in the sample of 11 GCs therefore unsurprisingly yields consistent values with a mean offset and standard deviation of −0.03 dex and 0.05 dex respectively.

Comparing our empirical calibration with theoretical predictions for both RRL and T2Cep may provide useful constraints for the underlying models. Marconi et al. (2021) predicted the following LL in the Gaia Wesenheit magnitude (DR2 photometric system) based on theoretical pulsation models of RRab and RRc: $$\begin{array}{rl}{M}_{G,\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}^W=& (-0.936\pm 0.049)+(-2.296\pm 0.032)\log P\\ & +(0.124\pm 0.005)[\mathrm{F}\mathrm{e}/\mathrm{H}],\end{array}$$(24) $$\begin{array}{rl}{M}_{G,\mathrm{R}\mathrm{R}\mathrm{c}}^W=& (-1.344\pm 0.033)+(-2.440\pm 0.028)\log P\\ & +(0.112\pm 0.005)[\mathrm{F}\mathrm{e}/\mathrm{H}].\end{array}$$(25) Analogously, theoretical LL for T2Cep (in the GDR3 photometric system) were recently presented by Marconi et al. (2025): $$\begin{array}{rl}{M}_{G,\text{T2Cep}}^W=& (-0.9917\pm 0.0032)+(-2.4655\pm 0.0055)\log P\\ & +(0.0065\pm 0.0022)[\mathrm{F}\mathrm{e}/\mathrm{H}].\end{array}$$(26) The above theoretical predictions were reported for different pivot points compared to the present work, namely at log P₀=0 and [Fe/H]₀=0. The most straightforward comparisons therefore concern the LL slopes β₀ (assumed linear) and the interceptmetallicity effect α_m. The slopes β₀ of RRc and T2Cep stars agree with our empirical calibration to within 1.0σ and 0.6σ, respectively. However, the predicted slope of RRab differs substantially (by 5.0σ) from our result. This could be particularly interesting in the context of a possible impact of helium abundance on LL, cf. Sect. 5.4. Significant discrepancies are seen for the intercept-metallicity coefficients α_m, which the models predicted to be significantly positive for both RRab and RRc, and close to 0 for T2Cep. By contrast, our analysis yielded virtually no metallicity dependence for RRab and RRc, and ambiguous results for T2Cep, which were generally dominated by elements related to the calibration of RRab or RRc. Formally, our results thus differ from the predictions by 16.38σ, 9.18σ, and 2.36σ for RRab, RRc, and T2Cep, respectively. However, these seemingly large deviations do not include the systematic uncertainties inherent to the models (whose small stated uncertainties we interpreted to reflect only formal fitting uncertainties) and should therefore be interpreted with caution.

The comparison of the predicted and observed LL intercepts is somewhat complicated by the disagreement in the metallicity effects and their uncertainties. Translating the (very precisely stated) theoretical relations to the pivot points adopted in this work yields significant differences for the intercepts of RRab, RRc, and T2Cep stars at the level of 4.6σ, 7.0σ, and 2.1σ, respectively. However, the significance levels of these differences may be misleading due to (a) the missing systematic uncertainties of the model predictions and (b) the need to translate to a common pivot metallicity using the aforementioned discrepant results in terms of metallicity effects.

5.3 Comparison of GC distances with the literature

Baumgardt & Vasiliev (2021, BV21) presented a catalog containing GC parameters, including distances featuring small uncertainties. The distances of these GCs were determined by averaging literature distance estimates from a variety of sources, often including highly correlated information, such as multiple estimates based on the same RRL member stars (e.g., analyzed in different studies). For example, the BV21 distance of NGC 5139 was estimated based on 27 published values, 14 of which were derived from RRL distances that share underlying assumptions and datasets. For each GC in the BV21 sample, the details of this averaging should be carefully considered as the available literature information varies from case to case. Compared to the GC distances in BV21, our homogeneously determined distances yield a mean difference of Δ μ=0.15 mag among anchors and Δ μ=0.17 mag among hosts, with dispersions of 0.08 and 0.15 mag, respectively. Figure 6 illustrates this comparison.

Based on the GC distances in BV21, Vasiliev & Baumgardt (2021) raised concerns about a possible distance-dependent bias of Gaia parallaxes in the low-S/N regime. However, our joint fit to the 93 GCs in our sample does not confirm such a trend. Instead, our mini-distance ladder in Fig. 3 shows no trend when comparing measured GC parallaxes to parallaxes obtained by inverting the precise luminosity distances by our fit. Thus, we find GDR3 GC parallaxes to remain reliable even at low S/N.

The RRL PLR calibrations presented by Bhardwaj et al. (2023), based on near-infrared photometry and anchored to field RRL, yield GC distances consistent with BV21. However, when applied to the LMC, these calibrations predict distances that are 0.06–0.12 mag larger than the detached eclipsing binary (DEB) distance of Pietrzyński et al. (2019). A similar offset was reported by Sicignano et al. (2024), who found Δ μ ≈ 0.09−0.15 mag between GC distances derived from LMC-calibrated T2Cep and those of BV21. Taken together, these results indicate that the DEB distance to the LMC is in tension with the BV21 GC distances inferred from both RRL and T2Cep.

The methodology adopted by Bhardwaj et al. (2017) differs from ours in several key respects. Their analysis relied on distances from the Bailer-Jones et al. (2021) catalog, which employs a Bayesian framework combining GDR3 parallaxes, Galactic priors, and bias corrections following Lindegren et al. (2021). Moreover, their sample is dominated by bright stars: 76% of RRab and 87% of RRc have G<13 mag, a regime where residual parallax offsets are known to be significant (see Khan et al. 2023, and references therein). For a typical GC anchor at 5 kpc (ϖ=200 μ as), a systematic parallax shift of +10 μ as would translate into a distance modulus shift of ≈−0.10 mag, thereby reconciling the GC distances with the LMC DEB distance. This suggests that the discrepancy between our results and those of Bhardwaj et al. (2023) is plausibly driven by residual parallax bias.

In summary, the systematic differences between our homogeneously determined GC distances, anchored to the parallaxes of Cruz Reyes et al. (2024), and those reported in the recent literature are best explained by (i) the averaging of heterogeneous and correlated values in the case of BV21, and (ii) uncorrected residual systematics in GDR3 parallaxes in the case of Bhardwaj et al. (2023).

Fig. 6 Comparison of distance moduli derived in this work for host (blue) and anchor (green) with literature values from Baumgardt & Vasiliev (2021). The top panel shows the direct comparison, while the bottom panel presents the residuals.

5.4 Metallicity effect and possible GC inhomogeneity

Our analysis reveals some limitations with respect to RRL metallicities. However, these limitations did not significantly affect the absolute calibrations of the standard candle intercepts.

Firstly, Fig. D.1 shows a trend between α_0,RRab, α_M,RRab, and the outlier-rejection threshold κ (Sect. 4.2.1). While we have adopted Chauvenet’s as an objective criterion for outlier rejection, we cannot distinguish whether inaccurate photometric metallicities or chemical inhomogeneities create the excess residual scatter. The former possibility is corroborated by alternative models (Sect. 4.2.4) involving a slope-metallicity dependence for RRab stars, which results in a putatively significant (5σ) detection of β_M,RRAb ≠ 0. Although a similar effect is predicted for classical Cepheids based on stellar evolution models (Khan et al. 2025), we stress that we consider the β_M, R RAb ≠ 0 result spurious, as it appears only for RRab stars with photometric metallicities and disappears when GC metallicities are considered instead. A period-dependent systematic in the calibration of photometric RRL metallicities based on P, R₂₁, and φ₃₁ by Li et al. (2023) could potentially explain this effect and appears to be a reasonable possibility given the inherent correlation between light curve shape and pulsation period. However, it is also possible that chemical inhomogeneities within GCs lead to excess LL residuals. In favor of this possibility is the observation that many GCs exhibit a larger spread in photometric metallicities than the reported scatter of their calibration. Specifically, Li et al. (2023) reported typical scatters of 0.24 dex (RRab) and 0.16 dex (RRc) for their photometric [Fe/H] calibration. In contrast, of the 20 GCs hosting >10 RRab, 5 exhibit σ_[Fe/H] > 0.3 dex; among the 9 GCs with >10 RRc, 6 show σ_[Fe/H]> 0.3 dex. Extreme cases include NGC 6715 (M 54), which exhibits a 0.44 dex spread in photometric [Fe/H] across 43 RRab, and NGC 6266 (M 62) with a dispersion of 0.42 dex across 24 RRc. These dispersions exceed the expected calibration scatter, suggesting that at least some GCs may host genuine internal metallicity spreads, as independently established for both M54 (Carretta et al. 2010, based on spectroscopic abundances) and M62 (Milone 2015, where a helium spread of Δ Y=0.08 was inferred to explain horizontal branch morphology). Evidence for chemical inhomogeneity among RRL in GCs was presented by Cruz Reyes et al. (2024) who showed that different helium abundances (Δ Y ≈ 0.06) were required to reproduce the instability-strip edges of RRab and RRc stars. This creates an exciting link to the study of multiple populations in GCs (Bastian & Lardo 2018). Helium and light-element variations would shift stars across the instability strip (e.g., Marconi et al. 2015, 2018), altering pulsation amplitudes and effective temperatures at fixed P.

These effects would naturally increase LL scatter when abundances are not determined individually (e.g., when GC-wide values for [Fe/H] are adopted) or if light-curve shapes depend on elements not considered in the spectroscopy-based calibration, such as helium. They may also contribute to the difference in LL slopes between our results and theoretical models (Sect. 5.2).

The two above paragraphs thus lead to the interesting possibility that inhomogeneities in helium abundance might introduce a period-dependent bias among the photometric iron abundances reported by Li et al. (2023). At present, this evidence remains circumstantial. However, detailed spectroscopic investigations of RRL stars in GCs would allow to directly test for intra-cluster inhomogeneity.

6 Summary and conclusions

This paper presents the first joint absolute calibration of the Leavitt Laws of RRab, RRc, and T2Cep stars. Our approach combines a linear least-squares framework, inspired by the SH0ES distance-ladder methodology, with an MCMC validation pipeline to assess uncertainties and robustness. The joint calibration of RRab, RRc, and T2Cep bring several advantages, notably through the ability of including a large number of GCs.

Our baseline analysis simultaneously calibrates LL of RRab, RRc, and T2Cep stars using reddening-free GDR3 Wesenheit magnitudes. The fiducial intercepts of RRab and RRc stars were thus calibrated to a distance accuracy of better than 1.0%, and T2Cep to 1.3%. We demonstrated robustness of this calibration using 31 fit variants exploring various systematics, including different parallax selections and metallicity treatments. A key result of our analysis are accurate and homogeneous distances to 93 GCs. These distances demonstrate that GC parallaxes based on GDR3 are not biased at the low-S/N end (i.e., as a function of distance). This is in contrast to previous claims in the literature, which we discuss in detail.

Our results were less clear in terms of the metallicity dependencies of RRL and T2Cep. For RRL stars, we used photometric metallicities from Li et al. (2023). While it is clear that the effect is very small, it is barely significant (3.5σ) for RRab stars and not significant (1.3σ) for RRc stars. For T2Cep, the metallicity effect was considered based on spectroscopic [Fe/H] for each cluster, resulting in a 2σ signal. As discussed in Sect. 5.4, these issues do not substantially affect the calibration of our LL intercepts. However, we point out that further work is required to understand the metallicity dependencies of population-II pulsators, in particular in GCs where we find tentative signs of chemical inhomogeneities among RRL in individual clusters.

Our work shows the high potential of population-II pulsators for high accuracy distance determination. Future work will incorporate parallaxes of field RRL stars and measure the LMC’s distance in order to strengthen the extragalactic distance scale by anchoring it to a large number of additional high-quality trigonometric parallaxes from Gaia.

Data availability

The results presented in Tables E.2 and E.3 are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/706/A87.

Acknowledgements

This research has received support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement No. 947660). R.I.A. is funded by a Swiss National Science Foundation Eccellenza Professorial Fellowship (award PCEFP2_194638). The authors thank Z. Prudil and M. Skarka for valuable discussions regarding the Blazhko effect. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://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 research has made use of NASA’s Astrophysics Data System; the SIMBAD database and the VizieR catalog access tool provided by CDS, Strasbourg. This work made use of the following open-source software: Python, NumPy (Harris et al. 2020), Pandas (pandas development team 2020), Matplotlib (Hunter 2007), SciPy (Virtanen et al. 2020), Astropy (Astropy Collaboration 2013, 2018, 2022), and emcee (Foreman-Mackey et al. 2013). We are grateful to the developers and maintainers of these packages for making their work available to the community.

References

Appendix A Joint fit

Equation (13) can be expressed in block form as this: $$\left[\begin{array}{c}{\mathit{m}}_{\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}^h\\ {({\mathit{m}}^a-{\mathit{\mu }}_0)}_{\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}\\ {\mathit{m}}_{\mathrm{R}\mathrm{R}\mathrm{c}}^h\\ {({\mathit{m}}^a-{\mathit{\mu }}_0)}_{\mathrm{R}\mathrm{R}\mathrm{c}}\\ {\mathit{m}}_{\mathrm{T}2\mathrm{C}}^h\\ {({\mathit{m}}^a-{\mathit{\mu }}_0)}_{\mathrm{T}2\mathrm{C}}\\ \mathbf{0}\end{array}\right]=\left[\begin{array}{ccccc}{\mathbb{I}}_{\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}^h& 0& {\mathbb{P}}_{\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}^h& 0& 0\\ 0& {\mathbb{I}}_{\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}^a& {\mathbb{P}}_{\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}& 0& 0\\ {\mathbb{I}}_{\mathrm{R}\mathrm{R}\mathrm{c}}^h& 0& {\mathbb{B}}_{\mathrm{R}\mathrm{R}\mathrm{c}}^h& {\mathbb{P}}_{\mathrm{R}\mathrm{R}\mathrm{c}}^h& 0\\ 0& {\mathbb{I}}_{\mathrm{R}\mathrm{R}\mathrm{c}}^a& {\mathbb{B}}_{\mathrm{R}\mathrm{R}\mathrm{c}}^a& {\mathbb{P}}_{\mathrm{R}\mathrm{R}\mathrm{c}}^a& 0\\ {\mathbb{I}}_{\mathrm{T}2\mathrm{C}}^h& 0& {\mathbb{B}}_{\mathrm{T}2\mathrm{C}}^h& 0& {\mathbb{P}}_{\mathrm{T}2\mathrm{C}}^h\\ 0& {\mathbb{I}}_{\mathrm{T}2\mathrm{C}}^a& {\mathbb{B}}_{\mathrm{T}2\mathrm{C}}^a& 0& {\mathbb{P}}_{\mathrm{T}2\mathrm{C}}^a\\ 0& \mathrm{I}\mathrm{d}& 0& 0& 0\end{array}\right]\cdot \left[\begin{array}{c}\mathit{\mu }\\ \mathit{\mu }\mathit{\mu }\\ {\mathit{p}}_{\mathrm{R}\mathrm{R}\mathrm{a}\mathrm{b}}\\ {\mathit{p}}_{\mathrm{R}\mathrm{R}\mathrm{c}}\\ {\mathit{p}}_{\mathrm{T}2\mathrm{C}}\end{array}\right],$$(A.1) where bold symbols denote vectors (e.g., m^h, μ₀), and blackboard bold letters represent matrices (e.g., 𝕀^h, 𝕇^h, Id). The superscripts h and a refer to hosts and anchors, respectively, while subscripts distinguish the different types of variable stars.

In this representation, the model structure is straightforward to visualize. The first row of the matrix corresponds to the hosts for RRab as given by Eq. (7), the second row to the anchors for RRab as described in Eq. (9). Rows 3 and 4 follow the same logic for RRc (Eqs. (8) and (10)), rows 5 and 6 for T2Cep (Eqs. (8) and (10)), and the last row represents the constraints on the known geometric distances of the anchors as defined in Eq. (11).

Appendix A.1 Observation vectory

The vector y consists of two components, m^h and (m^a−μ₀), repeated for each type of variable star, and a null vector 0. The component m^h contains the magnitudes of stars in all hosts and is further subdivided into sub-vectors m_i corresponding to each cluster i. Each sub-vector m_i contains the magnitudes m_{i, j} of the stars j within cluster i: $${\mathit{m}}^h={\left[\begin{array}{c}{\mathit{m}}_1^W\\ ⋮\\ {\mathit{m}}_n^W\end{array}\right]}_{({\sum }_i^n{n}_i\times 1)}\text{with}{\mathit{m}}_i^W={\left[\begin{array}{c}{m}_{i,1}^W\\ ⋮\\ m_{i,n_i}^W\end{array}\right]}_{(n_i\times 1)},$$(A.2) where n is the total number of host GC and n_i is the number of stars of a given type (RRab, RRc, or T2Cep) in host i.

The component (m^a−μ₀) represents the magnitudes of stars in all anchors, corrected by the known distance modulus μ_{0, i}. Like the previous component, this vector is composed of sub-vectors m_i−μ_{0, i} for each anchor i, containing the corrected magnitudes m_{i, j}−μ_{0, i} for each star j in anchor cluster i: $$({\mathit{m}}^a-\mathit{\mu })={\left[\begin{array}{c}{\mathit{m}}_1^W-{\mathit{\mu }}_{0,1}\\ ⋮\\ {\mathit{m}}_{n^{\text{'}}}^W-{\mathit{\mu }}_{0,n^{\text{'}}}\end{array}\right]}_{({\sum }_i^{n^{\text{'}}}{n}_i^{\text{'}}\times 1)}\text{with}{\mathit{m}}_i^W={\left[\begin{array}{c}{m}_{i,1}^W-{\mu }_{0,i}\\ ⋮\\ m_{i,n_i^{\text{'}}}^W-{\mu }_{0,i}\end{array}\right]}_{(n_i^{\text{'}}\times 1)},$$(A.3) where n^′ is the total number of anchor GC and $n_i^{\text{'}}$ the number of stars of a given type in anchor i. The last component of y is the null vector 0, with one entry for each anchor.

Appendix A.2 Design matrix L

The matrix 𝕃 is composed of two main components, 𝕀 and 𝕇, repeated for each type of variable star. The matrix 𝕀^h is similar to an identity matrix but with its diagonal filled with vectors 1 (containing only ones), while the remaining elements are zeros: $${\mathbb{I}}^h={\left[\begin{array}{ccc}{\mathbf{1}}_{n_1}& \dots & {\mathbf{0}}_{n_1}\\ ⋮& \ddots & ⋮\\ {\mathbf{0}}_{n_n}& \dots & {\mathbf{1}}_{n_n}\end{array}\right]}_{(\sum _i^n{n}_i\times n)}\text{\hspace{0.17em}with\hspace{0.17em}}{\mathbf{1}}_{n_i}={\left[\begin{array}{c}1\\ ⋮\\ 1\end{array}\right]}_{(n_i\times 1)}\text{\hspace{0.17em}and\hspace{0.17em}}\mathbf{0}={\left[\begin{array}{c}0\\ ⋮\\ 0\end{array}\right]}_{(n_i\times 1)}.$$(A.4)

The matrix 𝕇 encodes information about properties (period and metallicity) for each host (or anchor) cluster. It consists of sub-vectors: 1 filled with ones, P_i containing the logarithms of periods P_i,j relative to the pivot P₀, and [Fe/H]_i containing the metallicities relative to the pivot [Fe/H]₀: $${\mathbb{P}}^h={\left[\begin{array}{ccc}{\mathbf{1}}_{n_1}& {\mathit{P}}_1& {[\mathrm{Fe}/\mathrm{H}]}_1\\ ⋮& ⋮& ⋮\\ {\mathbf{1}}_{n_n}& {\mathit{P}}_n& {[\mathrm{Fe}/\mathrm{H}]}_n\end{array}\right]}_{({\sum }_i^n{n}_i\times 3)},$$ with $${\mathit{P}}_i={\left[\begin{array}{c}log\left(P_{i,1}\right)-log\left(P_0\right)\\ ⋮\\ log\left(P_{i,n_i}\right)-log\left(P_0\right)\end{array}\right]}_{(n_i\times 1)},$$ and $${[\mathrm{Fe}/\mathrm{H}]}_i={\left[\begin{array}{c}{[\mathrm{Fe}/\mathrm{H}]}_{i,1}-{[\mathrm{Fe}/\mathrm{H}]}_0\\ ⋮\\ {[\mathrm{Fe}/\mathrm{H}]}_{i,n_i}-{[\mathrm{Fe}/\mathrm{H}]}_0\end{array}\right]}_{(n_i\times 1)}.$$(A.5)

The matrix 𝔹^h represents the absolute calibration of the primary candles (α_0,RRab) when propagated to other types of variables, and adopts a similar block structure: $${\mathbb{B}}^h={\left[\begin{array}{ccc}{\mathbf{1}}_{n_1}& 0& 0\\ ⋮& ⋮& ⋮\\ {\mathbf{1}}_{n_n}& 0& 0\end{array}\right]}_{({\sum }_i^n{n}_i\times 3)}$$(A.6)

For all anchor components (𝕀^a, 𝕇^a, 𝔹^a), the number of hosts n and stars per host n_i are replaced by the number of anchors n^′ and stars per anchor $n_i^{\text{'}}$. The final component, Id, is an identity matrix of size (n^′ × n^′), representing the geometric constraint.

Appendix A.3 Parameter vector q

The parameter vector q consists of three components: μ, Δ μ, and p, with the last one repeated for each type of variable star. The component μ contains the fit parameters for the distances μ_i to the hosts, while Δ μ includes the correction parameters Δ μ_i for the anchors. Finally, p contains the parameters describing the properties of a given type of variable star: the absolute magnitude α₀ for a fiducial star at the pivot period and metallicity, the slope β of the LL with respect to period, and the slope α_M of the PLR with respect to metallicity: $$\mathit{\mu }={\left[\begin{array}{c}{\mu }_1\\ ⋮\\ {\mu }_n\end{array}\right]}_{(n\times 1)},\mathit{\Delta }\mathit{\mu }={\left[\begin{array}{c}\Delta {\mu }_{1^{\text{'}}}\\ ⋮\\ \Delta {\mu }_{n^{\text{'}}}\end{array}\right]}_{(n^{\text{'}}\times 1)},\mathit{p}={\left[\begin{array}{c}{\alpha }_0\text{or}\Delta {\alpha }_0\\ {\beta }_0\\ {\alpha }_M\end{array}\right]}_{(3\times 1)}.$$(A.7)

Appendix A.4 Covariance Matrix C

The covariance matrix C captures the uncertainties associated with the measurements in y and is approximated as diagonal. For the host and anchor stars, the uncertainties in m and (m−μ) are given by the photometric uncertainties $\sigma \left(m_{i,j}^W\right)$. For anchors, the geometric constraints are expressed as 0 in y, with their uncertainties represented by σ(μ_{0, i}) in the corresponding diagonal entries of C: $$\mathit{C}=diag\left[\begin{array}{c}\sigma \left({\mathit{m}}_{\text{RRab}}^h\right)\\ \sigma \left({\mathit{m}}_{\text{RRab}}^a\right)\\ \sigma \left({\mathit{m}}_{\text{RRc}}^h\right)\\ \sigma \left({\mathit{m}}_{\text{RRc}}^a\right)\\ \sigma \left({\mathit{m}}_{\text{T2Cep}}^h\right)\\ \sigma \left({\mathit{m}}_{\text{T2Cep}}^a\right)\\ \sigma \left({\mathit{\mu }}_0\right)\end{array}\right].$$(A.8)

This diagonal approximation neglects possible correlations between stars but greatly simplifies the fitting process. In Sect. 4.1.1, we assess the validity of this assumption by comparing the baseline fit with an MCMC-based implementation.

Appendix A.5 Least-Squares Solution and Parameter Estimation

The system of equations in Eq. (13) is solved with a least-squares minimization. The goodness of fit is quantified with the χ² statistic: $${\chi }^2={(\mathit{y}-\mathit{L}\mathit{q})}^T{\mathit{C}}^{-1}(\mathit{y}-\mathit{L}\mathit{q}).$$(A.9)

Minimization of χ² yields the maximum-likelihood parameters: $${\mathit{q}}_{\text{best}}={\left({\mathit{L}}^T{\mathit{C}}^{-1}\mathit{L}\right)}^{-1}{\mathit{L}}^T{\mathit{C}}^{-1}\mathit{y}.$$(A.10)

The associated covariance matrix, which encodes the uncertainties of the fitted parameters, is $${\mathit{\Sigma }}_{\text{best}}={\left({\mathit{L}}^T{\mathit{C}}^{-1}\mathit{L}\right)}^{-1}$$(A.11)

Appendix B Outlier rejection algorithm

This appendix presents the detailed formulation of the outlier rejection algorithm used in Sect. 3.4. The procedure is summarized in Algorithm 1, which generalizes the classical σ-clipping approach to account for the population-dependent scatter of different types of pulsating stars.

Appendix C MCMC corner plot

For completeness, we provide in Fig. C.1 the corner plot of the MCMC sampling described in Sect. 3.3, with results discussed in Sect. 4.1.1. The figure illustrates the posterior distributions of selected parameters and their correlations.

Appendix D Variant summary plots

Figure D.1 provides a visual summary of the absolute calibration parameters (α₀, Δ α₀) and metallicity terms (α_M) for the standard candles across all model variants (excluding the ’No RRab’ case, where the absolute scale and reference candle change).

Appendix E Result tables

This appendix compiles the tables referenced throughout the analysis. Table E.1 summarizes the numerical results from all fits, including the baseline solution and the variants described in Sect. 3.5. Reported are the best-fit parameters for the different standard candles, together with the number of stars included in each fit and the corresponding reduced chi-squared values.

Tables E.2 and E.3 list the GC classified as anchors and hosts, respectively. Each table provides the number of stars used, the distances derived from the baseline fit, the mean parallaxes of the clusters, and the literature distances from Baumgardt & Vasiliev (2021).

Fig. C.1 MCMC sampling of the χ2 statistic for the global fit, showing posterior distributions of selected parameters. Shown are the candle parameters α0, β, and αM for RRab and RRc stars. The corresponding parameters for T2Cep exhibit similar distributions and are omitted here for clarity. The plot also includes one representative host distance and one anchor correction factor. Posterior distributions (black) are well approximated by Gaussians, with means and dispersions closely matching the least-squares solution from the baseline (red). Contours indicate the 1σ, 2σ, and 3σ confidence regions. Correlations between period slopes and metallicity terms are discussed further in Sect. 4.2.2.

Fig. D.1 Visual summary of the absolute calibration parameters (α0, Δ α0) and metallicity terms (αM) for the standard candles across all model variants (excluding the ’No RRab’ case, where the absolute scale and reference candle change). The red dashed line and shaded band indicate the baseline values and their associated uncertainties, white the blue markers and error bars represent the values obtained from each variant.

All Tables

All Figures

Fig. 1 Mean residual parallax difference as a function of G-band magnitude, computed with rolling bins of width 0.5 mag and step size 0.05 mag. The x-axis marks the bin midpoints, while the y-axis gives the corresponding mean residual. Points are color-coded by the number of stars per bin, and the shaded area shows the 1σ uncertainty of the mean. The inset zooms into the range populated by most RRL and T2Cep stars, where a typical residual trend of ∼3−4 μas is seen.

In the text

Fig. 2 Baseline joint fit to fundamental-mode, first-overtone RR Lyrae, and type-II Cepheids in Galactic GC. Colors indicate different stellar types: blue for RRab, green for RRc, yellow for T2Cep, and red for stars excluded by the outlier rejection. Circles mark stars in hosts, while downward-pointing triangles denote stars in anchors. Top: absolute magnitudes, corrected for metallicity, plotted against log P. For hosts, $M_{i,j}^W=m_{i,j}^W-{\mu }_{0,i}-{\alpha }_M({[\mathrm{Fe}/\mathrm{H}]}_{i,j}-{[\mathrm{Fe}/\mathrm{H}]}_0)$; for anchors, an additional −Δ μ0, i term applies. Colored dashed lines show the fitted PLR. Bottom: residuals between the corrected magnitudes and the fitted PLR, with dashed lines marking the residual standard deviation. The inset (bottom right) shows the residual distributions.

In the text

Fig. 3 Visualization of the baseline mini-distance ladder. The lower-left panel shows anchor GC parallaxes on the x-axis against the joint-fit distances on the y-axis, including residuals. The mean difference between Gaia parallax and fit parallax for anchor GCs is −0.87 ± 1.52 μ as. The upper-right panel shows the host GC fit distances on the x-axis against the host GC parallaxes on the y-axis, with a mean difference of −5.9 ± 4.3 μ as affected by outliers. The two outlier hosts Terzan 1 and BH 229 are due to sparse membership and do not affect the fit. The comparison reveals no evidence for a distance dependent bias of GDR3 parallaxes.

In the text

Fig. 4 Correction parameters Δ μ for each anchor GC in the baseline. Blue symbols show the initial values (set to zero; Eq. (11)), and red symbols the fitted values (Eq. (A.10)). Clusters are ordered by increasing distance, so the typical uncertainty in Δ μ grows along the x-axis, although the number of member stars also contributes. Overall, 70% of anchors remain within 1σ of their Gaia parallax values, and all lie within 2σ.

In the text

Fig. 5 Residuals of the baseline fit for RRab (bottom) and RRc (top). Blue histograms show host stars, orange histograms anchor stars. Overplotted are Gaussian distributions with identical mean and dispersion (blue for hosts, red for anchors). For anchor RRab, a Kolmogorov-Smirnov test reveals significant deviations from Gaussianity, with a sharper central peak and extended tails. RRc residuals are both more Gaussian-like and systematically tighter than those of RRab.

In the text

	Fig. 6 Comparison of distance moduli derived in this work for host (blue) and anchor (green) with literature values from Baumgardt & Vasiliev (2021). The top panel shows the direct comparison, while the bottom panel presents the residuals.
In the text

Fig. C.1 MCMC sampling of the χ2 statistic for the global fit, showing posterior distributions of selected parameters. Shown are the candle parameters α0, β, and αM for RRab and RRc stars. The corresponding parameters for T2Cep exhibit similar distributions and are omitted here for clarity. The plot also includes one representative host distance and one anchor correction factor. Posterior distributions (black) are well approximated by Gaussians, with means and dispersions closely matching the least-squares solution from the baseline (red). Contours indicate the 1σ, 2σ, and 3σ confidence regions. Correlations between period slopes and metallicity terms are discussed further in Sect. 4.2.2.

In the text

Fig. D.1 Visual summary of the absolute calibration parameters (α0, Δ α0) and metallicity terms (αM) for the standard candles across all model variants (excluding the ’No RRab’ case, where the absolute scale and reference candle change). The red dashed line and shaded band indicate the baseline values and their associated uncertainties, white the blue markers and error bars represent the values obtained from each variant.

In the text