Cepheid Metallicity in the Leavitt Law (C- MetaLL) survey: IV. The metallicity dependence of Cepheid Period-Luminosity relations

Classical Cepheids (DCEPs) play a fundamental role in the calibration of the extra-galactic distance ladder which eventually leads to the determination of the Hubble constant($H_0$) thanks to the period-luminosity ($PL$) and period-Wesenheit ($PW$) relations exhibited by these pulsating variables. Therefore, it is of great importance to establish the dependence of $PL/PW$ relations on metallicity. We aim at quantifying the metallicity dependence of the Galactic DCEPs' $PL/PW$ relations for a variety of photometric bands ranging from optical to near-infrared. We gathered a literature sample of 910 DCEPs with available [Fe/H] values from high-resolution spectroscopy or metallicities from \gaia\ Radial Velocity Spectrometer. For all these stars, we collected photometry in the $G_{BP},G_{RP},G,I,V,J,H,K_S$ bands and astrometry from the \gaia\ DR3. These data have been used to investigate the metal dependence of both intercepts and slopes of a variety of $PL/PW$ relations at multiple wavelengths. We find a large negative metallicity effect on the intercept ($\gamma$ coefficient) of all the $PL/PW$ relations investigated in this work, while present data still do not allow us to draw firm conclusions regarding the metal dependence of the slope ($\delta$ coefficient). The typical values of $\gamma$ are around $-0.4:-0.5$ mag/dex, i.e. larger than most of the recent determinations present in the literature. We carried out several tests which confirm the robustness of our results. As in our previous works, we find that the inclusion of global zero point offset of \gaia\ parallaxes provides smaller values of $\gamma$ (in an absolute sense). However, the assumption of the geometric distance of the LMC seems to indicate that larger values of $\gamma$ (in an absolute sense) would be preferred.


Introduction
Classical Cepheids (DCEPs) are the most important standard candles of the extra-galactic distance scale thanks to the Leavitt Law (Leavitt & Pickering 1912), which is a relationship between period and luminosity (PL) of DCEPs.Once calibrated using independent distances based on geometric methods such as trigonometric parallaxes, eclipsing binaries and water masers, these relations constitute the first step for the cosmic distance scale, as they calibrate secondary distance indicators.These include Type Ia Supernovae (SNIa), which in turn allow us to measure the distances of distant galaxies, located in the steady Hubble flow.The calibration of this three-step procedure (also called cosmic distance ladder) allows us to reach the Hubble flow where the constant (the Hubble constant H 0 ) which connects 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 0 is an important quantity in cosmology because it sets the dimension and the age of the universe.Therefore, measuring the value ⋆ E-mail: etrentin@aip.de of the constant with an accuracy of 1% is one of the most important quests of modern astrophysics.However, there is currently a well-known discrepancy between the values of H 0 obtained by the SH0ES 1 project through the cosmic distance ladder (H 0 =73.01±0.99km s −1 Mpc −1 , Riess et al. 2022) and the ones measured by the Planck Cosmic Microwave Background (CMB) project using the flat Λ Cold Dark Matter (ΛCDM) model (H 0 =67.4±0.5 km s −1 Mpc −1 , Planck Collaboration et al. 2020).This 5σ discrepancy has so far no solution and, if confirmed, could point towards the need for a revision of the ΛCDM model.
In this context, one of the residual sources of uncertainty in the cosmic distance ladder is represented by the debated metallicity dependence of DCEP PL relations used to calibrate the secondary distance indicators.Indeed, a metallicity variation is predicted to affect the shape and width of the DCEP instability strip (e.g.Caputo et al. 2000) that in turn affects the coefficient of the PL relations (Marconi et al. 2005(Marconi et al. , 2010;;De Somma et al. 2022, and references therein).The dependence on metallicity of the PL relations and of the reddening-free Wesenheit  From top to bottom, the different panels show the comparison for the low, medium and top quality Gaia data, respectively (see text for details).In each panel, the average difference ∆ and its dispersion is displayed.
magnitudes2 even if small when involving near infra-red colours (NIR, see e.g.Fiorentino et al. 2013;Gieren et al. 2018), must be taken into account to avoid systematic effects in the calibration of the extra-galactic distance scale (e.g.Romaniello et al. 2008;Bono et al. 2010;Riess et al. 2016).Direct empirical evaluations of the metallicity dependence of PL relations using Galactic DCEPs with sound [Fe/H] measurements based on high-resolution (HiRes hereafter) spectroscopy have been hampered so far by the lack of accurate independent distances for a significant number of Milky Way (MW) DCEPs.The advent of the Gaia mission (Gaia Collaboration et al. 2016) has completely changed this scenario, providing accurate parallaxes starting from the data release 2 (DR2 Gaia Collaboration et al. 2018), which were further improved with the early data release 3 (Gaia Collaboration et al. 2021).In addition, the Gaia mission secured the discovery of hundreds of new Galactic DCEPs (Clementini et al. 2019;Ripepi et al. 2019Ripepi et al. , 2022) ) which, together with those discovered by other surveys such as the OGLE Galactic Disk survey (Udalski et al. 2018) and the Zwicky Transient Facility (ZTF Chen et al. 2020), constitute a formidable sample of DCEPs useful not only in the context of the extra-galactic distance scale but also for Galactic studies (e.g.Lemasle et al. 2022;Trentin et al. 2023, and references therein).However, until a few years ago the number of DCEPs with measures of metallicity from HiRes spectroscopy was mainly restricted to the solar neighbourhood where the DCEPs span a limited range in [Fe/H], centred on solar or slightly super-solar values, with a small dispersion of 0.2-0.3dex (e.g.Genovali et al. 2014;Luck 2018;Groenewegen 2018;Ripepi et al. 2019).This makes it very difficult to measure with sound statistical significance the metallicity dependence of PL relations using Galactic Cepheids.
In this context, a few years ago we started a project named C-MetaLL (Cepheid -Metallicity in the Leavitt Law, see Ripepi et al. 2021a, for a full description), with the specific goal to measure the chemical abundance of a sample of 250-300 Galactic DCEPs through HiRes spectroscopy, expressly aiming at enlarging the iron abundance range towards the metal-poor regime, i.e.
[Fe/H]< −0.4 dex, where only a few stars have abundance measurements in the literature.In the first two papers of the series (Ripepi et al. 2021a;Trentin et al. 2023), we published accurate abundances for more than 25 chemical species for a total of 114 DCEPs.In particular, in Trentin et al. (2023), we obtained measures for 43 objects with [Fe/H]< −0.4 dex, reaching abundances as metal-poor as −1.1 dex.The scope of this paper is to study the metallicity dependence of the PL relations using a [Fe/H] range of more than 1 dex.In this way, we aim to be able to discern between the two scenarios that came out in the recent literature

