© The Authors 2026

1. Introduction

Classical Cepheids (DCEPs) are crucial standard candles in the extragalactic distance scale due to the Leavitt law (Leavitt & Pickering 1912), which defines a correlation between their period and luminosity – a PL relation. Once calibrated using independent distances derived from geometric methods such as trigonometric parallaxes, eclipsing binaries, and water masers, these relations form the first step in constructing the cosmic distance scale. The first rung provides the basis for calibrating secondary distance indicators, such as Type Ia supernovae (SNIa), which enable us to measure distances to distant galaxies within the steady Hubble flow. The calibration of this three-step process (commonly referred to as the cosmic distance ladder) allows us to reach the Hubble flow, where the Hubble constant (H₀) – the value linking the distance to the recession velocity of galaxies – can be estimated (e.g. Sandage & Tammann 2006; Freedman et al. 2012; Riess et al. 2016, and references therein).

The value of H₀ is a key parameter in cosmology as it determines both the scale and age of the Universe. Consequently, achieving a 1% precision in its measurement is one of the primary goals of modern astrophysics. However, there is an ongoing and well-known discrepancy between the H₀ values obtained by the SH0ES¹ project using the cosmic distance ladder (H₀ = 73.01 ± 0.99 km s⁻¹ Mpc⁻¹, Riess et al. 2022b) and those measured by the Planck mission from the cosmic microwave background under the flat Λ cold dark matter (ΛCDM) model (H₀ = 67.4 ± 0.5 km s⁻¹ Mpc⁻¹, Planck Collaboration VI 2020). The tension seems to be confirmed by new James Webb Space Telescope observations by the SH0ES team (Riess et al. 2024), albeit Freedman et al. (2024) results, again based on JWST data but using a different sample of SNIa, appear to favour a reduction in the size of the discrepancy between the early and late estimates of H₀. On the other hand, if confirmed, the 5σ discrepancy would signal the need to revise the ΛCDM model.

One of the remaining sources of uncertainty in the cosmic distance ladder is the debated metallicity dependence of the DCEP PL relations used to calibrate secondary distance indicators. A variation in metallicity is expected to influence the shape and width of the DCEP instability strip (e.g. Caputo et al. 2000), which in turn alters the coefficients of the PL relations (Marconi et al. 2005, 2010; De Somma et al. 2022, and references therein). Although it has been demonstrated that the metallicity dependence of the PL relations alone cannot solve the Hubble tension (e.g. Riess et al. 2022b), its precise value could nevertheless be crucial to establish the actual size of the discrepancy in σs between the cosmic ladder and the CMB+ΛCDM measurements of H₀. Furthermore, accurately knowing the metallicity dependence of PL relations allows us to constrain pulsation models and validate their physical assumptions.

Direct empirical tests of metallicity effects on the PL and PW relations, based on Galactic DCEPs with reliable [Fe/H] abundances from high-resolution (HiRes) spectroscopy, have long been limited by the lack of precise independent distances for enough Milky Way (MW) DCEPs. The Gaia mission (Gaia Collaboration 2016) has radically changed this. Starting with data release 2 (DR2) (Gaia Collaboration 2018) and further improved in early data release 3 (EDR3) (Gaia Collaboration 2021), Gaia has provided accurate parallaxes. It has also led to the discovery of hundreds of new Galactic DCEPs (Clementini et al. 2019; Ripepi et al. 2019, 2023), alongside other surveys such as OGLE (Udalski et al. 2018) and ZTF (Chen et al. 2020), creating a large, valuable DCEP sample for both extragalactic distance scale and Galactic studies (e.g. Lemasle et al. 2022; Trentin et al. 2023, 2024b).

Yet until recently, DCEPs with HiRes-based metallicities were mostly confined to the solar neighbourhood, covering a narrow [Fe/H] range, around the solar value, with small scatter (0.2–0.3 dex) (e.g. Genovali et al. 2014; Luck 2018; Groenewegen 2018; Ripepi et al. 2019). The limited range makes it difficult to assess metallicity effects on Galactic Cepheid PL relations with strong statistical significance.

To tackle this issue, a few years ago we launched the C-MetaLL project² (Cepheid – Metallicity in the Leavitt Law; see Ripepi et al. 2021a, for details). Its goal is to measure chemical abundances for approximately 400 Galactic DCEPs using HiRes spectroscopy, and specifically extending the iron abundance range into the metal-poor regime, particularly [Fe/H] < −0.4 dex. In the project’s papers I, II, V, and VI (Ripepi et al. 2021a; Trentin et al. 2023; Bhardwaj et al. 2024; Trentin et al. 2024b), we provided accurate abundances of over 25 chemical species for 290 DCEPs spanning a wide range in metallicity (+0.3 < [Fe/H] < −1.1 dex; see Sect. 2.3) and with an approximately uniform number of DCEPs at all [Fe/H] values.

In Trentin et al. (2024b) we exploited this sample to study the chemical composition of the Galactic disc and spiral arms. In this paper, we investigate the metallicity dependence of the DCEPs’ PW relations using our updated and homogeneous sample alone, i.e. without using literature abundances as in our previous papers, due to insufficient proprietary data. In our previous works (Ripepi et al. 2020, 2021a, 2022a; Bhardwaj et al. 2024; Trentin et al. 2024a), we have found a generally large metallicity dependence of the DCEPs PL and PW relations, of the order of γ ∼ −0.3/−0.5 mag/dex, where γ is the usual notation for the metallicity term of the intercept in the PL and PW relations (see Sect. 3 for details). These findings are in contrast to those of the SH0ES team, who identified a smaller dependence (∼ − 0.2 mag/dex Riess et al. 2021) based on a local sample of 75 DCEPs with a narrow [Fe/H] range. Moreover, our results provide values of γ that are larger (in the absolute sense) than those reported in the past by several authors (e.g. Scowcroft et al. 2009; Gieren et al. 2018; Cruz Reyes & Anderson 2023) and, in particular, by Breuval et al. (2021, 2022, 2024), who, similarly to what was reported in SH0ES works, found a mild metallicity dependence from the analysis of DCEP samples in the MW and the Large and Small Magellanic Clouds (LMC and SMC). These systems show different metallicities, decreasing from about solar to [Fe/H] ∼ –0.75 dex. Bhardwaj et al. (2023) showed that the two approaches of using spectroscopic metallicities of individual Galactic Cepheids with Gaia parallaxes and using the MW and Magellanic Clouds as representative samples provide different quantifications of the metallicity dependence of PL relations. Even more discrepant are the results from Madore & Freedman (2025), who performed five empirical tests across a broad metallicity range and various wavelength regimes (optical to mid-infrared), concluding that no statistically significant metallicity dependence could be detected in the available PL relations. In a subsequent paper, Madore et al. (2025) reinforced their conclusions using theoretical static stellar atmosphere models to evaluate the differential impact of metallicity on Cepheid spectral energy distributions. However, a recent and detailed re-examination of the methods used by Madore & Freedman (2025) has suggested potential uncertainties regarding the robustness of their conclusions (Breuval et al. 2025). Finally, in a recent paper, Wang & Chen (2025), adopting basically the same input data and a very similar approach to that of our previous paper C-MetaLL-IV (Trentin et al. 2024b), found very similar results to ours in the same filters.

Intending to shed some light on this intricate matter, in this paper we analyse the metallicity dependence of the PW relations for Galactic DCEPs. We do not address the PL relations because of the significant uncertainties on the interstellar extinction for a sizeable fraction of our targets, which are placed far away in the disc and can easily reach values of E(B − V) larger than one magnitude. The use of Wesenheit magnitudes substantially mitigates this problem, especially in the NIR bands, where the coefficients multiplying the colour are small.

The paper is structured as follows. In Sect. 2 we present the sample used in our analysis, in Sect. 3 we describe the analysis procedure, and in Sect. 4 we describe our results. We discuss them and give our conclusion in Sect. 6.

2. Dataset

In this section, we describe in detail the properties of the DCEP sample used in this work. To calculate the desired period–Wesenheit–metallicity (PWZ) relations, we require additional ingredients besides the abundances, such as the multi-band photometry and parallaxes. The collection of this information is described in the following sections. All data used in this work are provided in Table A.1.

2.1. Periods of pulsation and photometry

As in our previous works, the relevant data for the DCEP sample used in this work have mainly been taken from the list published in the Gaia Data Release 3 (DR3, see Ripepi et al. 2023). From the same source, we adopted periods, epochs of maximum light, amplitudes in the G band, and the intensity-averaged G, G_BP, and G_RP magnitudes, derived according to a specific treatment of Cepheid pulsating stars (see Clementini et al. 2016, 2019; Ripepi et al. 2019, 2023). For the 15 stars without specific treatment in Gaia DR3, we adopted periods from Pietrukowicz et al. (2021) and average magnitudes from the general Gaia DR3 catalogue³ (see Gaia Collaboration 2023).

The Gaia photometry was used to calculate homogeneous V,  I magnitudes for 248 stars using the precise transformations by Pancino et al. (2022), which were corrected for a small shift as is explained in Trentin et al. (2024a). For the remaining 42 (bright) stars, we adopted the photometry by Bhardwaj et al. (2015), which, as has been demonstrated by Bhardwaj et al. (2024), is consistent with the V,  I magnitudes obtained with the Pancino et al. (2022) transformations corrected as mentioned above.

Photometry in the NIR bands J, H, K_s was collected from the literature using intensity-averaged magnitudes from full light curves when possible. In particular, we used the data by (i) Bhardwaj et al. (2015) and Groenewegen (2018) for 36 and 6 bright DCEPs, respectively; (ii) Bhardwaj et al. (2024) data collected in the context of the C-MetaLL survey for an additional 61 stars (paper V of the series). For the remaining 187 stars, we had to rely on the single-epoch photometry from the 2MASS survey (Skrutskie et al. 2006), which was used together with template light curves to estimate average magnitudes in the NIR bands. To this aim, we proceeded exactly as is described in Ripepi et al. (2021a) and Trentin et al. (2024a). We recall that the periods and epochs of maximum light have been taken from Gaia DR3 when available and from the ASAS-SN survey (All-Sky Automated Survey for Supernovae Shappee et al. 2014; Christy et al. 2023) for seven stars for which epochs were missing in Gaia DR3.

To calculate the Wesenheit magnitudes in the HST bands, we proceeded as in Trentin et al. (2024a), i.e. we adopted the Riess et al. (2021) transformations between Johnson-Cousins H,  V , I and the F160W,  F555W,  F814W bands. In the following, to avoid confusion, we call cF160W,  cF555W,  cF814W the HST-like magnitudes obtained starting from ground-based data, where the ‘c’ stands for ‘calculated’. The same prefix will be applied to the Wesenheit magnitudes obtained with these data (see Table 1). To verify the reliability of the calculated HST magnitudes, we carried out a test which is described in detail in Appendix B. As a result, we find that apart from some scattering, the calculated HST magnitudes (and Wesenheit functions) do not present any relevant systematic difference with the native HST photometry.

The adopted photometry is listed in Table A.1, where different labels point to the sources of the distinct values. We note that, although the spectroscopy used in this paper is fully homogeneous, as is shown in this section, the photometry is not yet. The analysis of a fully homogeneous sub-sample of C-MetaLL data can be found in Bhardwaj et al. (2024).

2.2. Fundamentalisation of 1O pulsators

In a recent work, Pilecki (2024) discussed in detail the fundamentalisation of the periods of first-overtone (1O) pulsators. They provided new empirical relationships for different environments such as the MW, LMC, and SMC. However, they also specified that for the aim of using the fundamentalisation for the PL and PW relations, it is preferable to use the relations they had published in their previous work (Pilecki et al. 2024), in which they had provided new relations for LMC and SMC. Therefore, in this work, we decided to use for all the 1O (and also for first over second overtone modes – 1O2O) DCEPs the following relation from Pilecki et al. (2024):

$$\begin{array}{c}\hfill P_F=P_{1O}(1.418+0.114log{P}_{1O}),\end{array}$$(1)

which is valid for the LMC, i.e. for [Fe/H] ∼ –0.4 dex (with σ = 0.07 dex Romaniello et al. 2022). This is justified by the fact that: (i) this relation has been accurately derived using the best data available (the LMC DCEPs are almost all at the same distance, contrarily to those in the MW and SMC); and (ii) the average metallicity of the LMC is close to the mean metallicity of our sample of 1O pulsators, which is exactly ⟨[Fe/H]⟩_1O ∼ −0.4 dex, with a dispersion of ∼0.3 dex.

2.3. Abundances

The dataset utilised in this work comprises DCEPs whose abundances have been published in our previous papers. In more detail, we adopted data for one star from both Catanzaro et al. (2020) and Ripepi et al. (2021a), and 47, 65, 42, and 134 stars from Ripepi et al. (2021b), Trentin et al. (2023), Bhardwaj et al. (2023), and Trentin et al. (2024b), respectively⁴. This data comprises 290 DCEPs, for which we have provided abundances for the iron and more than 20 chemical species in a homogeneous way. The sample’s period and iron abundance ranges are shown in Fig. 1. To obtain a good fit of the PLZ and PWZ relations, the entire space of these two parameters should be filled evenly. Although a larger number of long-period metal-poor DCEPs would be desirable (we are collecting additional observations for this aim in the context of the C-MetaLL project), the sample appears already to be rather equilibrated. In particular, a comparison with the sample adopted in Trentin et al. (2024b), displayed in Fig. C.1, reveals a substantial improvement in this respect.

Fig. 1. Periods and iron abundances spanned by the investigated sample. The different pulsation modes are labelled in the figure.

2.4. Extinction law and definition of the Wesenheit magnitudes

In this work, for homogeneity with our previous papers, we adopted the standard Cardelli et al. (1989) extinction law with R_V = 3.1 as the baseline. In addition, we adopted the Fitzpatrick (1999) results, when the difference with the previous one was significant (e.g. in the V , I and J, K_S colours). Some specific Wesenheit magnitudes were taken directly from the literature. All the quantities used in this work are listed in Table 1.

2.5. Parallaxes

All the 290 objects in our sample have a valid parallax value from Gaia DR3. In more detail, three objects have negative parallaxes, while 141 and 237 objects have a relative error on the parallax of σϖ/ϖ< 0.1 and 0.2, respectively. Figure 2 shows the variation in σϖ/ϖ with the GaiaG magnitude. As expected, fainter objects have less precise parallaxes. The figure also shows that more metal-rich objects generally have more precise parallaxes, because, on average, they are closer to the Sun. The farthest objects, which are also the most metal-poor, have in general less precise parallaxes, which is why in the C-MetaLL project we aim to further enlarge the sample of metal-poor objects, i.e. to compensate for the lower precision of individual sources with statistics. A comparison with the sample adopted in Trentin et al. (2024b) (see Fig. C.2) confirms the above considerations and testifies to the improvement of the present work sample in this context.

Fig. 2. Relative parallax error as a function of the GaiaG magnitude. The points are colour-coded according to their iron abundance.

Concerning the goodness of the parallaxes, 19 objects have a Gaia RUWE parameter (Renormalised Unit Weight Error Lindegren et al. 2021) beyond the generally accepted good threshold of 1.4. On the other hand, we have 15 objects with astrometric_gof_al below the threshold of 12.5 adopted by Riess et al. (2021) for their DCEP sample. The distribution of these parameters for our sample is shown in Fig. 3⁵. The slight gap on the y axis between 12 and 14 seems to justify a preliminary cut at astrometric_gof_al = 12.5, and thus keeping four objects with RUWE slightly exceeding the value 1.4. In any case, our fitting procedure allows for outlier rejection; hence, the occurrence of a few DCEPs with possibly incorrect parallaxes will not impact the results. The figure also shows that there is no dependence on metallicity among the rejected objects. The list of the 15 rejected objects is shown in Table C.1 in Appendix C.

Fig. 3. Reliability of parallaxes based on the two parameters RUWE and astrometric_gof_al (in absolute value). For reference, the typical limits at 1.4 and 12.5 have been reported with dashed lines.

3. Analysis

The fitting method used in this paper is the same as in Ripepi et al. (2022a), which, in turn, is based on the approach by Riess et al. (2021). We first defined the photometric parallax (in mas) as

$$\begin{array}{c}\hfill {\varpi }_{\mathrm{phot}}={10}^{-0.2(w-W-10)},\end{array}$$(2)

where w is a generic apparent Wesenheit magnitude (see Table D.1), while W is the absolute Wesenheit magnitude, which can be written as

$$\begin{array}{c}\hfill W=\alpha +\beta ({log}_{10}P-1.0)+\gamma [\mathrm{Fe}/\mathrm{H}].\end{array}$$(3)