Bands
Notes.The Wesenhit magnitude names are contained in column 1, while the ξ values, from the Cardelli law (Cardelli et al. 1989), are in column 2. Note that for brevity, in the following the Wesenheit in the HST bands will be referred to as W HS T H,V−I .
concerning the metallicity dependence of PL relations.In fact, in our previous works using Galactic DCEPs and Gaia parallaxes, we have found a rather large dependence in the NIR bands, on the order of ∼ −0.4 mag/dex (Ripepi et al. 2020(Ripepi et al. , 2021a)), or even more when using the Gaia bands (∼ −0.5 mag/dex  Breuval et al. (2021Breuval et al. ( , 2022)), which obtained a similar result using three Cepheid samples in the MW, Large and Small Magellanic Clouds (LMC and SMC, respectively) as three representative objects with different mean metallicities.
The paper is organized as follows: in Sect. 2 we describe the sample of DCEPs and their properties.In Sect. 3 the method to derive the period-luminosity-metallicity (PLZ) and period-Wesenheit-metallicity (PWZ) relations is discussed; in Sect. 4 and Sect. 5 we describe and discuss our results and in Sect.6 we close the paper with our conclusions.

Description of the data used in this work
In this section, we describe the sample of DCEPs used in this paper and their photometric, spectroscopic and astrometric properties.All the data employed in our analysis are listed in Table 1.

Photometry
Optical V, I 3 photometry is available in the literature for about 488 and 364 DCEPs of our sample, respectively (Groenewegen 2018;Ripepi et al. 2021a).For the remaining stars, we decided to use the homogeneous and precise Gaia G, G BP , G RP photometry transformed into the Johnson-Cousins V, I bands by means of the relations published by Pancino et al. (2022).We calculated these aforementioned magnitudes for all the stars having full G, G BP , G RP values (all DCEPs except one), using the intensityaveraged magnitudes from the Gaia Vizier catalogue I/358/vcep (Ripepi et al. 2022) or the simple average magnitudes from the Gaia source catalogue (Vizier I/355/gaiadr3) for the few stars not present in the quoted catalogue.Figure 1 shows the comparison between the literature V, I magnitudes and those calculated from the Gaia photometry for the stars having both values.While the transformed V photometry appears perfectly compatible with that from the literature (V Lit − V Gaia ∼ 0.0 with dispersion σ ∼ 0.04 mag), the transformed I photometry appears to be slightly too faint: I Lit − I Gaia ∼ 0.035 (mag) with dispersion σ ∼ 0.035 mag.Therefore we used the transformed V bands with 3 The I photometry is in the Cousins system no modification, while we corrected the transformed I bands by increasing their value by 0.035 mag.As for the uncertainties, we assumed 0.02 mag in both V and I for the literature sample when the number is not available in the original publication, while for the remaining stars, we propagated the errors, also taking into account (summing in quadrature) the uncertainties in the transformations provided by Pancino et al. (2022).For a handful of stars, the DCEPs (G BP -G RP ) colours were beyond the validity limit of Pancino et al. (2022) relations.For these stars, we adopted the synthetic V,I magnitudes calculated by Gaia Collaboration et al. (2023) based on Gaia DR3 photometry.
Near-infrared J, H, K S band photometry has been taken from van Leeuwen et al. (2007); Gaia Collaboration et al. (2017); Groenewegen (2018); Ripepi et al. (2021a) for the literature sample, and derived from single epoch 2MASS photometry (Skrutskie et al. 2006) for the remaining stars.To this aim, we adopted the same procedure outlined in Ripepi et al. (2021a), using the ephemerides (periods and epochs of maximum light) from Gaia Vizier catalogue I/358/vcep.As in our previous work, the uncertainties on the mean magnitudes were calculated using Monte Carlo simulations, varying the 2MASS magnitude and the phase of the single epoch photometry within their errors, where the last quantity was calculated using the errors on the periods given in the Gaia catalogue.
In this work, we also considered the Wesenheit indices, which are reddening-free quantities by construction (Brodie & Madore 1980).They are obtained by combining the standard magnitude, in a given photometric band, with a colour term according to the following equation: with X i indicating the generic band and the coefficient ξ 2,3 1 coincides with the total-to-selective absorption and is obtained by assuming an extinction law.The photometric band combinations, adopted in this work, together with the coefficient of the colour term are listed in Table 2.In order to transform the Johnson-Cousins-2MASS ground-based H, V, I photometry into the HST correspondent F160W, F555W, F814W filters, we considered the photometric transformations by Riess et al. (2021).However, we have to caution the reader that these transformations have been derived for metallicities not too far from the solar one; therefore they are not tested for values of [Fe/H] as low as those displayed by many DCEPs of the sample used in this work.In Sect.5.4, we discuss this point more in detail.
In analogy with photometry, the reddening for the literature sample was taken from van Leeuwen et al. (2007); Gaia Collaboration et al. (2017); Groenewegen (2018); Ripepi et al. (2021a), while for the remaining stars, we used the same period-colour relations used in Ripepi et al. (2021a) involving (V − I) colour (their eq.4), thus obtaining E(V − I) excesses that were converted in the corresponding E(B − V) ones using the relation (Tammann et al. 2003).The uncertainties on these reddening values were calculated by summing in quadrature the rms of eq. 4 in Ripepi et al. (2021a), with the errors on V and I magnitudes.However, to be conservative, for values of E(B − V) > 1.0 mag, we assumed an uncertainty of 10% on the reddening.

Metallicity
Metallicities were taken from various literature sources.A large part of the sample is the same as in Trentin et al. (2023, see their sect. 4.1).More in detail, we considered the large compilation Notes.The ID identifies each different fit to the data; α, β, γ and δ are the coefficients of the PLZ/PWZ relations; µ LMC 0 represents the true distance modulus of the LMC; rms is the root mean square; Band specifies the photometric band considered in the fit; Mode identifies the sample adopted; N dat is the number of DCEPs adopted.number of DCEPs with metallicities from HiRes spectroscopy is 635.Concerning homogeneity, apart from the few stars from GALAH, PASTEL or GC17, the main samples adopted in this work are the G18, K22 and combined C-MetaLL data.As remarked above, the G18 sample is already homogeneous, while the K22 abundances were homogenized with those of the C-MetaLL sample in Trentin et al. (2023).As discussed in Ripepi et al. (2021b), the C-MetaLL data have only two stars in common with G18, namely X Sct and V5567 Sgr.For these two stars, the abundances agree well within 0.5 σ.
In addition to this already large sample, we decided to exploit the recent results released by the Gaia DR3.To this aim, we cross-matched the DCEPs catalogue published by Gaia DR3 (Vizier catalogue I/358/vcep) as amended by (Ripepi et al. 2022) with that of the astrophysical parameters (Vizier catalogue I/355/paramp), retaining only matching stars having a global  5), then we assigned a data quality flag to every object.Assigned flags are 1,2,3 for top, medium and low quality respectively, according to Recio-Blanco et al. (2023, see sect.9).The Gaia abundances were derived using the stacked RVS spectra over all the epochs of observations.The resulting spectra are therefore an average of many spectra over the pulsation cycle.As the DCEPs vary their effective temperature, surface gravity and microturbulent velocity along the pulsation cycle, we can in principle expect an impact on the derived abundances.In addition, the abundances published in the Gaia Astrophysical parameters are the mean metallicities [M/H], which may differ from the [Fe/H] scale used for the HiRes sample.To investigate these potential concerns, we compared the corrected Gaia metallicities with those from the HiRes sample for the 475 stars in common.The result is shown in Fig. 2.There is an overall good agreement for all the Gaia subsamples (low, medium and top quality).In all the cases the average difference between HiRes and Gaia results is very low, i.e. ∆ ∼ −0.03 dex, with a comfortable moderate dispersion of the order of 0.11-0.13dex (worst for the low quality Gaia data), indicating that overall the Gaia spectroscopic results are of comparable quality of the HiRes ones and the data can be used together.This is especially true for [Fe/H] HiRes > −0.3 dex.Below this value, the Gaia sample is dominated by lowquality data and the dispersion becomes significantly larger (top panel).We corrected for the small offset of the Gaia data and decided to use in the following only the top and medium-quality data, to avoid including some unreliable abundance values from the low-quality data.The resulting sample is then composed of 910 DCEPs divided into 282 DCEP_1Os, 22 DCEP_1O2Os, 581 DCEP_Fs and 25 DCEP_F1Os, where DCEP_1O2Os and DCEP_F1Os represent multi-mode pulsators.For these two last cases, we adopted the longest period, i.e. the 1O period for the 1O2O pulsators and the F one for the F1O DCEPs.Therefore, our sample includes the equivalent of 304 1O and 606 F mode DCEPs.