Note that, contrary to our previous work, given the reduced number of objects compared with Trentin et al. (2024a), we did not try to study the metallicity dependence of the slope, i.e. of the log P coefficient. This is also justified because we calculated the zero-point (ZP) counter-correction directly from the data, and thus already added a fourth parameter to the fit. In any case, the analysis of the slope dependence on metallicity is postponed to a future paper, which will include the complete C-MetaLL sample, consisting of more than 400 DCEPs.

Indicating with ϖ_DR3 the DR3 parallax already corrected for the ZP bias (see Sect. 2.5), we wish to minimise the following quantity:

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum \frac{{({\varpi }_{DR3}-{\varpi }_{\mathit{phot}}+ϵ)}^2}{{\sigma }^2}}.\end{array}$$(4)

Here σ is the total error obtained by summing up in quadrature the uncertainty on ϖ_DR3 and ϖ_phot: $\sigma =\sqrt{{\sigma }_{{\varpi }_{DR3}}^2+{\sigma }_{{\varpi }_{\mathrm{phot}}}^2}$, while ϵ is the parallax ZP counter-correction. In more detail, the astrometric σ_ϖDR3 consists of two contributions: (i) the standard error of the parallax as reported in the DR3 catalogue, which we conservatively increased by 10%; and (ii) the uncertainty on the individual ZP corrections, assumed to be 5 μas (see e.g. Lindegren et al. 2021; Riess et al. 2021). The uncertainty on ϖ_phot is more complex to evaluate. Using the equivalence δϖ/ϖ=δD/D, where D denotes the distance, and recalling the definition of the distance modulus μ = −5 + 5log₁₀D, error propagation and straightforward algebra yield

$$\begin{array}{c}\hfill {\sigma }_{{\varpi }_{\mathrm{phot}}}=0.46\times {\sigma }_{\mu }\times {\varpi }_{\mathrm{phot}},\end{array}$$(5)

where ${\sigma }_{\mu }=\sqrt{{\sigma }_w^2+{\sigma }_W^2}$. The term σ_w is readily computed by propagating the uncertainties on the single magnitudes in the Wesenheit definitions of Table 1. In contrast, estimating σ_W is more challenging, as it requires prior knowledge of the intrinsic scatter of the adopted relation. To this end, we followed the approach by Riess et al. (2021), adopting a scatter of 0.06 mag in the NIR bands, while in the optical bands we assumed a more conservative value of 0.1 mag. As is discussed in Ripepi et al. (2022a), this is justified on theoretical grounds (e.g. De Somma et al. 2020, and references therein). Note also that Eq. (5) requires the knowledge of ϖ_phot, and, in turn, of w and W. While w is a measured quantity, W includes the unknowns α,  β,  γ. Therefore, the minimisation procedure detailed below proceeded iteratively:

• To obtain a reasonable starting point for the parameter estimation (priors), we first minimised the standard χ² function described in Eq. (4). The resulting values of α,  β,  γ, and ϵ were used as inputs in the subsequent steps. During this first step, we checked the sample for large outliers, which could indicate undetected problems with the astrometry or photometry. To do this, we applied a sigma-clipping procedure, identifying stars whose weighted difference between the calculated and observed parallaxes exceeded 3σ. This method flagged eight problematic stars (listed in Table C.2): four of them were deviant in almost all PW relations, likely due to incorrect parallaxes, and the other four were outliers only in the optical relations, likely due to poor photometry in those bands. These stars were subsequently removed from the derivation of any PW relation in which they were identified as outliers (see Table 2).

• To further investigate the statistical properties of our dataset, we examined the residuals of the standard χ² fit described in the previous point. Taking as an example the WHVI magnitude, we find that the residuals exhibit significant excess variance throughout the sample, with a reduced chi-squared of χ²/d.o.f.⁶ ≈ 2.6 and AIC ≈ 700. This indicates that the scatter in the data is approximately 1.6 times larger than the formal uncertainties⁷. Furthermore, the distribution of residuals is non-Gaussian, displaying tails with significantly more high-significance deviations than expected. To address this over-dispersion from biasing our parameter estimation, we adopted a robust loss function based on the Cauchy distribution⁸. The Cauchy loss is defined as

  $$\begin{array}{c}\hfill {\mathcal{L}}_{\text{Cauchy}}=\frac{c^2}{2}{\displaystyle \sum _i}log(1+{\left(\frac{r_i}{c}\right)}^2),\end{array}$$(6)

  where r_i = (ϖ_DR3 − ϖ_phot + ϵ)²/σ² is the normalised residual, and c is a scale parameter that controls the influence of outliers. The parameter c was estimated from the residuals of the initial chi-square fit using the median absolute deviation (MAD), scaled to be consistent with a normal distribution: c = 1.4826 × MAD. This formulation provides a smooth, non-quadratic penalty for large residuals, limiting the influence of potential outliers without requiring explicit clipping.

• To estimate the posterior distributions of the parameters α,  β,  γ, and ϵ, we adopted a Bayesian approach and sampled from the posterior using the affine-invariant Markov Chain Monte Carlo (MCMC) sampler emcee (Foreman-Mackey et al. 2013). The log-likelihood was defined as

  $$\begin{array}{c}\hfill log\mathcal{L}(\theta )=-\frac{1}{2}{\mathcal{L}}_{\text{Cauchy}}(\theta ),\end{array}$$(7)

  where θ = (α,  β,  γ, ϵ) is the vector of model parameters. We initialised a set of 50 walkers around the initial solution and evolved them for 5000 steps, discarding the burn-in phase and thinning the chain to obtain a clean sample of the posterior. Tests using a larger number of walkers (up to 1000) showed no significant improvement, so we adopted 50 to optimise computational efficiency. Our approach enables the estimation of parameter uncertainties and degeneracies, as well as robust inference in the presence of non-Gaussian errors.

Our best estimation of each parameter is represented by the median of the posterior distribution, sampled via the MCMC. We set the 16th and 84th percentiles as uncertainties.

To test the soundness of the procedure outlined above, we tried to reproduce the results by Riess et al. (2021). Indeed, we used the same method as in that work, the only difference being a different minimisation technique. As is shown in detail in Appendix D, our procedure gives the same numbers as in the literature.

4. Results

Based on this positive test, we obtained the PWZ relations for the selected Wesenheit magnitudes defined in Table 1. Two examples of the results for the WJK_S and WcHVI are shown in Fig. 4, where we show the posterior distribution of the four (α,  β,  γ, ϵ) parameters. The comparison between the Gaia DR3 and the photometric parallaxes for the same selected Wesenheit magnitudes is shown in Fig. 5. Here, the objects considered outliers by the Cauchy-loss functions are highlighted in red. We recall that the adopted procedure assigns less weight to the outlying points but does not remove them. The figure shows the robustness of our fits, and in particular that the distribution of residuals shows no significant, unmodeled dependence on metallicity, further supporting the consistency of the baseline fit.

Fig. 4. Examples of corner plot with the posterior distributions for the output parameters and the best-fit solution obtained for the WJKS (top) and WcHVI (bottom) Wesenheit magnitudes. The units of α,  β,  γ,  ϵ and of their uncertainties are mag, mag/dex, mag/dex, and mas, respectively.

Fig. 5. Comparison between final photometric parallax (i.e. the one calculated with the final result of the fitting procedure) and that from Gaia. In filled circles and filled stars are the objects not considered or considered as outliers by the Cauchy loss function. All the data are colour-coded according to their [Fe/H] ratio. The top and bottom panels show the results for the WJKS and WcHVI Wesenheit magnitudes, respectively. The units of α,  β,  γ,  ϵ and ϖ are mag, mag/dex, mag/dex, mas, and mas, respectively.

Applying the procedure to all the Wesenheit magnitudes leads to the derivation of the parameters, which are shown in Fig. 6 and listed in Table 2. The figure displays the results for Wesenheit magnitudes with increasing wavelength from left (optical) to right (NIR). We note that:

• α and β have similar behaviour. As expected, their values tend to increase (in absolute value) towards the NIR. However, the trend is reversed when considering the Gaia and VI (OGLE) Wesenheits. In particular, there is a significant difference between the α and β values for the VI Wesenheits defined by OGLE and those calculated from the Fitzpatrick (1999) extinction law. This does not occur for the JK_S owing to the smaller difference between the coefficients that multiply the colour in this case.

• γ does not vary significantly from the optical to the NIR, even though it approaches the literature values in WJK (see shaded region in Fig. 6), but remains more negative in WVK and WcHVI.

• The Gaia parallaxes ZP counter-correction tends to vary systematically from ∼16 μas to a few μas from the optical to the NIR, having the lowest value for the WJK_S. The average value of about 10 μas is marked in the figure for reference. For instance, this is the value obtained for WcHVI.