Astrometry
To carry out our analysis, we adopted the parallaxes from Gaia DR3, which were individually corrected for the zero point offset adopting the recipe by Lindegren et al. (2021, L21) 4 .Note the 4 https://www.cosmos.esa.int/web/gaia/edr3-codecharacteristic shape of the correction, with a sharp hump around G ∼ 12 mag whose explanation can be found in Lindegren et al. (2021).The individual corrections are shown in Fig. 3 as a function of the magnitude and of the ecliptic latitude.As criteria of the goodness of the astrometry, we adopted two indicators: i) the fidelity_v2 index as tabulated by Rybizki et al. (2022), retaining only objects with values larger than 0.5 and ii) the RUWE parameter published in Gaia DR3.Even if the Gaia documentation recommends using a threshold below 1.4 for good astrometric solutions, we choose to be less conservative and use a slightly larger threshold equal to 1.5.This choice allows us to retain possibly good objects falling just outside the 1.4 threshold.In any case, our robust outlier removal procedure (see Sect. 3) will take care of possible deviating stars introduced by the adopted slightly enlarged threshold.
The simultaneous use of the fidelity_v2 and RUWE parameters led us to reject 78 objects from our sample.A large fraction of the rejected DCEPs, namely 35, have G ≤ 7 mag.This is not surprising, as it is known that Gaia parallaxes are still uncertain for such bright stars (e.g.Lindegren et al. 2021).

Derivation of the PLZ/PW Z relations
In this section, we describe the procedure adopted to calibrate the PLZ/PWZ relations.The approach is the same as in Ripepi et al. (2020) and Ripepi et al. (2021a).To avoid any bias, the whole DCEPs sample was considered, including negative parallaxes and without any selection on the parallax relative errors.Moreover, we adopted the astrometry-based-luminosity (ABL) formalism (Feast & Catchpole 1997;Arenou & Luri 1999), which allows to treat the parallax ϖ as a linear parameter: ABL = ϖ10 0.2m−2 = 10 0.2(α+(β+δ[Fe/H])(log P−log P 0 )+γ[Fe/H]) (2) where the parallax ϖ, the apparent generic magnitude m, the period P and the metallicity [Fe/H] are the observables, while the unknowns are the four parameters α, β, γ, δ (intercept, slope, metallicity dependence of the intercept and the slope, respectively).The P 0 quantity is a pivoting period (10d) adopted to reduce the correlation between the α and β parameters which dominate the PL/PW relations.
To take into account the presence of possible outlier measurements, we applied multiple sigma-clipping removals but limited the number of rejected data to ∼ 10% of the full sample.Based on the results described in the previous papers of this series, we decided to exclude the case with no metallicity dependence at all (i.e.γ 0) in Eq. 2. Finally, the 1O pulsators are included in the fitting sample by fundamentalising their periods according to the relation by Alcock et al. (1995) and Feast & Catchpole (1997).In the following sections, we will consider two data samples: i) Lit+DR3 sample including all the selected sources introduced in Sect.2; ii) Lit sample obtained by excluding the sources having only Gaia DR3 metallicity estimate, therefore retaining only DCEPs with metallicity measures from ground-based HiRes spectrographs.
The uncertainties on the parameters calculated in the following analysis are obtained through the bootstrap technique.A set of 1000 resampling experiments is performed and the coefficients of the Eq. 2 are calculated for each obtained random data set.Finally, the quoted errors are estimated by considering the robust standard deviation of the obtained coefficient distributions.

Literature and Gaia DR3 sample
The results of the fitting procedure obtained for this case are listed in Table 3.We studied the dependence of the fitted PLZ/PWZ coefficients on the central wavelength of the considered filters.To associate a characteristic wavelength to the Wesenheit magnitudes, we decided to consider the central wavelength of the main band (e.g.G band for the W G,BP−RP case).The results of this analysis are shown in Fig. 4 for the F+1O sample (i.e. the sample including both F and fundamentalised 1O pulsators) and assuming a metallicity dependence also in the slope of the PLZ/PWZ relation (i.e.δ 0).Looking at the PLZ coefficients (blue filled circles), a strong linear dependence as a function of the λ −1 is evident only for the α and β coefficients (correlation coefficient R 2 ≃ 1), while for γ and δ the slope of the fit is consistent with zero.The relations between α , β and λ −1 are the following: α = (0.618 ± 0.003) • λ −1 + (−3.428 ± 0.006) with rms = 0.019, R 2 = 0.99, and: with rms = 0.059, R 2 = 0.99.The strong dependence of the slope on the wavelength is a well-known feature, as widely discussed by Madore & Freedman (2011).
Comparing these results with those in the previous case, we found good agreement between all the parameters, except for γ in the W HS T H,V−I case, whose values have slightly decreased but are still in agreement within 2σ.To easily compare our γ values with those in the recent literature, we traced two reference lines in the third panel of Fig. 4 and Fig. 5 corresponding to γ = −0.2 and −0.4 dex.Indeed, this range of values includes almost all the recent estimates for the value of γ (see table 1 of Breuval et al. 2022, and Sect. 5.5 for a comprehensive list of recent results).Concerning the δ parameter, the bottom panel of Fig .shows a significant metallicity dependence of the slopes only in the cases of the V-band PL and W HS T H,V−I PW.We note that in both cases, the corresponding γ coefficients are well consistent with the average literature values.
For the sake of completeness, in Appendix A we also plot the results obtained by fitting the Eq. 2 to the sample including only the F pulsators (see Fig. A.1 and Fig. A.2).The exclusion of the fundamentalised 1O pulsators increases the errors on all the coefficients, especially for δ.This parameter assumes also in general higher values and is more scattered compared with the F+1O case.In more detail, we notice that the α and the β coefficients result, within the errors, slightly increased and decreased respectively, while the γ values are smaller (in an absolute sense), except for cases of the G, G bp and W HS T H,V−I magnitudes.The δ values are sensitively higher as well (V band and W HS T H,V−I are again exceptions), but the large uncertainties do not allow us to draw firm conclusions.
To help the visualisation of the fitted relations listed in Table 3, we considered δ = 0, because in this case, the dependence on log(P) in Eq. 2 can be separated from that on metallicity.Figure 6 and Fig. 7 show the projections of the PLZ relations in two dimensions, using P and [Fe/H] on the x-axis respectively.Similarly, Fig. 8 and 9 display the same projections for the PWZ relations.Figures 6 and 8 show no colour trend with the metallicity, indicating that the correction for this parameter was effective.On the other hand, Fig. 7 and 9 display in a direct way the metallicity dependence of the intercepts of the PL and PW relations, respectively.In these figures, the slopes of the solid lines, representing the fits for the different magnitudes, are a direct visualisation of the values of γ.It is possible to appreciate the relevance and the importance of the extension towards the metal-poor regime for the determination of this parameter.Indeed, even a small number of objects with [Fe/H]< −0.4 dex can affect and constrain significantly the slope in these plots.