• The values of γ and ϵ are correlated. This means that, for example, for the WcHVI used in the distance scale, passing from ϵ = 10 μas as found from our fitting procedure to ϵ = 14 μas (Riess et al. 2021) would increase γ to ∼ − 0.35 mag/dex. Naturally, since other correlations are present, this would also imply α ∼ −6.04 mag, increasing the discrepancy of the ZP compared with Riess et al. (2021, see also Appendix B)

Fig. 6. Values of the four parameters, α,  β,  γ, and ϵ, for the different Wesenheit magnitudes considered in this work. Red dots show the results obtained by adopting the Fitzpatrick (1999) reddening law (F99 in the label). The shaded area comprises the range −0.4 < γ < −0.2 mag/dex, where most of the recent literature results are located. The units of α,  β,  γ,  ϵ and of their uncertainties are mag, mag/dex, mag/dex, and mas, respectively.

As in our previous papers, we validated the distances derived from these PWZ relations by comparing them with the 1% precise geometric distance of the LMC derived by Pietrzyński et al. (2019). To this aim, we adopted a sample of about 4500 DCEPs in the LMC that have optical photometry from Gaia and OGLE IV surveys as well as J, K_S photometry from the Vista Magellanic Cloud survey Cioni et al. (2011, VMC); Ripepi et al. (2022b, VMC) and H data from Inno et al. (2016). Similarly to the Galactic sample, we used these H,  V,  I data to calculate the corresponding HST Wesenheit magnitudes. Then, we applied the derived PWZ relations to this LMC sample, adopting [Fe/H] = −0.41 dex (Romaniello et al. 2022), and thus derived a distance modulus from each star.

The median of the resulting distance moduli distribution provides us with an estimate of the LMC distance. The calculation of its error is discussed in Appendix E. The derived LMC distances for each PWZ relation are displayed in Fig. 7 and listed in Table 2. The resulting μ_LMC values are typically about 2–5% smaller than the geometric distance by Pietrzyński et al. (2019) except for WG, which provides a ∼4% farther distance. In any case, given the typical 5.5% uncertainty, almost all the Wesenhit magnitudes analysed in this work provide LMC distances well within 1σ compared to the geometric distance of the LMC. We note that the best results are provided by the NIR Wesenheit magnitudes, especially the WVKs one. Overall, this validation procedure suggests that the results presented in this section are reliable. We consider the PWZ in Table 2 to be our baseline results.

Fig. 7. LMC distance moduli calculated from our best-fitting PWZ relations. Red dots show results relative to the Wesenheit definition according to the Fitzpatrick (1999) reddening law. The shaded region corresponds to the geometric distance of the LMC by Pietrzyński et al. (2019) ±1σ.

5. Variations from the baseline

To explore possible biases in our results, which could lead to the large metallicity dependence of the PWZ relations in the previous section, we re-ran the calculations by applying possible variations described in the following sections. In the next sections, we use the Akaike information criterion (AIC) and the Bayesian information criterion (BIC) to compare different model calculations (see Appendix F for a detailed definition of these quantities) that involve samples of approximately the same size.

5.1. Use of pure σ-clipping method

To be consistent with the techniques adopted in previous works, in this variation of the baseline, instead of using the Cauchy-loss technique, we applied a pure σ-clipping method. The results of this procedure for the case with 3σ and 5 iterations are listed in the top part of Table G.1 and shown in Fig. 8. We experimented by lowering the sigma-clipping to 2.5 and 2σ, and/or increasing the number of iterations, without significant changes in the coefficient values, but with increasingly smaller errors, as expected.

Fig. 8. Comparison between the PWZ parameters obtained with our baseline sample (F+1O pulsators, using [Fe/H], black squares) and those obtained: (i) using a pure σ-clipping algorithm (red circles); (ii) using only F-mode pulsators with [Fe/H] as abundance indicator (blue triangles); (iii) using F+1O pulsators with the total metallicity [M/H] as abundance indicator (green diamonds); and (iv) using only objects with G < 12.5 mag (violet stars). The shaded area comprises the range −0.4 < γ < −0.2 mag/dex, where most of the recent literature results are located. The units of α,  β,  γ,  ϵ and of their uncertainties are mag, mag/dex, mag/dex, and mas, respectively.

In general, the main difference from the baseline is a modest variation (both increasing or decreasing) in the γ terms at different wavelengths, accompanied by a general reduction in errors. The distance of LMC is essentially the same as in the baseline, albeit with smaller errors. On the contrary, the AIC and BIC values are significantly larger. However, these goodness of fit values are not directly comparable with the baseline, because in this variant we minimised the χ², while in the baseline we minimised the Cauchy-loss function. Overall, this experiment shows that the introduction of the Cauchy-loss algorithm does not introduce significant differences with more standard procedures.

5.2. Use of F-mode pulsators only

To investigate the impact of using 1O-mode pulsators with fundamentalised periods together with F-mode DCEPs, we re-ran the PWZ calculations using only the latter. Of the 183 F-mode DCEPs in our sample, only 173 have usable parallaxes. The results of the calculations are shown in Fig. 8 and in the middle-top part of Table G.1. The adoption of the reduced sample essentially produces: (i) errors larger by a factor of 1.5 or more on all the parameters; (ii) γ values much more negative (larger in absolute value) compared with the baseline calculation; (iii) parallax ZP counter-correction much larger (even positive) than the baseline PWZs; and (iv) low values for the LMC distance modulus. We interpret these findings as resulting from the reduced sample size. Indeed, since the overall relative error on the parallaxes for most of the C-MetaLL sample stars is of the order of 10–15%, with an increase towards the most metal-poor objects (see Fig. 2), we need to use a larger sample to compensate for the reduced precision of the parallaxes with the statistics. In other words, currently, using the F-mode pulsators alone is not feasible.

5.3. Use of the total metallicity [M/H] instead of [Fe/H]

It is usual to express the metal abundance in terms of [Fe/H], i.e. the difference of the logarithm of the iron over hydrogen abundance ratios measured in the atmosphere of the target star and the Sun. This is based on the assumption that the solar heavy element distribution is universal. However, we have known for a long time that this is not always the case. For example, metal-poor stars in the Galactic halo show various degrees of enhancement of the so-called α elements (e.g. O, Ne, Mg, Si, S, Ca, and Ti) compared with the Sun (see e.g. Chapter 8 of Salaris & Cassisi 2005, for a general discussion). The amount of this enhancement can be of the order of [α/Fe] ∼ 0.3–0.4 dex. In these cases, one should consider using the total metallicity [M/H] instead of [Fe/H]. As was demonstrated by Salaris et al. (1993), the net effect of including the [α/Fe] enhancement is to simulate a larger metallicity that, according to Salaris et al. (1993), can be computed from the relation

$$\begin{array}{c}\hfill [\mathrm{M}/\mathrm{H}]\sim [\mathrm{Fe}/\mathrm{H}]+{log}_{10}(0.638\times {10}^{[\alpha /\mathrm{Fe}]}+0.362).\end{array}$$(8)

Typically, in the MW halo, the [α/Fe] enhancement becomes significant at [Fe/H]< –0.6 dex (this number varying in different environments, see e.g. Tolstoy et al. 2009). Since our sample includes objects more metal-poor than the threshold, we tested the effect of the [α/Fe] enhancement on the determination of the PWZ relations for DCEPs using the [α/Fe] values published by Trentin et al. (2024b) for all the 290 stars in the C-MetaLL sample (their fig. 6). Although our DCEPs are disc and not halo stars, their [α/Fe] enhancement can reach values of ∼+0.4 dex. We therefore calculated a new series of PWZ relations using [M/H] instead of [Fe/H] as an abundance indicator. To this aim, we converted [Fe/H] to [M/H] using Eq. (8) and re-ran the calculations on the F+1O sample. The result of this test is shown in Fig. 8 and Table G.1 (middle-bottom part). The effect of using [M/H] provides results and errors close to the baseline. However, we obtain γ and ϵ values which are smaller (more negative) and larger than the baseline, respectively. Thus, the large (in the absolute sense) values of γ obtained from our baseline calculations cannot be caused by the use of [Fe/H] instead of [M/H]. Since the changes of γ and ϵ compensate for each other, the calculated LMC distance moduli are about 1% larger than the baseline⁹, and thus more in agreement with the geometric value. The AIC and BIC values are close to those of the baseline in the optical but slightly lower for the NIR Wesenheit magnitudes. In the latter cases, therefore, the present variation suggests a slightly better fit to the data.

5.4. Use of the brightest DCEPs only

In this variation, we explore the impact of restricting the sample to only the brightest DCEPs, i.e. those with the most precise parallaxes. To avoid reducing the sample size excessively, we adopt a magnitude limit of G < 12.5 mag, which yields 142 usable objects, slightly more than half of the full usable sample of 275 stars. Applying such a magnitude cut introduces well-known observational biases, notably the Malmquist bias¹⁰. Therefore, this test should be regarded as a diagnostic exercise aimed at better understanding the properties of the DCEP sample used in this work.

A further consequence of the magnitude cut is a significant reduction in the metallicity range covered by the sample, as the fainter – typically more distant – DCEPs also tend to be more metal-poor. The result of this variation is shown in Fig. 8 and in Table G.1 (bottom part). While large uncertainties are introduced due to the smaller sample size, it is evident that the magnitude cut leads to a general (absolute) decrease in the γ values, whereas the other parameters remain nearly unchanged and the average distance modulus of the LMC approaches the geometric one by Pietrzyński et al. (2019). It is unclear whether this shift in γ, bringing the values closer to those found by the SH0ES collaboration, is driven by parallax precision or the reduced metallicity range. However, a combination of both effects is likely.

We also attempted the same test on the complementary faint sample. Still, the results were inconclusive due to the low precision of the parallaxes and the limited metallicity range, this time biased towards more metal-poor stars. Some experiments in this direction for the WcHVI magnitudes are listed in Appendix I.

5.5. Binning of the data in [Fe/H]

To further investigate the origin of the large metallicity dependence of our PWZ relations, we decided to estimate the γ term similarly to that adopted by Breuval et al. (2024). These authors assume a universal slope for the HST WHVI magnitudes and then determine the intercept of the simple PW relations in different environments using geometric distances; for example, Gaia parallaxes in the MW, and eclipsing binaries in the LMC (Pietrzyński et al. 2019) and SMC (Graczyk et al. 2020). Finally, they plot the values of these intercepts against [Fe/H] in the three environments. The slope of the linear regression provides the value of γ (see e.g. Fig. 8 in Breuval et al. 2024).

Besides the assumption about the universal slope, the subsequent approximation in this approach is to consider each environment as ‘mono-metallic’. This is justified, based on the present data, in the LMC and SMC, where the dispersion in [Fe/H] of the measured DCEPs is below 0.1 dex (see Romaniello et al. 2022, 2008, for the LMC and SMC, respectively). However, this is less accurate in the MW, where a wide range of abundances is present, thanks to the disc metallicity radial gradient (e.g. Luck & Lambert 2011; Genovali et al. 2014; Kovtyukh et al. 2022; Trentin et al. 2024b). In any case, to try to simulate this approach, we binned our C-MetaLL DCEP sample with usable parallaxes (excluding the outstanding outliers, as explained in Sect. 3) into three metallicity bins, including almost the same number of pulsators (about 90) with average [Fe/H] = −0.66,   − 0.27,   + 0.0 dex and dispersions of 0.16, 0.10, and 0.12 dex, respectively. For each bin, we then carried out the photometric parallax χ² minimisation using the Cauchy loss function, as described above (Eq. 4), removing the γ term and fixing both the slope, β, and the Gaia parallax ZP counter-correction, ϵ. In this way, only α is constrained by the χ² minimisation. We carried out the calculations using the value of ϵ = 10 μas found from our baseline results in the WcHVI case (see Table 2) and two values for β; namely, 3.21 and 3.25. The former comes from our baseline calculation, while the latter is the one used by the SH0ES team. The result of this exercise is displayed in Fig. 9, where we plot the estimated values of α (which still includes the effect of metallicity) versus the binned metallicity [Fe/H] with the respective dispersions. The three data points in the figure have been fitted with a regression line using the ODR (orthogonal distance regression) algorithm as implemented in the Scipy package (Virtanen et al. 2020). This method was chosen because it is robust, and errors on both axes can be accounted for.

Fig. 9. Estimate of the γ parameter for the WcHVI Wesenheit magnitude obtained by binning our data in three [Fe/H] intervals (filled blue circles). The top and middle panels show the results for two different choices of the β parameter (see labels). The bottom panel shows the results for the same slope as in the middle panel but for a different value of ϵ. For comparison, filled magenta circles show the results from Breuval et al. (2024).

The γ values derived this way (ranging from γ ∼ −0.36 to ∼ − 0.44 mag/dex) are smaller in absolute value than those obtained by fitting all PWZ parameters simultaneously (∼ − 0.5 mag/dex). This makes them more similar to the value from Breuval et al. (2024) (∼ − 0.234 ± 0.052 mag/dex), although ours are still more negative by about 1–2σ. Unsurprisingly, the least discrepant γ value (a difference of slightly more than 1σ) was obtained using the β and ϵ values adopted by the SH0ES group (bottom panel).

Using the three selected pairs of β and ϵ, the average differences between our results and those of Breuval et al. (2024) are consistent with zero within the errors for ϵ = 10 μas (see top and middle panels of Fig. 9). This difference increases to a barely significant 2% for ϵ = 14 μas (bottom panel of Fig. 9).

Therefore, we seek to understand the cause of the discrepancy in γ. Using the regression lines in the three panels of Fig. 9 as references, we can evaluate the difference between our solution and the Breuval et al. (2024) data. The top and middle panels (obtained with ϵ = 10 μas) show that the largest deviation between Breuval et al. (2024) and the regression line is about ±0.07 mag. This occurs at the most metal-rich and metal-poor bins (representing MW and SMC galaxies), while the LMC bin aligns almost perfectly with our solution. This alignment explains the near-zero average difference but is difficult to interpret. While the difference in the metal-poor regime could plausibly be ascribed to less accurate Gaia parallaxes for these fainter pulsators, the discrepancy in the MW data point is harder to explain, as parallaxes in this metallicity range should be the least affected by systematics. The result for ϵ = 14 μas and β = −3.25 mag/dex (the values adopted by Breuval et al. 2024) is even more complex. In this case, the discrepancy in the most metal-poor (SMC) bin nearly vanishes, while the disagreement for the other two bins appears to be enhanced, even though the magnitude of the deviation is smaller than before (∼0.05 mag). This accounts for the small, barely significant non-zero difference between our solution and the Breuval et al. (2024) data, and for the reduced scatter (∼2%) compared to the previous cases (∼4%).

The cause of this complex behaviour is not entirely clear. The main conclusion we can draw is that accurately measuring γ is a difficult task. Indeed, Fig. 9 demonstrates how changes of just a few hundredths of a magnitude at the extremes of the metallicity distribution are sufficient to significantly alter the value of γ. We conclude that a decisive step forward in establishing the correct value of γ requires two things: first, the addition of more homogeneous DCEP data across the entire metallicity range (which we aim to provide in future C-MetaLL data releases), and second, and more importantly, the expected reduction in systematic errors in the upcoming Gaia DR4 parallaxes.

5.6. Variable parallax zero point counter-correction

Some literature papers (e.g., Riess et al. 2022a; Khan et al. 2023) have reported a possible dependence of the Gaia parallax ZP counter-correction – ϵ in our case – on the brightness of the specific sample. The general result is that there is a sort of linear decrease (in an absolute sense) in ϵ going towards fainter magnitudes, up to approximately zero at G ∼ 16 mag. It is worth noting that there is no consensus regarding this occurrence, as has been discussed in other papers (e.g. Molinaro et al. 2023; Groenewegen 2024). In any case, to not leave anything unexplored, considering that our targets span a range of about six magnitudes, we substituted in Eq. (4) the constant ϵ with a linear equation as a function of the magnitude G: z + p × G and re-ran the calculation. The idea was to verify if the release of the assumption of a constant ϵ could be responsible for the large γ values we find in the C-MetaLL survey. The results of this procedure are reported in the first part of Table H.1. There is no effect on the values of γ, and the values of the slopes, p, are small in all the cases and barely significant at a 2σ level in the NIR bands only. In any case, the solution appears valid, as the distance of the LMC, albeit systematically larger than the geometric one in all the cases, is consistent with this within the errors. Concerning the AIC and BIC values, we do not observe any significant decrease compared with the baseline, as expected in the case of better modelling. On the contrary, the BIC value, which tends to penalise the overfitting of the data, is slightly increased. We conclude that the dependence of the Gaia parallaxes ZP counter-correction on brightness, if any, is small and not able to explain the large metallicity dependence we find in our work.

5.7. Non-linear dependence on [Fe/H]