Literature sample only
The results of the fitting procedure obtained including the literature sample only are listed in Table 4.As for the previous case, we fitted α and β coefficients as a function of λ −1 of the specific band (only for PLZ).These fits are shown in the top two panels of Fig. 10 and their equations are: with rms = 0.04, R 2 = 0.99, and: with rms = 0.06, R 2 = 0.99.The fits shown in Fig. 5 (δ = 0 case for the Lit.sample) for α and β are, respectively: with rms = 0.035, R 2 = 0.99, and: with rms = 0.057, R 2 = 0.99.The inclusion of the δ parameter in the fit of the Lit.sample only, determines a general decrease (in an absolute sense) of the γ values, which are closer to the -0.2,-0.4mag/dex range typical of previous works.However, at the same time, the δ parameters show a considerable dependence of the PLZ/PWZ slopes on metallicity.If instead the δ parameter is neglected, as shown in Fig. 11, we find again an increase of the absolute value of γ, as already found for the Lit.+Gaia case.Figures A.3 and A.4 in Appendix A report the above analysis considering only the F pulsators.A comparison of the results obtained for the F+1O and the F samples reveals no significant differences.

Comparison between Lit.+Gaia and Lit. samples
It is interesting to study the impact of the inclusion of the Gaia DCEP sample together with the HiRes one.To this aim, we compared the estimated parameters for both cases in Fig. 12 and Fig. 13 for the PLZ and PWZ relations, respectively.While the α and β parameters seem to be quite robust, with the exception of the W HS T H,V−I case, for γ and δ the Lit.sample presents slightly smaller and higher values (in an absolute sense) compared with the Lit.+Gaia data set, respectively.In almost all the cases, however, the differences are insignificant within 1σ.For this reason, we expect to find similar trends in the following discussions.Therefore, unless specified, we will refer to the Lit.+Gaia sample only.The reader can find in App.B all the analogues of the figures shown.

Parallax correction
A well-known feature of Gaia parallaxes in DR3 is the need for a global offset in the parallaxes after the L21 correction (e.g.Riess et al. 2021Riess et al. , 2022;;Molinaro et al. 2023, and references therein).The exact amount of this quantity is debated and appears to depend on the properties (e.g.location on the sky, magnitude and colour) and kind of stellar tracer adopted for its estimate (see discussion in Molinaro et al. 2023).To take into account this feature, we also conducted the PLZ/PWZ fit after adding the global offset to the EDR3 parallax values of −14 µas (Riess et al. 2021) or −22 µas (Molinaro et al. 2023).A graphical comparison between these different cases for the Lit.+Gaia sample is shown in Fig. 14 and 15 for the PLZ and PWZ, respectively.The introduction of the parallax correction increases the absolute values of α and β by 1-4%, while γ decreases monotonically (in an absolute sense) for larger values of the parallax offset.This is particularly visible for the PWZ relations which show much less scatter than the PLZ ones.In any case, the maximum variation is barely larger than 1σ.The δ coefficient varies less than γ with a slightly decreasing trend which is however insignificant at the 1σ level.Overall, these effects have already been found and are coherent with those discussed in Ripepi et al. (2021a).In practice, the introduction of a significant global parallax zero point offset tends to reconcile our metallicity dependence with the typical literature values, as γ tends to diminish at increasing offset (both in an absolute sense).But, as discussed in the next section, the adoption of large offset values affects the absolute distance scale.

The LMC distance
To test the goodness of our PLZ/PWZ relations and to understand which parameter affects the calibration, we used our relations to estimate the distance to the LMC and compared the results with the currently accepted geometric value µ 0 = 18.477 ± 0.026 mag (Pietrzyński et al. 2019).To this aim, we used the same method described in detail in Ripepi et al. (2021a).The distance moduli µ LMC 0 , obtained in this way are listed in Table 3 and 4 for the Lit.+Gaia and Lit.sample, respectively.The samples including only F pulsators generally give worse results than the F+1O sample, therefore, in the following discussion we use only this last data set.Figure 16 shows the resulting µ LMC 0 for the PLZ/PWZ relations for Lit.+Gaia samples.The introduction of the global parallax correction worsens the distance estimation with the exception of the W HS T H,V−I .In more detail, the results with no correction and with Riess et al. (2021)'s correction embrace approximately the lower and upper allowed limit respectively, suggesting an intermediate value.Therefore, if we use the Pietrzyński et al. (2019) distance of the LMC as a reference, we have to conclude that the rather large values of γ found in this work are favoured compared to the smaller ones from the recent literature.