To explore a non-linear dependence of the metallicity term (γ), we can introduce higher-order terms for the metallicity into the model. The most straightforward approach within the current polynomial-fitting framework is to include a quadratic term for the metallicity, [Fe/H]². Therefore the model is now the following:

$$\begin{array}{c}\hfill W=\alpha +\beta ({log}_{10}P-1.0)+\gamma 1[\mathrm{Fe}/\mathrm{H}]+\gamma 2{[\mathrm{Fe}/\mathrm{H}]}^2,\end{array}$$(9)

where γ1 and γ2 are the linear and quadratic terms in metallicity. As in the previous case, we now have five parameters to fit. The result of this exercise is reported in the second part of Table H.1. We note that the linear term remains large and of the order of −0.5 mag/dex in all the cases, while the quadratic term is only significant at 1σ. The values of AIC and BIC are increased compared with the baseline, indicating that this quadratic model does not describe the data in a better fashion.

Further investigations about the value of γ for the WcHVI have been extensively carried out in Appendix I and are not reported here for brevity. The main results are discussed in the next section.

6. Discussion and conclusions

In this work, we have exploited the homogeneous spectroscopic abundances provided for 290 DCEPs in the context of the C-MetaLL project (Trentin et al. 2024b). We collected or calculated intensity-averaged magnitudes for these stars in various bands, which we used to calculate several optical and NIR Wesenheit magnitudes. We adopted both the Cardelli et al. (1989) and Fitzpatrick (1999) reddening laws. Empirical relations from Riess et al. (2021) were employed to transform Johnson-Cousins V,  I and H,  V,  I Wesenheit magnitudes in the respective HST F555W,  F814W and F160W,  F555W,  F814W kins. Our database was completed by: (i) periods and modes of pulsation from the literature (Gaia and OGLE IV surveys); and (ii) parallaxes from the Gaia mission corrected for the individual ZP bias Lindegren et al. (2021). This wealth of data was then used to derive new PWZ relations from 275 DCEPs with usable parallaxes, utilising the photometric parallax technique and a new minimisation technique that can handle outlier measurements without removing the data. The fitting procedure also calculates the unknown global ZP counter-correction for the Gaia parallaxes (ϵ).

Compared with our previous work, here we adopted a DCEP sample with homogeneous abundance determinations, an almost even population of objects over the entire metallicity range spanning approximately +0.3 < [Fe/H] < −1.1 dex, the widest ever used in such calculations. Moreover, we adopted a new, more accurate fundamentalisation relation for the 1O pulsators and a new robust photometric parallax technique based on the MCMC algorithm, which allowed us to determine all fitting parameters simultaneously, including the value of ϵ directly from the data. In particular, we adopted a Cauchy likelihood, which enabled us to safely handle the substantial excess variance observed in our full dataset.

Our findings can be summarised as follows:

• The metallicity dependence, γ, is confirmed to be larger in absolute value compared with the recent findings (Breuval et al. 2022, 2024), with typical values of γ ∼ −0.5 mag/dex in the optical and HST Wesenheit bands and γ ∼ −0.4 mag/dex in the NIR bands (e.g. using J,  K_S). On the other hand, our results agree well in the NIR bands with those by Wang & Chen (2025), who used a DCEP sample and fitting technique similar to that of the present and our previous papers.

• The Gaia parallax ZP counter-correction (ϵ) varies almost monotonically from the Gaia bands towards the NIR, spanning from about 16 to 7 μas with typical errors of about 4 μas. In particular, for the WcHVI magnitude, we obtain 10±4 μas, which agrees within 1σ with Riess et al. (2021).

• The application of our PWZ relations to 4500 DCEPs in the LMC provides independent distances to this galaxy, which are typically within 1σ of the geometric measurement by Pietrzyński et al. (2019), with the remarkable exceptions of the WVI and WcVI magnitudes. Overall, this consistency demonstrates the soundness of our work.

• The inclusion of the Fitzpatrick (1999) reddening law in place of Cardelli et al. (1989)’s in the V,  I and J,  K_S Wesenheit definition does not impact the results significantly.

• The use of only F-mode pulsators yields significantly larger errors in the parameters and unreliable results. A larger number of F-mode pulsators is needed to obtain results comparable with the entire F+1O sample.

• The adoption of the total metallicity [M/H] instead of [Fe/H] by correcting the latter with the [α/Fe] measurements provided by Trentin et al. (2024b) has the effect of increasing in an absolute sense the values of γ, while reducing those of ϵ.

• The selection of only the brightest DCEPs (G < 12.5 mag), which have on average better parallax precision, provides values of α,  β,  ϵ that follow the baseline’s behaviour for the different Wesenheit magnitudes, while γ generally decreases in absolute value to ∼ − 0.3; −0.4 mag/dex, i.e. closer to SHOES’s team results. This occurrence may be driven by the improved parallax precision or the reduced metallicity range, though a combination of both effects is probable. This procedure, however, besides reducing the metallicity range, also introduces Malmquist bias, and therefore has to be regarded merely as an experiment.

• We investigated the potential dependence of the Gaia parallax ZP counter-correction (ϵ) on stellar brightness by using a linear correction (z + p × G). This procedure was implemented to check if releasing the assumption of a constant ϵ could account for the large metallicity dependence (γ) found in our C-MetaLL survey. The results show that the values of γ are unaffected, and the slope (p) is negligible, confirming that any brightness dependence cannot explain our findings.

• The introduction of a quadratic metallicity term, γ₂[Fe/H]², into the period-luminosity relation (Eq. (9)) did not significantly improve the model fit, as is indicated by the increased AIC and BIC values compared to the baseline linear model. While the linear metallicity term (γ₁) remains large (≈ − 0.5 mag/dex), the quadratic term is only significant at the 1σ level. This suggests that a non-linear (at least a quadratic) dependence on metallicity is not strongly supported by the current dataset, and the simpler linear model provides a better description of the data.

• The binning of our data into three metallicity intervals to simulate the data used by Breuval et al. (2024) and adoption of their fitting technique yield values of γ (∼ − 0.36 to ∼ − 0.44 mag/dex) that are closer to, but still 1–2σ more negative than, the results from Breuval et al. (2024). This discrepancy is complex to understand, as small magnitude differences (0.05–0.07 mag) at the metallicity extremes (SMC and MW) are sufficient to significantly change the fitted γ. We conclude that accurately measuring γ is extremely difficult and requires more homogeneous data and improved Gaia DR4 parallaxes to reduce systematic errors.

• As is shown in Appendix I, smaller (absolute) values of γ (∼ − 0.3 mag/dex) for the WcHVI relation are obtained for a sub-sample of DCEPs that are either brighter than G ≈ 11.5–12.5 mag, located at distances of < 3/4 kpc (and thus benefitting from more precise parallaxes), or exhibit metallicities characteristic of solar neighbourhood stars (approximately −0.5/0.7 < [Fe/H] < +0.1/0.2 dex). If the larger value of γ is not attributable to increased parallax uncertainties in more distant (d > 3/4 kpc), metal-poor ([Fe/H] < −0.7 dex) stars, this could suggest either an intrinsically higher metallicity coefficient or a non-linear dependence of DCEP luminosities on metallicity at the low-metallicity end. Quantifying this effect with high precision remains challenging due to the significant covariance between the PWZ zero point, γ, and the adopted parallax ZP offset correction.

Overall, the results presented in this paper suggest that the exact metallicity dependence of PW relations for DCEPs remains uncertain. From a theoretical point of view, pulsational predictions are in favour of a mild value of γ ∼ −0.2 mag/dex (e.g. Anderson et al. 2016; De Somma et al. 2022; Khan et al. 2025), similar to that provided by the empirical technique adopted by Breuval et al. (2021, 2022, 2024). However, in this paper, we have shown that our fitting procedure is solid (see Appendix D); therefore, the problem seems to be in the derivation of the γ value using different geometric distance indicators, such as in Breuval et al. (2024) or using uniquely parallaxes from Gaia. Both approaches have pros and cons and are prone to different systematic uncertainties, such as the use of a universal slope and an average metallicity in each galaxy in the former approach or the bias on the Gaia parallax in the latter.

As was stressed in the introduction, the size of the metallicity dependence of DCEPs’ PL and PW relations has important astrophysical consequences. Perhaps the most notable one is that a larger (absolute) value of γ than the one currently used in the cosmic ladder would go in the direction to diminish the value of the inferred H₀ by 1–2% for e.g. γ ∼ −0.4 mag/dex. However, this reduction, far from solving the Hubble tension, would just reduce the discrepancy by about 1σ.

In the context of the C-MetaLL project, we expect to improve our results shortly based on the measurement of abundances for about additional 100 DCEPs for which the analysis is still ongoing as well as the release of homogeneous g,  r,  i,  z,  J,  H,  K_S rapid-eye-movement telescope¹¹ photometry for more than 120 DCEPs (see Bhardwaj et al. 2024, for early results). However, a significant improvement in this topic is expected with the next release of the Gaia mission, Data Release 4, which will provide parallaxes that are not only more precise but also most likely less affected by systematic effects.

Data availability

Full Table A.1 is available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/708/A216.

Acknowledgments

We thank our anonymous referee for their insightful comments. We also warmly thank A. Riess and L. Breuval for very useful discussions and suggestions, which helped us to improve our work. We acknowledge funding from: INAF GO-GTO grant 2023 “C-MetaLL – Cepheid metallicity in the Leavitt law” (P.I. V. Ripepi); PRIN MUR 2022 project (code 2022ARWP9C) ‘Early Formation and Evolution of Bulge and Halo (EFEBHO)’, PI: Marconi, M., funded by the European Union – Next Generation EU; Large Grant INAF 2023 MOVIE (P.I. M. Marconi). AB thanks funding from the Anusandhan National Research Foundation (ANRF) under the Prime Minister Early Career Research Grant scheme (ANRF/ECRG/2024/000675/PMS) This research has made use of the SIMBAD database operated at CDS, Strasbourg, France. G.D.S. acknowledges funding from the INAF-ASTROFIT fellowship (PI G.De Somma), from Gaia DPAC through INAF and ASI (PI: M.G. Lattanzi), and from INFN (Naples Section) through the QGSKY and Moonlight2 initiatives. 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 was supported by the Munich Institute for Astro-, Particle and BioPhysics (MIAPbP), which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. This research was supported by the International Space Science Institute (ISSI) in Bern/Beijing through ISSI/ISSI-BJ International Team project ID #24-603 – “EXPANDING Universe” (EXploiting Precision AstroNomical Distance INdicators in the Gaia Universe).

References

Appendix A: Data used in this paper

This appendix includes an extract of the table including all the measurements used in this paper.

Appendix B: Comparison between the native and ground-based HST magnitudes

In this paper, we converted ground-based magnitudes in the H, V, I bands into the corresponding HST ones F106W,  F555W,  F814W, adopting the transformations provided by Riess et al. (2021). To verify the accuracy of our transformed magnitudes cF106W, cF555W, and cF814W and relative Wesenheit cWVI and cWHVI, we considered a sample of 69 DCEPs with native HST photometry published by Riess et al. (2021)¹² and proceeded as follows: 1) calculated V, I from the Gaia DR3 G, G_BP and G_RP through the transformation by Pancino et al. (2022); 2) obtained H band from Trentin et al. (2024a); 3) calculated cF106W, cF555W, and cF814W and relative Wesenheit cWVI and cWHVI with the quoted transformation; 4) compared each calculated magnitudes and Wesenheit with the native HST data. The results of this test are reported in Fig. B.1 for cF106W, cF555W, and cF814W and Fig. B.2 for cWVI and cWHVI. This test shows the absence of significant discrepancies between our way of calculating the corrected HST photometry and the native HST one.

Fig. B.1. Comparison between the F106W, F555W, and F814W magnitudes observed with HST by Riess et al. (2021) and those calculated based on ground-based H, V, I magnitudes transformed into the corresponding HST magnitudes using the conversions by Riess et al. (2021).

Fig. B.2. As in Fig. B.1 but for the WVI and WHVI Wesenheit magnitudes.

Appendix C: Properties of the sample in comparison with literature

Table C.1 lists the stars which have been excluded from the PW calculations due to the poor quality of the astrometric solution (see Sect. 2.5 for details).

Figure C.1 shows the different locations of the bulk data in our previous paper C-MetaLL IV (Trentin et al. 2024b), including literature values, and the present sample, only based on C-MetaLL spectroscopic data. Similarly, Fig. C.2 displays the distribution of parallax error as a function of both GaiaG mag and [Fe/H] for the C-MetaLL IV (Trentin et al. 2024b) sample. A comparison with Fig. 2 reveals again the more homogeneous distribution in metallicity of the present work sample.

Fig. C.1. Comparison of the period-[Fe/H] distributions between the samples adopted in C-MetaLL IV paper (Trentin et al. 2024b) (blue density distribution) and in the present work (red filled circles).

Fig. C.2. As in Fig. 2, but for the sample adopted in C-MetaLL IV paper (Trentin et al. 2024b).).

Appendix D: Test of the photometric parallax procedure

In this section, we test the ability of our implementation of the photometric parallax method to give the correct results. To this aim, we adopted the Riess et al. (2021) sample of Galactic DCEPs observed with HST (their table 1)¹³ and applied the same selection procedures in terms of Gaia astrometry (through the astrometric_gof_al parameter), including the exclusion of the outliers S Vul and SV Vul. We also used the same width of the instability strip for the WHVI-pagination magnitude equal to 0.06 mag. We then fitted Eq. 4 as explained in Sect. 3. The resulting posterior distributions for the output parameters are shown in Fig. D.1. The results for the four parameters shown in the figure are also listed in Table D.1, with the relative errors. They do not match exactly those by Riess et al. (2021) because these authors adopted a prior for the value of the global counter-correction, while we did not. However, our results match perfectly those obtained by Högås & Mörtsell (2025), which also re-analysed the data by Riess et al. (2021) but left the four fitting parameters free to vary, as we did here (see Table D.1). In addition, this fit provides a χ²/dof ≈ 1.04 (and AIC = 71), as expected for a sample with correctly estimated errors, and in agreement with the value of χ²/dof by Riess et al. (2021). We conclude that our implementation of the photometric parallax fitting is reliable and free of technical errors.

Fig. D.1. Corner plot with the posterior distributions for the output parameters and the best-fit solution obtained for the Riess et al. (2021) sample and all parameters free to vary.

Appendix E: Calculation of the error on the distance of LMC

The total error is calculated by summing the statistical and systematic errors in quadrature: σ_total² = σ_stat² + σ_sys².

The statistical error represents the uncertainty on the mean distance modulus derived from the sample of N DCEPs in the LMC. It is calculated as ${\sigma }_{\text{stat}}=\frac{{\sigma }_{\mu }}{\sqrt{N}}$, where: σ_μ is the standard deviation (1.48 × MAD) of the distance modulus values for the N stars in the sample. Since N is large (> 4000 stars), σ_stat is of the order of 0.005 mag in all the bands, hence it is the systematic uncertainty that dominates the error budget.

The systematic error is the quadrature sum of three main components: σ_sys² = σ_PW² + σ_ϵ² + σ_⟨[Fe/H]⟩², where the different contributions are: i) the uncertainty from the PW relation’s calibrated coefficients errors (σ_α, σ_β, σ_γ), propagated in the usual way; ii) the uncertainty on the W absolute magnitude due to the error on the value of ϵ (σ_ϵ), evaluated at the distance of the LMC and iii) the uncertainty from the LMC’s assumed mean metallicity, which is however very small (Romaniello et al. 2022) and can be safely neglected.

Appendix F: Model Comparison Metrics

We utilise the AIC and the BIC to formally compare the statistical performance of different PWZ models. These criteria balance the goodness-of-fit against the complexity (number of free parameters, k).

For a model where the maximum log-likelihood, lnℒ_max, is estimated by the maximum of the log-probability function ln𝒫_max from the MCMC sampler, the criteria are defined as:

$$\begin{array}{cc}\hfill \text{AIC}& =2k-2ln{\mathcal{L}}_{max}\hfill \end{array}$$(F.1)

$$\begin{array}{cc}\hfill \text{BIC}& =kln(N)-2ln{\mathcal{L}}_{max}\hfill \end{array}$$(F.2)

where k is the number of free parameters in the model (e.g., k = 4 for the baseline model), N is the number of DCEPs used in the fit and lnℒ_max is the maximum log-likelihood, calculated using the median parameters from the MCMC chain: lnℒ_max ≈ −0.5∑(Cauchy Loss)

The model with the lowest AIC or BIC value is preferred. We use the result derived from the MCMC sampling as the basis for these metrics, as it is generally more robust than a simple minimisation result.

Appendix G: Summary of the main variations of the baseline

In this section, we report the main variations of the baseline carried out for all the PW relations investigated in this work. The results of these tests is shown in Table G.1.