Metallicity sampling
Even if the results of our C-MetaLL project enlarged considerably the number of DCEPs having [Fe/H] < −0.5 dex, these objects represent less than 10% of the pulsators used in this work, while the bulk of the sources is distributed around the solar value.To test whether this unbalanced metallicity distribution could affect our determination of the PLZ/PWZ relations, we have conducted an experiment whose goal was to recalculate these relationships by resampling the input data such that its metallicity distribution is uniform over the whole range.To this aim, we divided the sample into five bins in metallicity imposing the same number of stars in each bin.As there are few stars in the low metallicity regime, the sources with [Fe/H] < −0.5 dex  were all kept in the same bin and their number, N [FeH<−0.5]= 44 (considering the astrometric selection described in Sect.2.3) set the number of pulsators included in each of the four remaining bins, whose size is approximately 0.25 dex in metallicity.The DCEPs populating each bin, apart from the lowest metallicity one, were randomly picked.Having resampled the data, we car-ried out the fit as in Sect. 3 considering the Gaia+Lit.data set only, given that the other cases provide similar results.We also restricted the test to the F+1O sample, because the number of F pulsators in the most metal-poor bin is less than half of the F+1O combined sample, meaning that the total number of DCEPs used for the fit would be too small to achieve meaningful results.We repeated 10,000 times the fitting procedure, obtaining a distribution of values for the α, β, γ, δ coefficients.Their median values and relative uncertainties (scaled median absolute deviation or MAD) are compared with the results obtained using the entire sample in Fig. 17 and 18 for PLZ and PWZ, respectively.Both figures show an excellent agreement for all the coefficients.Therefore, we can conclude that the results obtained in this work are not affected by the sample's unbalanced distribution in the metallicity.

The case of the W HS T H,V−I
The PW relations involving the W HS T H,V−I magnitude presented in the previous sections show a peculiar behaviour compared with the other bands.We have already mentioned in Sect.2.1 that a possible explanation for this occurrence is that the photometric transformations between the Johnson-Cousins H, V, I bands could be inadequate at low metallicities, given that they were calculated for pulsators with abundances close to the solar one.To test this working hypothesis, we re-fitted the period-W HS T H,V−I relation by using a metallicity range similar to that used to calculate the photometric transformations.In particular, we adopted the metallicity interval −0.21 <[Fe/H]< 0.39 dex covered by the DCEPs having HST photometry (see Riess et al. 2022), thus excluding the metal-poor tail of the distribution.Because of the narrower range in metallicity, the number of useful pulsators is reduced significantly and we decided to consider only the Lit.+Gaia case, with fundamentalized 1O sources, in order to have a fitted sample as large as possible.Finally, the fitted sample reduces to 686 pulsators, and we obtain α = −5.430± 0.017, β = −2.977± 0.040, γ = −0.142± 0.096, neglecting the δ term, and µ LMC 0 = 18.310 ± 0.027 mag.Whilst the former two parameters agree with the results obtained using the entire sample (see Table 3), the γ term and the LMC distance diverge by more than 1σ.This suggests that the transformation between Johnson-Cousins and HST magnitudes, not tested for low metallicities values, might have an impact on the final results for this band, although we cannot exclude that this difference is due to the cut made in metallicity, spanning now a range less than 1 dex.

Comparison with the literature
Figure 19 shows several empirical estimations of γ throughout the last 20 years.Results from this work are also plotted, considering the case of the Lit.+Gaia data sample.In more detail, different techniques have been used in order to analyse the metallicity effect on PL and PW relations.Early studies adopting distances from the Baade-Wesselink (BW) analysis (e.g.Storm et al. 2004Storm et al. , 2011;;Groenewegen 2013) reported very small values for the γ parameter in the V,I,K bands and in the W JK ,W V K ,W V I Wesenheit magnitudes.More recently, a stronger effect (γ ∼ −0.2 mag/dex), albeit smaller than the findings of this work, have been found in the same bands by Gieren et al. (2018) from the BW analysis of DCEPs in the Galaxy, LMC and SMC, which were used to extend the range of metallicity of the pulsators adopted in the analysis.The advent of Gaia DR2 parallaxes permitted us to obtain the first reliable evaluation of the γ term using only Galactic DCEPs having metallicities from HiRes spectroscopy.In particular, Groenewegen (2018) and later Ripepi et al. (2019Ripepi et al. ( , 2020) ) found γ ∼ −0.1 − 0.4 mag/dex in a variety of bands and Wesenheit magnitudes (see Fig. 19), therefore closer to the values found in this paper, but with a significance generally lower than 1σ, owing to the still too poor precision of DR2 parallaxes.
The improved Gaia EDR3 parallaxes, instead, allowed us to obtain larger (in an absolute sense) γ value in the first paper of the C-MetaLL project (Ripepi et al. 2021a) and for the W G magnitude (Ripepi et al. 2022).
Other studies (Groenewegen et al. 2004;Fouqué et al. 2007;Wielgórski et al. 2017;Breuval et al. 2021Breuval et al. , 2022) ) together with the already quoted Gieren et al. (2018), compared the properties (basically the zero points of the PL or PW relations) of the DCEPs in the MW and in the more metal-poor galaxies LMC and SMC to estimate the extent of the γ value.In particular, Breuval et al. (2021Breuval et al. ( , 2022)), used geometric distances for the DCEPs in the MW (from Gaia EDR3 parallaxes) and in the Magellanic Clouds (from eclipsing binaries) to estimate γ in the same bands/Wesenheit magnitudes treated in this work, after fixing the slope of the PL and PW relations to that of the LMC, with results ranging between −0.178 and −0.462 mag/dex.In the context of the SH0ES project (Riess et al. 2016(Riess et al. , 2019(Riess et al. , 2021(Riess et al. , 2022)), the metallicity effect is one of the outputs of the process for the estimation of H 0 .In these cases, the γ coefficient in the W HS T H,V−I magnitude used by the SH0ES team, is of the order of −0.2 mag/dex, while we obtain a slightly larger (in an absolute sense) value in our best case, i.e. with the Lit.+Gaia sample without the δ term (see e.g.Fig. 5).Recent theoretical studies (Anderson et al. 2016;De Somma et al. 2022) also predict mild effects of the metallicity, with γ ranging from −0.27 to −0.13 mag/dex, albeit the models by De Somma et al. ( 2022) only deals with PW relations.We notice that in almost all the cases we find larger γ values (in an absolute sense) than the literature, albeit most of them are distributed along the literature's lower limit of −0.4 mag/dex (see also Fig. 4).On the other side, for the V-band and the W HS T H,V−I magnitude, smaller coefficients are found, in agreement with the upper limit of −0.2 mag/dex.

Conclusion
In this fourth paper of the C-MetaLL series, we studied the metallicity dependence of the PL relations in the following bands: G BP , G RP , G, I, V, J, H, K S and of the PW relations in the following Wesenheit magnitudes: W G,BP−RP , W H,V−I , W HS T H,V−I ,W J,J−K , W V,V−K .To this end, we exploited the literature to compile a sample of 910 DCEPs having [Fe/H] measured from HiRes spectroscopy and complemented it with a number of stars having metallicity measurements based on the Gaia RVS instrument as released in DR3.For all these stars, we provide a table with photometry in the G BP , G RP , G, I, V, J, H, K S bands, metallicity and astrometry from Gaia DR3.
We carried out our analysis adopting two different samples: one composed only of literature data and one including both the latter and Gaia DR3 results.In order to estimate the parameters of the PLZ/PWZ relations, we used the ABL formalism which allowed us to treat the parallax ϖ linearly and to preserve the statistical properties of its uncertainties.We considered two functional forms for the PLZ/PWZ relations including i) the metallicity dependence on the intercept only (three-parameter case); ii) the metallicity dependence on both intercept and slope (fourparameter case).The main findings of our analysis are: 1. Regarding PLZ relations, both α and β show a linear dependence on the wavelength.We provide the linear relation-ships between these quantities and λ −1 .These relations do not change considerably either when we use the Lit.+Gaia or Lit.samples or if we carry out a three or four-parameter fit.Because of the low wavelength coverage, no clear dependence can be discussed for the PWZ relations; 2. A clear negative dependence of the intercept on the metallicity (γ-coefficient) is found for all the PL and PW relations.For the three-parameter solutions, the values of γ are around −0.4 : −0.5 dex with no clear dependence on the wavelength.Thus, this work confirms our previous results reported by Ripepi et al. (2021a); Molinaro et al. (2023).
For the four-parameter solutions, the values of γ generally slightly decrease (in an absolute sense), especially for the Lit.only sample.In general, the γ coefficients found in this work are larger (in an absolute sense) than those in the literature which range between −0.2 and −0.4 dex; 3. The dependence of the slope on the metallicity (δ coefficient) remains undetermined.Indeed, in about half of the cases, δ assumes values comparable with 0 within 1σ for the Gaia+Lit.sample, while, especially in the Lit.case it generally assumes a positive value comprised between 0 and +0.5.It is worth noticing that the W HS T H,V−I usually shows a different behaviour compared with all the other magnitudes.A possible explanation for this occurrence is that, as stated in Sect.2.1, the transformations between Johnson-Cousins and HST magnitudes are not tested at low metallicity.4. The main difference between using the F and F+1O samples consists of larger error bars for the former case, especially for the γ and δ coefficients; 5.The inclusion of global zero point offset to the individually corrected parallaxes according to L21 has a larger impact on the γ coefficients rather than on the δ ones.More in detail, the adoption of offsets by −14 µas (Riess et al. 2021) and −22 µas (Molinaro et al. 2023) implicates smaller and smaller values of γ (in an absolute sense), therefore more in agreement with recent literature.However, if we use the geometric distance of the LMC by Pietrzyński et al. (2019) as a reference and calculate the distance of this galaxy using our relations with and without the global offset, we find that a good agreement for the distance of the LMC is found for values of the offset in between null and the 14 µas values.These results support larger values (in an absolute sense) of the γ value.6.The possible effect of the uneven distribution in the metallicity of the sample used in this work was investigated by resampling our data set in order to have a balanced number of DCEPs at every metallicity value.We have carried out the fitting procedure on 10,000 samples extracted from the total sample obtaining a value for the four coefficients of the fit for every data set.The medians of the obtained distributions for all the coefficients appear in excellent agreement with those obtained using the entire sample.Therefore, we can conclude that the results obtained in this work are not affected by the sample's unbalanced distribution in the metallicity.
The results presented in this paper show the importance of the extension of the metallicity range when carrying out the analysis.Further observations with the aim of gathering HiRes spectroscopy of MW metal-poor DCEPs will offer the opportunity to better constrain the dependence on the metallicity of both the intercept and the slope.In particular, regions in the Galactic anticentre direction are the most promising targets of future observations, where DCEPs are expected to have [Fe/H]< −0.3 : −0.4 dex, thus further populating the metal-poor tail of DCEPs distribution.