Appendix H: Test of variable parallax zero point counter-correction

To test the possible dependence on the GaiaG magnitude of the parallax counter-correction, we modified Eq. 4, substituting ϵ with the linear equation z + p × G, where z and p are the new two unknowns which replace ϵ. The introduction of an additional unknown impacts the stability and robustness of the solution, given that the number of data points remains the same as before. The result of this test is shown in Table H.1.

Appendix I: Further variations for the cHVI Wesenheit magnitude

In this appendix, we present the results of the photometric parallax fit for the cHVI Wesenheit magnitude, including additional variations of the baseline. This procedure is applied only to the cHVI Wesenheit because this is the Wesenheit magnitude used in the extragalactic distance scale. Thus, it is crucial to understand the origin of the larger γ value in our baseline compared with the values used by the SH0ES team. Specifically, we carried out the photometric parallax solution making selections as follows: i) in G magnitude, using only stars brighter than 10.5, 11.5, 12.5, 13.5, 14.5 and 15.7 mag (i.e. the whole sample); ii) in G magnitude, using only stars fainter than 10.5, 11.5 and 12.5 mag (fainter thresholds would comprise too few stars) iii) in parallax signal-to-noise, using 10 and 15 as thresholds; iv) in distance, using 3 and 4 Kpc as thresholds; v) in metallicity, using minimum values of -0.7 and -0.5 dex in [Fe/H] and +0.1, +0.2 and +0.3 dex for the maximum; vi) in period, SNR and metallicity for F-mode pulsators only. In this way, we cut the extreme values of metallicity that may cause extreme values of γ. The results of the entire procedure are reported in the first part of Table I.1. The variations involving only low-quality parallaxes, namely SNR(ϖ) < 10; SNR(ϖ) < 15 and Dist> 4Kpc, or G > 11.5 produce non-realistic results, with intercepts ∼0.3-0.5 mag larger than the acceptable range of values. Values of γ > −0.40 mag/dex, i.e. significantly smaller in an absolute sense than the baseline, can only be obtained with samples including the brightest objects (i.e. G< 12.5 mag) or significantly restricting the metallicity range, particularly using objects with [Fe/H]<  + 0.1 dex.

However, due to the high correlation among the four parameters used in the calculations, it is difficult to determine which sample characteristic produces the large values of γ. Therefore, we decided to fix a couple of parameters, leaving the other two free to vary. We started by fixing α = −5.95 mag and β = −3.25 mag/dex, according to the typical values obtained by the SH0ES team for the LMC. The result of the calculations in this case is listed in the second part of Table I.1. All the selections give high values of γ (with large errors) apart from the same three cases mentioned above. In almost all cases, the values of ϵ are significantly smaller than when all four parameters are free to vary, and compatible with zero in most cases. Remarkably, the resulting LMC distance moduli are consistently smaller than the geometric value, indicating that this choice of alpha and beta is inadequate to describe the data.

The second test consisted of fixing β as before and ϵ instead of α. We choose two values for ϵ, namely −10 μas, as resulted from our baseline calculation, and −14 μas, according to Riess et al. (2021). The results of the test in the two cases are reported in the top and bottom parts of Table I.2, respectively. We note that only the options including low-quality parallaxes give unreliable values of α. Apart from this, both fixed values of ϵ provide a large number of cases with γ < −0.4 mag/dex, though the case of −14 μas provides generally lower (in absolute sense) values of γ, up to ∼ − 0.2 mag/dex, albeit often with significances only slightly larger than 1 σ owing to the reduced samples (about 100 objects or fewer). Notably, the large majority of the solutions give reliable values for the distance to LMC. Concerning the goodness of fit, the values of AIC and BIC are generally moderately larger in the case of the four free-to-vary parameter case. The other three cases show comparable results with a slightly better fitting for β = 3.25 and ϵ = 0.014 mas.

Overall, the experiments discussed above tell us that utilising DCEPs close to the Sun, or, equivalently, bright, or with a reduced range in metallicity, produces values of γ generally significantly lower than those obtained using the entire sample. There are multiple explanations for this occurrence: i) the value of γ is actually as high as -0.5 mag/dex; ii) γ has a non-linear dependence on metallicity (but see Sect. 5.3); iii) the high values of γ are introduced by unknown systematic errors linked to the low-quality parallaxes. With the present dataset, we are not able to determine which scenario is the right one. To this end, increasing the sample of low-metallicity (i.e. fainter objects, with lower-quality parallaxes) DCEPs could help to increase the statistics. Alternatively, the improved Gaia DR4 parallaxes may help in reducing the systematics.

All Tables

All Figures

	Fig. 1. Periods and iron abundances spanned by the investigated sample. The different pulsation modes are labelled in the figure.
In the text

	Fig. 2. Relative parallax error as a function of the GaiaG magnitude. The points are colour-coded according to their iron abundance.
In the text

	Fig. 3. Reliability of parallaxes based on the two parameters RUWE and astrometric_gof_al (in absolute value). For reference, the typical limits at 1.4 and 12.5 have been reported with dashed lines.
In the text

	Fig. 4. Examples of corner plot with the posterior distributions for the output parameters and the best-fit solution obtained for the WJKS (top) and WcHVI (bottom) Wesenheit magnitudes. The units of α,  β,  γ,  ϵ and of their uncertainties are mag, mag/dex, mag/dex, and mas, respectively.
In the text

Fig. 5. Comparison between final photometric parallax (i.e. the one calculated with the final result of the fitting procedure) and that from Gaia. In filled circles and filled stars are the objects not considered or considered as outliers by the Cauchy loss function. All the data are colour-coded according to their [Fe/H] ratio. The top and bottom panels show the results for the WJKS and WcHVI Wesenheit magnitudes, respectively. The units of α,  β,  γ,  ϵ and ϖ are mag, mag/dex, mag/dex, mas, and mas, respectively.

In the text

Fig. 6. Values of the four parameters, α,  β,  γ, and ϵ, for the different Wesenheit magnitudes considered in this work. Red dots show the results obtained by adopting the Fitzpatrick (1999) reddening law (F99 in the label). The shaded area comprises the range −0.4 < γ < −0.2 mag/dex, where most of the recent literature results are located. The units of α,  β,  γ,  ϵ and of their uncertainties are mag, mag/dex, mag/dex, and mas, respectively.

In the text

	Fig. 7. LMC distance moduli calculated from our best-fitting PWZ relations. Red dots show results relative to the Wesenheit definition according to the Fitzpatrick (1999) reddening law. The shaded region corresponds to the geometric distance of the LMC by Pietrzyński et al. (2019) ±1σ.
In the text

Fig. 8. Comparison between the PWZ parameters obtained with our baseline sample (F+1O pulsators, using [Fe/H], black squares) and those obtained: (i) using a pure σ-clipping algorithm (red circles); (ii) using only F-mode pulsators with [Fe/H] as abundance indicator (blue triangles); (iii) using F+1O pulsators with the total metallicity [M/H] as abundance indicator (green diamonds); and (iv) using only objects with G < 12.5 mag (violet stars). The shaded area comprises the range −0.4 < γ < −0.2 mag/dex, where most of the recent literature results are located. The units of α,  β,  γ,  ϵ and of their uncertainties are mag, mag/dex, mag/dex, and mas, respectively.

In the text

Fig. 9. Estimate of the γ parameter for the WcHVI Wesenheit magnitude obtained by binning our data in three [Fe/H] intervals (filled blue circles). The top and middle panels show the results for two different choices of the β parameter (see labels). The bottom panel shows the results for the same slope as in the middle panel but for a different value of ϵ. For comparison, filled magenta circles show the results from Breuval et al. (2024).

In the text

	Fig. B.1. Comparison between the F106W, F555W, and F814W magnitudes observed with HST by Riess et al. (2021) and those calculated based on ground-based H, V, I magnitudes transformed into the corresponding HST magnitudes using the conversions by Riess et al. (2021).
In the text

	Fig. B.2. As in Fig. B.1 but for the WVI and WHVI Wesenheit magnitudes.
In the text

	Fig. C.1. Comparison of the period-[Fe/H] distributions between the samples adopted in C-MetaLL IV paper (Trentin et al. 2024b) (blue density distribution) and in the present work (red filled circles).
In the text

	Fig. C.2. As in Fig. 2, but for the sample adopted in C-MetaLL IV paper (Trentin et al. 2024b).).
In the text

	Fig. D.1. Corner plot with the posterior distributions for the output parameters and the best-fit solution obtained for the Riess et al. (2021) sample and all parameters free to vary.
In the text