Fig. 1 .
Fig. 1.Comparison between the literature V, I photometry and that from Gaia through the transformation byPancino et al. (2022).The top and bottom panels show the comparison in V and I bands, respectively.

Fig. 2 .
Fig. 2. Comparison between the [Fe/H] from HiRes data and Gaia [M/H].From top to bottom, the different panels show the comparison for the low, medium and top quality Gaia data, respectively (see text for details).In each panel, the average difference ∆ and its dispersion is displayed.

Fig. 3 .
Fig. 3. Individual parallax corrections from L21 as a function of the G magnitude.The points are colour coded according to the ecliptic latitude.

Fig. 4 .
Fig. 4. Coefficients of PLZ/PWZ relations, obtained from the fit to the Lit.+ Gaia sample, plotted against the λ −1 parameter.From top to bottom, the shown panels correspond to the α, β, γ and δ coefficients.In each panel, PLZ and PWZ coefficients are plotted with blue and red symbols respectively, while the solid line shows the linear fit only for the PLZ relations.Grey dashed lines in the γ panel delimit the range of results in the literature (−0.2 to −0.4 mag/dex, seeBreuval et al. 2022,  and Sect.5.5), while in the δ panel, it corresponds to the value 0.

Fig. 6 .
Fig. 6.PL relations for the bands studied in the paper where the magnitudes have been subtracted by the metallicity contribution.The colour bar represents the metallicity.

Fig. 7 .
Fig. 7. Visualization of the PL dependence on the metallicity.The y-axis is the magnitudes subtracted by the period contribution and the x-axis is the metallicity.The colour bar represents the logarithm of the period.

Fig. 12 .Fig. 13 .
Fig. 12.Comparison between the PLZ coefficients, obtained with the Lit.+Gaia data set and those obtained with the Lit.data set.From topleft to bottom-right, the panels show the results for the α, β, γ and δ coefficients respectively.For each coefficient, different photometric bands are plotted using different symbols, as labelled in the top-left panel.Grey-filled and empty symbols indicate the fit results including or excluding the δ parameter, respectively.

Fig. 14 .
Fig. 14.Comparison between the PLZ coefficients, obtained for different global offset values.The coefficient values obtained without offset assumption, are chosen as reference and are on the x-axis, while those obtained by assuming the global parallax offset by Riess et al. (2021) and Molinaro et al. (2023) are on the y-axis and are plotted respectively with green and pink symbols.From top-left to bottom-right, the panels show the results for the α, β, γ and δ coefficients respectively.For each coefficient, different photometric bands are plotted using different symbols, as labelled in the top-left panel, while grey-filled and white-filled symbols indicate the fit results including or excluding the δ parameter, respectively.

Fig. 16 .
Fig. 16.The LMC distance modulus, obtained using the relations of Table 3 for the F+1O cases, is plotted against the central wavelength of the considered photometric bands.To avoid label and symbol overlapping, the I band result has been shifted by -0.1µm −1 .The top panel shows the results for the PLZ relations, while the bottom panel shows those for the PWZ relations.Different colours indicate the results obtained for different parallax shift values: L21 correction (black squares), L21 + Riess et al. (2021) offset (green squares), L21 + Molinaro et al. (2023) (violet squares).Empty and grey-filled squares indicate the results obtained excluding and including the δ coefficient, respectively.

Fig. 17 .
Fig. 17.Comparison between the PLZ coefficients obtained with the resampling in metallicity with those obtained using the entire data set.From top-left to bottom-right, the panels show the results for the α, β, γ and δ coefficients respectively.For each coefficient, different photometric bands are plotted using different symbols, as labelled in the top-left panel.In contrast, red-filled and white-filled symbols indicate respectively the fit results including or excluding the δ parameter.

Fig. A. 3 .
Fig. A.3.Same as Fig. 10 but for the Lit.data set including only F pulsators.

Fig. A. 4 .
Fig. A.4. Same as Fig. 11 but for the Lit.data set including only F pulsators.

Fig. B. 1 .
Fig. B.1.The same as Fig. 14 but for the Lit.sample.

Table 1 .
Photometric, astrometric and spectroscopic data for the DCEP sample used in this work.Only the first ten lines of the table are shown here to guide the reader to its content.The machine-readable version of the table will be published at the CDS (Centre de Données astronomiques de Strasbourg, https://cds.u-strasbg.fr/).

Table 2 .
The photometric bands and the colour coefficients adopted to calculate the Wesenheit magnitudes considered in this work.