© The Authors 2025

1. Introduction

Classical Cepheids (CCs) are the most important primary distance indicators when calibrating the cosmic distance scale. The physical basis of the existence of period–luminosity (PL) and period–luminosity–colour (PLC) relations for these stars is well known. The existence of a period–luminosity–temperature–mass (PTLM) relation directly descends from the period–mean density relation and the Stefan-Boltzmann law (see e.g. Bono et al. 1999; Marconi et al. 2013a, and references therein). Moreover, it is well known from stellar evolution models that intermediate-mass stars, depending on their chemical composition, obey a mass–luminosity (ML) relation, during their central He-burning phase (see e.g. Chiosi et al. 1992, 1993; Bono et al. 1999). According to several authors (see e.g. Chiosi et al. 1993; Bono et al. 1999; Anderson et al. 2014, 2016, and references therein), the presence of non-canonical effects (e.g. core overshooting, mass loss, and rotation) results in a systematic increase in the luminosity obtained from the ‘canonical’ ML relation. For example, De Somma et al. (2020, 2022) considered two over-luminosity levels above the canonical relation (+0.2 dex and +0.4 dex) to take non-canonical effects into account. By combining the PTLM with the ML relation, one obtains aperiod–luminosity–temperature (PLT) or PLC relation. The widely adopted PL relations result from the projection of the PLC relation onto a PL plane (see Madore & Freedman 1991), which implies neglecting the finite colour extension of the instability strip. Another relation, the so-called period–Wesenheit (PW) relation (Madore 1982), is obtained by replacing the magnitude with a function defined as the magnitude corrected for a colour term, multiplied by a coefficient that is exactly the ratio between the total and the selective extinction. This function is reddening-free by definition, provided that an extinction law is assumed a priori. Thus, the coefficients of the ML relation affect the PL, PLC, and PW relations and, in turn, the Cepheid distance scale. In this context, a combination of accurate geometric distance information and independent pulsation modelling results is needed to constrain the ML relation. The third ESA Gaia Data Release (Gaia DR3) includes a sample of 799 Cepheids and 1096 RR Lyrae (see e.g. Katz et al. 2023; Clementini et al. 2023; Ripepi et al. 2023, and references therein) with parallaxes, epoch photometry, and radial velocity (RV) measurements. The adopted dataset is the fundamental basis forinvestigating Cepheid stellar parameters.

From a theoretical point of view, the use of non-linear convective pulsation models to compare predicted and observed light and RV curves represents a powerful tool for simultaneously constraining the intrinsic stellar parameters and the individual distances of the investigated pulsating stars. Several applications in the literature support the predictive capabilities of this method for both the CCs (see e.g. Wood et al. 1997; Bono et al. 2002; Keller & Wood 2002, 2006; Natale et al. 2008; Marconi et al. 2013a,b, 2017; Ragosta et al. 2019) and RR Lyrae (see e.g. Bono et al. 2000a; Marconi & Clementini 2005; Marconi & Degl’innocenti 2007).

These pulsation codes are not only capable of predicting stellar parameters, but they also have been able to describe dynamical phenomena of these stars, for example beat Cepheids (Kolláth et al. 2002), double mode RR Lyrae behaviour (Szabó et al. 2004), and period doubling and related non-linear phenomena (Kolláth et al. 2011; Plachy & Kolláth 2013; Smolec 2014, 2016). Meanwhile, PL relations have been extensively studied as well (see e.g. Marconi et al. 2015, 2018; Deka et al. 2022; Das et al. 2021, 2024)

In this work we used a new extensive set of non-linear convective pulsation models to fit a subsample of 46 fundamental CCs from the Gaia DR3 dataset. This represents the largest sample of CCs analysed so far through the model fitting technique.

The work is structured as follows: In Sect. 2 we present the adopted sample, and the new computed pulsation models are discussed in Sect. 3. The developed model-fitting approach is outlined in Sect. 4, together with the selection of the best-fit models based on the inclusion of Monte Carlo simulations. Section 6 presents the results of the fitting method, and Sect. 7 some closing remarks.

2. The sample selection

The initially adopted sample of sources was obtained from a query to the Gaia public archive¹. All the DR3 fundamental mode CCs, all of which have epoch RV data were extracted and amounted to 799 variables. This catalogue was cross-matched with the list of confirmed Galactic CCs by Ripepi et al. (2023). This cross-match left us with 558 variables. Since the pulsational model calculation is very time-consuming, we decided to reduce the sample by selecting only those having a good phase sampling in the RV curve². To this aim, we modelled the RV time series by fitting them with a truncated Fourier series (Appendix A) and selected only those without large and/or spurious model oscillations, possibly due to wide phase intervals with a dearth of measures. In particular, a selection was performed considering only sources according to the following conditions: (i) RV curves with at least 20 data points; (ii) a uniformity index³UI_RV > 0.9; (iii) visual inspection to remove the remaining few bad sources; and (iv) DR3 relative parallax error σ_ϖ/ϖ < 50%.

After this procedure, the sample was reduced to 198 selected sources. This sample was further reduced by looking at the residuals around the Fourier models, for both photometric and RV data. Specifically, we performed a statistical analysis of the Fourier fit scatter for all time series. In particular, we calculated the root mean square (rms) of residuals around the Fourier model for all the considered photometric and RV curves and removed all the sources with the rms of residuals deviating by more than 3.5 median absolute deviations (MADs) from the median of the rms distribution, for one or more observed time series. The remaining sample after this step contains 46 sources.

We then complemented the information with available metal abundance measurements, as metallicity is a key parameter affecting stellar pulsation and the shape of the light and RV variations. We adopted the work by Trentin et al. (2024, T24), which reports the iron abundance for more than 900 CCs obtained based on high-resolution spectroscopy.

The period and the [Fe/H] distributions of our final sample are shown respectively in Figs. 1 and 2, while Fig. 3 shows the distribution of Z values, obtained⁴ by assuming the elemental composition of the Sun Z_⊙ = 0.02 dex (see e.g. Buldgen et al. 2024, and references therein)

Fig. 1. Distribution of the selected classical Cepheids period values obtained from the DR3.

Fig. 2. Distribution of the selected classical Cepheids [Fe/H] values obtained through high-resolution spectroscopy (T24).

Fig. 3. Distribution of the Z values of all the selected variables computed by converting the spectroscopic [Fe/H] values from Trentin et al. (2024).

3. Pulsational model computations

The theoretical framework considered in this work is based on a non-linear radial pulsation code that adopts a non-local time-dependent treatment of the turbulent convection (see e.g. Stellingwerf 1982; Bono & Stellingwerf 1994; Bono et al. 1999, and references therein).

The procedure followed to obtain the hydrodynamic pulsational models reproducing the observed light and RV curves, consists of four main steps: (i) the definition of a multi-dimensional grid of possible structural parameter (M , T_eff , log(L/L_⊙) , α ⁵) values to be given as input to the hydrodynamical code to reproduce the observed period; (ii) the linear analysis that provides the static envelope structure to be perturbed through the subsequent non-linear analysis; (iii) the non-linear analysis aiming at producing the final full-amplitude pulsating models, including both absolute bolometric magnitude and pulsational velocity curves; and (iv) the transformation of the theoretical bolometric light curves into the Gaia photometric system.

3.1. Initial grid of structural parameters

Since chemical composition plays a key role in the computation of pulsation models, and given that the selected sample spans a wide range in Z, from Z < 0.01 dex to Z ≃ 0.045 , we decided to divide the sample into four metallicity bins centred around the values Z₁ = 0.01 dex, Z₂ = 0.02 dex, Z₃ = 0.03 dex, and Z₄ = 0.04 dex, each characterized by a width of 0.005 dex. The helium contents (Y), corresponding to the adopted Z bins, are Y₁ = 0.26 , Y₂ = 0.28 , Y₃ = 0.28, and Y₄ = 0.29 dex. Based on this division, each bin contains 13, 24, 5, and 4 sources, respectively. From this point forwards, the central value of each bin will be representative of the elemental composition of the stars within that bin.

For each fixed chemical composition, the input mass and effective temperature values for the linear analysis were varied on a rectangular grid, covering the ranges $3<\mathrm{M}/{\mathrm{M}}_{\odot }<13$ with a step of 0.5 M_⊙ and $3500<{\mathrm{T}}_{\mathrm{eff}}(\mathrm{K})<7500$ with a step of 50 K. Then, for each combination of mass and effective temperature, the luminosity value is found by satisfying the observed period, as described in Appendix B. In this procedure only luminosities satisfying a variation within −0.06 dex to 0.6 dex with respect to the canonical ML relation, as predicted by Bono et al. (2000b), are kept for the following analysis (see Appendix B). The lower limit corresponds to –3 times the standard deviation of the adopted ML relation (Bono et al. 2000b), while the upper limit takes into account the presence of full core overshoot (Chiosi et al. 1993). Taking advantage of the well-known relations between luminosity and the other structural parameters, our method allows us to define a complete grid of initial structural parameters, by significantly reducing the number of models to be computed.

As regards to the α parameter, it was varied between 1.4 and 1.9 in steps of 0.05, resulting in 11 possible α values. We note that when the α values do not match those provided in Table 3 of De Somma et al. (2020), we applied linear interpolation if our values fell between two of those considered by the authors, or we adopted their most extreme relations if our values lay outside their specified α range.

3.2. The hydrodynamical analysis

The linear survey is based on a non-adiabatic radiative code (see Bono et al. 1999, and references therein, for detail) that was run for all the M–T_eff –log(L/L_⊙) values defined above for each assumed chemical composition. The subsequent non-linear analysis takes the static structure and the predicted eigenfunctions obtained from the linear code as input and enables the determination of which models exhibit a full amplitude stable pulsational behaviour.

The non-linear code introduces an initial dynamical perturbation (typically 20 km/sec) to the static model, modulated within the envelope structure with the predicted radial eigenfunction obtained from the linear analysis. Subsequently, it monitors the time evolution of its characteristic parameters to assess the attainment of full amplitude behaviour. Integration of hydrodynamical equations is performed over a considerable number of periods, with the satisfaction of convergence conditionsdetermined by the total work approaching vanishing values and the peak-to-peak bolometric amplitude exhibiting negligible change (see below).

Typically, the code is initiated with a fixed number of computing cycles, ranging from 1000 to around 10 000, depending on the model’s mass and position within the instability strip, with models located close to the extreme edges of the strip typically requiring a larger number of cycles. Given the resource-intensive nature of the non-linear step, we implemented runtime tests within the code to assess the model’s convergence status. Specifically, the non-linear code computation is halted under the following circumstances over the last 40 period cycles:

(i) The normalized total work remains stable around vanishing values (≤5 ⋅ 10⁻⁴), and the peak-to-peak bolometric magnitude amplitude shows minimal variation, staying below 5 ⋅ 10⁻³ mag.

(ii) The peak-to-peak amplitude is consistently smaller than 0.02 mag and steadily decreasing.

(iii) The conditions described in point (i) are met, but the period converges to the expected value associated with another pulsation mode (the model reaches a stable pulsational regime but with a mode that is not the required one).

Upon meeting the conditions in point (i), the model achieves full amplitude behaviour, pulsating stably. Conversely, if conditions (ii) or (iii) are satisfied, the pulsation ceases or the model switches to a different pulsational mode. This strategy allowed us to speed up the model computation by stopping it before the assumed total number of cycles is reached, as soon as the convergence condition is achieved, and in particular when the pulsation behaviour switches off or changes mode.

3.3. Converting from bolometric to Gaia multi-filter light curves

The conversion of the bolometric magnitude curves into the Gaia photometric bands was achieved by adopting the bolometric correction grids based on the PHOENIX model atmosphere (Chen et al. 2019). We refer the reader to Marconi et al. (2021) for all the details of our adopted procedure. Here we just recall that we used a proprietary C code to interpolate the Chen et al. (2019) tables to obtain the bolometric corrections corresponding to the log(g) , T_eff, and Z of our modelled curves. For all the transformed models we assumed ${\mathrm{M}}_{\odot }^{\mathrm{bol}}=4.79$ mag (De Somma et al. 2022, and references therein).

4. The model fitting method

The adopted model fitting procedure is based on the minimization of the following χ² functions (see also Marconi et al. 2013a, 2017; Ragosta et al. 2019):

$$\begin{array}{c}\hfill {\chi }_{j=1\cdots N_{\mathit{bands}}}^2={\displaystyle \sum _{i=1}^{N_{\mathit{pts}}^j}}{\left(\frac{m_i^j\left({\varphi }_i^j\right)-M^j({\varphi }_i^j+\delta {\mathrm{\Phi }}^j)+{\mu }^j}{{\sigma }_i^j}\right)}^2+{\chi }_{Am{p}^j}^2\end{array}$$(1)

$$\begin{array}{c}\hfill {\chi }_{\mathit{RV}}^2={\displaystyle \sum _{l=1}^{M_{\mathrm{pts}}}}{\left(\frac{v_l\left({\varphi }_l\right)-a(V_{\mathrm{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})+\gamma )}{{\sigma }_l^{\mathit{RV}}}\right)}^2+{\chi }_{Am{p}_{\mathit{RV}}}^2,\end{array}$$(2)

where we included the term χ_Amp² for every fitted band (including RVs), which is equal to

$$\begin{array}{c}\hfill {\chi }_{Am{p}^j}^2={\left(\frac{Am{p}_{\mathit{Model}}^j-Am{p}_{\mathit{Fourier}}^j}{{\sigma }_{\mathit{Amp}}^j}\right)}^2.\end{array}$$(3)

To justify the inclusion of this additional term, we performed tests that revealed that, in some cases, models minimizing the χ² function without the ${\chi }_{\mathrm{Amp}}^2$ term, do not always reproduce the peak-to-peak amplitude of the observational curves. We therefore decided to give greater weight to models that also reproduce as good as possible the amplitude, whose reference value (Amp_Fourier^j) is derived from the Fourier model constructed as previously described.

The first equation above refers to the j-th photometric band, while the second equation is used to fit the RV curve. The sums are performed over the number of observational data points for both photometry (N_pts) and RV (M_pts). The lower case symbols, ϕ, m, and v indicate the pulsation phases, the measured magnitude, and the measured RV, respectively. The upper case functions M^j(ϕ) and V_puls(ϕ) give analytic expressions for the model j-th photometric band and the pulsational velocity curves⁶. The expected outputs of this best-fitting procedure are the distance modulus, μ^j, and the phase shift, $\delta {\mathrm{\Phi }}^{\mathrm{j}}$, in the j-th band for the photometry, while the model fitting of the RV curves provides information on the projection factor, $(\mathrm{p}=-\frac{1}{\mathrm{a}})$, the velocity of the centre of mass, γ, and the phase shift, $\delta {\mathrm{\Phi }}_{\mathrm{RV}}$.

The projection factor (p-factor) allows us to compare the modelled pulsational velocity to spectroscopically measured RV. It quantifies the projection of the pulsational velocity of the stellar atmosphere along the line of sight. In particular it is defined as p = V_puls/RV and assumes values larger than 1.0. The value of the p-factor is debated in the literature due to the difficulty of its determination. Indeed, it depends on the physical structure of the stellar atmosphere, on the velocity gradient in the atmosphere and on the way the RV is obtained experimentally (see e.g. Nardetto et al. 2011; Molinaro et al. 2011; Ragosta et al. 2019, and references therein). The first works addressing the problem of the p-factor from a theoretical perspective date back to Eddington (1926) and Carroll (1928). These authors, by including the effects of limb darkening and constant-velocity expansion of the stellar atmosphere in their theoretical models, obtained a p-factor value of p = 24/17 = 1.41. This value is close to the maximum limit, 1.5, which is derived when the Cepheid is approximated as a simple uniform disk (Nardetto 2017). A more complex model, which accounts not only for the geometric effect but also for the relationship between the stellar surface’s RV and that derived from spectral lines, was developed for example by Nardetto et al. (2004). This model also includes the effects of limb darkening and atmospheric velocity gradients. Several studies have investigated the existence of a relationship between the p-factor and the pulsation period, but no definitive conclusion has been reached so far (Groenewegen 2007, 2013; Kervella et al. 2017; Nardetto et al. 2009; Neilson et al. 2012; Storm et al. 2011).

In the current work, we forced the p-factor values to be in the range 1.10–1.5 to take current uncertainties in the literature evaluations into account. The upper limit was fixed according to the expected maximum theoretical value, as stated above, while the lower limit was adopted according to the recent results obtained by Trahin et al. (2021). They calibrated the p-factor value by combining the parallax-of-pulsation method, implemented through the SPIPS algorithm (Mérand et al. 2015), with the Gaia Early Data Release 3 (EDR3) parallaxes (see also Anderson et al. 2024; Breitfelder et al. 2016; Trahin et al. 2021). Finally the fitted phase shifts are introduced to take into account small uncertainties in the estimate of the phases of the observations, and are typically around ∼0.02.

In the current work, the χ² minimum value was obtained using an analytical approach, rather than relying on commonly used numerical minimization routines (e.g. Levenberg 1944). Specifically, we show that minimizing the χ² function in Eqs. (1) and (2) can be reduced to solving two equations for the unknowns δΦ^j and δΦ_RV, while the remaining fit parameters can be expressed as functions of these phase shifts (see Appendix C). This approach was feasible in our specific case because we were able to derive the first-order analytical derivatives and solve the resulting system of equations explicitly.

The minimization procedure is applied separately to each photometric band and to RV, and finally, the obtained minimum χ² values are combined into one cumulative minimum value according to the following relation:

$$\begin{array}{c}\hfill {\chi }_{\mathit{comb}}^2=1/N_{\mathit{bands}}{\displaystyle \sum _{j=1}^{N_{\mathit{bands}}}}{\chi }_j^2/{(d.o.f.)}_j,\end{array}$$(4)

where N_bands is the number of available bands (including RV), while (d.o.f)_j represents the number of degrees of freedom for the j-th band⁷.

Here we want to focus the attention on the fact that Eq. (4) was minimized at fixed α value, for all the possible values we used to divide the range 1.4–1.9 (see Sect. 3.1). Finally, the optimal α value was determined as described below. For each candidate ${\alpha }_{\mathrm{i}}\mathrm{i}=1\dots 11$, we: selected the corresponding subsample of pulsation models; performed independent fits by minimizing the ${\chi }_{\mathrm{comb}}^2$ function (Eq. (4)); recorded the resulting best-fit ${\chi }_{\mathrm{i}}^2$ value; and chose the final α such that it satisfies ${\chi }_{\mathrm{best}}^2=\min {\chi }_{{\alpha }_1}^2,\dots ,{\chi }_{{\alpha }_{11}}^2$.

5. Monte Carlo methods for best model selection and uncertainty quantification

Monte Carlo techniques encompass a broad class of numerical methods based on the random resampling of a dataset to derive statistical estimates (e.g. confidence intervals and error estimations). In this work we used two different kinds of simulations: (i) bootstrap simulations, obtained by resampling input time series; and (ii) Monte Carlo simulations with observational noise injection. The following subsections detail each method’s implementation and their specific roles in our analysis.

5.1. Bootstrap resampling

Among Monte-Carlo simulations, bootstrap resampling is a widely used approach. Bootstrap methods come in several variants (e.g. parametric, non-parametric, and jackknife), and we refer the reader to the review by Feigelson & Babu (2012) for a comprehensive discussion. In this study, we specifically applied the classical non-parametric bootstrap method. This approach involves generating mock datasets by randomly resampling the original dataset with replacement. Each simulated dataset has the same size as the original set of measurements, but individual data points can appear multiple times, while others may be absent in a given resampled dataset (Feigelson & Babu 2012).

In our work, we performed 250 random simulations. This number was chosen to strike a balance between computational time and the need for a sufficiently large number of trials. Each random experiment consists of two main steps:

• Resampling: The observational time series of both photometry and RV are resampled with replacement, ensuring that each mock dataset retains the same number of data points as the original one.

• Model fitting: The full fitting procedure is repeated for each mock dataset.

This approach serves a dual purpose. First, for a given n-tuple of structural parameters (M, T_eff, log L/L_⊙)⁸, corresponding to a fixed stellar model, the simulations allow us to quantify the uncertainties in the fitted parameters derived from the solutions of Eqs. (C.1), (C.2), (C.3), (C.4), and (C.5) (namely μ_i, ϕ_i, γ, ϕ_rv, and the p-factor). The uncertainty for each fitted parameter is computed as the robust standard deviation (1.4826 · MAD), which rescales the MAD to match the standard deviation of a normal distribution. These values are obtained from the 250 best-fit solutions generated by the Monte Carlo realizations.

Second, for each of the 250 simulations, the optimal pulsation model is identified as the one minimizing the total χ² function defined in Eq. (4). This yields a distribution of optimal models that typically span multiple nodes of the parameter-space grid. Since the location of the minimum χ² varies from simulation to simulation, this dispersion provides a direct measure of the parameter-space degeneracies introduced by observational uncertainties. This behaviour is a fundamental property of χ²: if the phase coverage of the data is uneven, the best-covered phases will dominate the χ² statistic, effectively weighting those regions more heavily during minimization. Our resampling strategy, combined with the inclusion of the peak-to-peak amplitude term (Eq. (3)), helps mitigate this effect.

Focusing on the last itemized point discussed above, we analysed the distribution of minimum χ² values (hereafter χ_min²) obtained from our 250 simulations. A statistical characterization of this distribution enabled us to estimate the expected minimum χ² value. We adopted the median as a robust central estimator, and the 5th–95th percentile range as a measure of its uncertainty. Prior to this analysis, significant outliers – defined as values beyond ±3.5 ⋅ MAD – were excluded to ensure statistical reliability. All models contributing to the clipped χ_min² distribution were then used to infer the optimal structural parameters and their associated uncertainties. This subset defines the best-fitting solutions for the observed source. For each structural parameter (mass, T_eff, and log(L/L_⊙)), the best estimate was taken as the median of its distribution, while the uncertainty was quantified via the semi-range (i.e. half the difference between the maximum and minimum values). This conservative approach was adopted to: (i) provide a robust estimate that reflects the full variability of viable solutions; (ii) avoid assumptions regarding the shape of the underlying distributions; and (iii) prevent the underestimation of the true uncertainties.

5.2. Measurement-error-driven Monte Carlo sampling

In addition to the bootstrap technique introduced in the previous section, this work employs another Monte Carlo-based approach. Specifically, to estimate the uncertainty associated with a given property in our sample of Cepheids – and to perform a statistical analysis that incorporates the measurement errors affecting each input data point (e.g. the p-factor and global parallax offset) – we adopted the following methodology: (i) we conducted 250 Monte Carlo simulations, perturbing each data point within its error bar under the assumption of Gaussian noise statistics, and (ii) for each realization, we computed the mean of the studied quantity. The resulting distribution of the 250 mean values was then analysed, and we defined the 5th-95th percentile interval as the best estimate of the uncertainty for the considered quantity.

6. Results of the fit

The multi-filter Gaia light curves and the velocity curves of the best models found (lines) for all the selected sources, are shown in Fig. D.1, over-imposed to the observed photometric and RV time series. Table D.1 lists the best-fit parameters derived from our analysis. For each star, we report the optimal estimate of the pulsation period (P_best), the corresponding value of the α parameter, and the central values of metallicity and helium content (Z_best, Y_best) characterizing the best-fitting chemical composition bin. We also provide the best-fit structural parameters, namely mass (M), effective temperature (T_eff ), and luminosity (log(L/L_⊙)), as well as the over-luminosity (Δ_ML) with respect to the canonical ML relation. Additionally, the table includes the best-fit p-factor, the γ velocity, the theoretical parallax estimate, and the global χ² value obtained by minimizing Eq. (4)

While not perfect, the theoretical framework provides a reasonable match to both light and RV curves, demonstrating its ability to predict multi-filter variations of Cepheids with different periods and metallicities, despite some systematic deviations. Here we list the main features and issues around the fit, hinting that the problem of stellar pulsation is still not solved:

• Systematic larger modelled amplitudes in the photometric Gaia bands than in the RV curves (CK Cam, V407 Cas, IN Aur, V432 Ori, OZ Cas, V737 Cen, and S Vul): These discrepancies could potentially be reduced by fine-tuning the α parameter, as increasing α might yield models with smaller photometric amplitudes, thereby improving the fit. However, raising α would also decrease the p-factor (see Eq. (2)), which – especially for stars like IN Aur, V432 Ori, OZ Cas, V737 Cen, and S Vul – could push it towards the lower limit of the physically expected range. The observed differences likely stem from limitations in the current treatment of convection-pulsation coupling, particularly the challenges associated with turbulent convection modelling (Marconi et al. 2017; Kovács et al. 2023).

• Minor phase shifts in secondary features (BF Cas, V384 CMa, AP Pup, WX Pup, and OGLE-GD-CEP-0339): These small discrepancies can rise from limitation in the treatment of the convection-pulsation coupling that affects the relative amplitude and position of the primary and secondary maximum.

• Severe mismatches (V339 Cen, HZ Per, SS CMa, AD Pup, RZ Vel, and S Vul): The point detailed above is more evident in these cases, which are closer to the red edge of the instability strip and in some cases affected by limitations in the light curve sampling (see e.g. S Vul).

All the cases listed are characterized by a global χ² value significantly exceeding the expected value of 1, reflecting notable discrepancies between the model and the observations. While these issues are essential to resolve for a complete understanding of the underlying physics, they have a minor impact on the determination of structural parameters. In the following sections, we therefore focus on presenting the results from the analysis of the best-fit model properties.

6.1. Predicted versus measured periods

In Fig. 4 we compare the pulsational periods based on the best-fit models, with the values from Gaia DR3. In the bottom panel we show the percentage differences between modelled periods and observational periods. The median and the robust standard deviation (1.4826 · MAD) of the differences are 1.3% and 1.6%, respectively, while the maximum and minimum differences do not exceed 10%.

Fig. 4. Top panel: Modelled period values (vertical axis) versus the DR3 periods (horizontal axis). The 1:1 line is also plotted to facilitate the comparison. Bottom panel: Percentage difference between the theoretical periods and the DR3 values as a function of the observational period.

6.2. The projection factor

As previously noted, the match between the theoretical and observed amplitudes of the RV curves provides the corresponding stellar p-factor.

Using the inferred results for the p-factor, we investigated its trend with the period by fitting the equation $\mathrm{p}=\alpha +\beta ·log(\mathrm{P})$. The results of this analysis are shown in Fig. 5. In the figure, we have highlighted the three different types of fits that were tested. Despite of the fact that the obtained slopes are significantly different from zero, the very small determination coefficient R² ∼ 0.22 indicates that the correlation is very weak. Although our analysis does not reveal a statistically significant correlation, large uncertainties associated with several p-factor estimates may hinder the detection of a possible underlying trend.

Fig. 5. p-factor as a function of the best-fit pulsational period. Three kinds of linear fits are also reported: the weighted linear fit (blue line), the unweighted linear fit (green line), and the robust linear fit (magenta line). The dashed black line indicates the median p-factor computed by considering the whole sample.

Neglecting the period dependence, we derived the mean p-factor value based on the discussed model fitting results with its associated statistical error: p = 1.22 ± 0.05. The uncertainty on the p-factor was determined through a set of Monte Carlo simulations. First, we varied the p-factor of each individual source within its error, assuming a Gaussian distribution. For eachsimulation, we calculated the mean of the resulting p-factors. Finally, the uncertainty was computed as the 5–95 percentile interval of all the mean values obtained.

The lack of a significant dependence of the p-factor on the pulsation period is consistent with the findings of Nardetto et al. (2009) and Neilson et al. (2012), where the p-factor was calibrated using theoretical models. Similarly, Groenewegen (2007) found no significant period dependence, deriving a constant value of p = 1.27 ± 0.05 by applying the Baade-Wesselink method in its surface brightness form. A comparable p-factor of 1.27 ± 0.06 was also obtained by Mérand et al. (2015) for δ Cep.

Our estimate of the p-factor is in excellent agreement with the values derived for the binary Cepheid OGLE-LMC-CEP-0227 by Marconi et al. (2013a, p = 1.20 ± 0.08) and Pilecki et al. (2013, p = 1.21 ± 0.03 ± 0.04). This binary Cepheid is particularly significant for studying the structural parameters of Cepheids, as its binary nature provides an independent estimate of the primary characteristics of the variable. In Marconi et al. (2013a), the variable was modelled using hydrodynamic pulsation models, while Pilecki et al. (2013) introduced a method specifically designed for analysing double-lined eclipsing binaries containing a radially pulsating star.

Our best estimate is close to the lower end of the assumed p-factor range, but we need to extend the sample in order to equally populate the various period ranges in order to provide more stringent results.

6.3. The α parameter

The determination of p-factor is strongly affected by the choice of the optimal α parameter. Turbulent convection models based on a single-equation formalism, commonly used in radial pulsation studies, often fail to simultaneously reproduce both the RV and photometric variations when a constant p-factor is assumed, as shown by Kovács et al. (2023). This limitation also affects our hydrodynamic code and has similarly been observed in Cepheid models by Marconi (2017).

We notice that, on the basis of our results, no clear trend is found for the variation of α with the effective temperature. However, the sample analysed here might be not large enough to investigate such a dependence, and the large scatter visible in the α versus T_eff plot (see Fig. 6) prevents us from drawing firm conclusions. Moreover, metallicity affects the predicted amplitudes and in turn the selected α values, so that differences in metallicity in some cases contribute to the plot dispersion. It is nevertheless interesting to note that, previous findings in the literature, noted an anti-correlation between the eddy viscosity and the effective temperature for RR Lyrae stars by Kovács et al. (2023, 2024). Similar results were reported for classical and type-II Cepheids by Deka et al. (2024), although without an in-depth analysis of this aspect. While one could argue that those previous results were obtained using different codes – namely the Budapest-Florida Code (Yecko et al. 1998) and MESA-RSP (Paxton et al. 2019; Smolec & Moskalik 2008) – similar trends were also found for RR Lyrae stars and bump Cepheids with the Stellingwerf code (Di Criscienzo et al. 2004; Bono et al. 2002). However, all the investigations adopted a single chemical composition and did not simultaneously take photometric and RV variations into account.

Fig. 6. Best-fit α parameter plotted against the best values of the Teff. Different colours indicate different metallicity bins.

6.4. Absorption

To obtain the reddening-free distance modulus μ₀, we corrected the distance moduli in the individual bands (μ^j parameter in Eq. 1) using the Cardelli extinction law (Cardelli et al. 1989) and the E(B-V) values published in T24. Specifically, we adopted the following equations,

$$\begin{array}{c}\hfill {\mathrm{A}}_{G_{\mathit{BP}}}=0.504·E(B-V)\end{array}$$(5)

$$\begin{array}{c}\hfill {\mathrm{A}}_G=0.582·E(B-V)\end{array}$$(6)

$$\begin{array}{c}\hfill {\mathrm{A}}_{G_{\mathit{RP}}}=0.762·E(B-V),\end{array}$$(7)

to correct the distance moduli μ_GBP, μ_G and μ_GRP, obtained from the fitting procedure, and finally, we computed the weighted average of the three values, which we consider as our best estimate of μ₀.

6.5. Parallax comparison

The reddening free distance moduli, obtained as stated in the previous section, are converted into parallax values (see Col. 14 of Table D.1) and are compared with the Gaia DR3 parallax values (including the correction by Lindegren et al. (2021), hereinafter L21) in Figure 7.

Fig. 7. Top panel: Comparison between the Gaia DR3 parallaxes, corrected according to the L21 recipe, and those obtained from our pulsation model fitting technique, shown together with 1:1 line (grey solid). Bottom panel: Residuals, $\mathrm{\Delta }={\varpi }_{\mathrm{models}}-{\varpi }_{\mathrm{DR}3}$, around the 1:1 line plotted against the DR3+L21 parallaxes.

The analysis presented in the figure indicates that, for parallax values below approximately 0.50 mas, the individual data points show good agreement with the 1:1 line. For parallaxes above this threshold, however, our fitting procedure tends to yield slightly larger values compared to those from DR3+L21. A similar underestimation of parallaxes was reported by Breuval et al. (2020) based on Gaia DR2 data. While the two results are not directly comparable due to differences between DR2 and DR3, the observed behaviour may suggest thepersistence of residual systematics at higher parallax values. According to our results, we find an increasing agreement between the two parallax samples as the parallax value decreases.

To quantify the offset in the Gaia DR3 + L21 parallaxes, we performed a statistical analysis using the theoretical parallaxes as a reference. This analysis involved a series of Monte Carlo simulations similar to those described in Sect. 6.2.

First, the differences (Δ) between the theoretical parallaxes and those from DR3 were computed. The uncertainty on each Δ value was determined as the quadrature sum of the uncertainties on the theoretical parallaxes and those on the DR3+L21 parallaxes. Monte Carlo simulations were then carried out by perturbing the Δ values within their uncertainties, assumed to follow a Gaussian distribution.

The global parallax offset derived from this analysis is 0.014${}_{-0.017}^{+0.022}$ mas, with the uncertainty determined using the 5–95% quantile range. These findings suggest that no significant offset exists in the DR3+L21 parallaxes. A comparable result is obtained by calculating the weighted mean and its associated error from the simulations instead of relying on quantile analysis. In this case, the global parallax offset is 0.005 ± 0.012 mas, further supporting the conclusion that the L21 correction improves the agreement between the Gaia DR3 parallaxes and the reference values.

Other authors have derived the global offset in the DR3+L21 parallaxes. Using eclipsing binaries as reference for distance, Stassun & Torres (2021) obtained a positive parallax offset equal to +0.015 ± 0.018 mas, even if even if consistent with zero, in agreement with our findings.

Wang et al. (2022) used a sample of giant stars with accurate known distance and found that the official L21 model largely improves the parallax determination, with a residualoffset amounting to +0.0026 mas, (+0.0029 mas) for the five-parameter (six-parameter) solutions. Fabricius et al. (2021) compared the EDR3 parallaxes with the distances obtained from other external catalogues (see their Table 1) and found that the L21 correction significantly improves the parallax comparison, in agreement with our results.Significant negative global offset values were found by other authors in the sense that L21 recipe seems to over-correct the parallaxes, returning parallaxes larger than the reference values. For instance, Riess et al. (2021) find a global offset Δϖ= − 0.014 ± 0.006 mas, while two other recent studies, Molinaro et al. (2023) and Bhardwaj et al. (2021), found offsets Δϖ= − 0.022 ± 0.004 mas and Δϖ= − 0.022 ± 0.002 mas, respectively. We refer the reader to Fig. 10 in Molinaro et al. (2023) for a more complete review of the most recent determinations of the global parallax offset.

6.6. Period–Wesenheit–metallicity relation

To minimize the effect of reddening on the distance determination, the reddening-free by definition PW relations and their metal-dependent versions (PWZ) are often adopted (see e.g. T24). The period and metallicity values characterizing the investigated Cepheid sample are not representative of the entire expected range for CCs. As a consequence, we did not fit the PWZ relation by using the mean parameters of the selected variables, but we studied how they distribute in comparison with the PWZ relation recently obtained by our group within the C-MetaLL project (see e.g Ripepi et al. 2021; Trentin et al. 2024).

To this end, we calculated the absolute mean Wesenheit magnitude⁹ in the Gaia bands using the parallax corrections from L21 and the mean magnitudes G, G_BP, and G_RP from the sample of CCs by T24. The obtained values were plotted as a function of period in the top panel of Fig. 8, by adopting a colour scale showing the trend of the [FeH] values. The data representing our sample are also displayed with magenta circles. Additionally, the best-fit PWZ relation from T24 (see their Table 3) is plotted as a dashed black line. The residuals of our sample relative to the best-fit relation are shown in the bottom panel of the same figure. The agreement appears very good. A statistical analysis hints that there is a marginal systematic residual with median value equal to 0.016 mag, compared with the fitted line. Moreover, the robust standard deviation of the distribution of the residuals amounts to 0.10 mag, indicating that our sample is reasonably distributed around the PWZ relation by T24.

Fig. 8. Absolute Wesenheit magnitude plotted against the pulsation period. The fundamental CCs from T24 are also plotted and coloured according to their [FeH] values, together with their best PWZ fitted relation (dashed line). The same quantities for the sample studied in this work are also plotted, together with their errors, with magenta circles. Their residuals around the plotted fit are also shown in the bottom panel.

6.7. Mass–luminosity relation

The inferred individual mass and luminosity values allow us to understand how the best-fit models compare with the ML relations predicted by stellar evolution. In Fig. 9 we show the distribution of the obtained best-fit models (violet circles), divided into the four selected metallicity bins, in the ML plane, compared with the relations predicted by the adopted stellar evolution scenario in the canonical assumption (black solid line), including a mild (+0.2 dex, red dashed line) or full (+0.4 dex, blue dashed line) over-luminosity due to overshooting or mass loss, or rotation. We notice that for each metallicity bin, the bulk of the model results are located within the canonical and the mildly noncanonical relations, with a smaller number of models between the mildly and the fully non-canonical lines. The inferred dispersion supports the idea that the observed Cepheid sample do not obey a unique universal ML relation, possibly reflecting the stochastic nature of mass loss, thus displaying a non-negligible spread in the ML plane. This dispersion is also expected to affect the coefficients of the PL, PLC, and PW relations.

Fig. 9. Mass and luminosity best values obtained by our pulsation model fit technique for all the selected sources (violet points), together with the expected canonical ML relation (solid black line), the mild (+0.2 dex) non-canonical ML relation (dashed red line), and the full (+0.4 dex) non-canonical ML relation (dashed blue line).

7. Conclusions

In this work we used hydrodynamic pulsation models to fit a sample of Galactic CCs for which photometry, RV, and parallax data from Gaia DR3 were available; this was complemented with recent high-resolution spectroscopic metallicities. Several criteria were considered for the sample selection: (i) the widest possible coverage of the expected period range for CCs; (ii) well-sampled light and RV curves in phase; (iii) sufficiently accurate parallaxes (σ_ϖ/ϖ < 50%); and (iv) metallicities obtained through recent high-resolution spectroscopic measurements.

Using the models generated by our hydrodynamic pulsation code, we performed a fitting procedure based on the minimization of the χ² function. By employing a set of bootstrap Monte Carlo simulations, we investigated how the best-fit model varied by randomly resampling the input data. This procedure allowed us to determine the uncertainties on both the fit parameters and the structural parameters estimated for each variable.

The main results of this best-fitting method are:

• The pulsational periods are accurately reproduced by hydrodynamical models, with percentage differences that are not larger than 10%.

• No significant trend was found for the inferred p-factor as a function of the pulsation period, although the relatively large uncertainties associated with several p-factor estimates may hinder the detection of a possible underlying trend. By averaging the obtained values, our best estimate of the p-factor is p = 1.22 ± 0.05. To provide more stringent results, we must extend the investigated sample by uniformly populating the period ranges.

• The parallaxes obtained from the inferred de-reddened distance moduli are in satisfactory agreement with the Gaia DR3 parallaxes corrected according to the L21 recipe, with a residual global offset that is consistent with zero.

• The inferred parallaxes also provide a very good agreement with the PWZ relation by T24.

• The stellar masses and luminosities of the derived best-fit model solutions suggest that, for each selected metallicity bin, the bulk of the investigated Cepheids obey an ML relation that is between the canonical and the mildly over-luminous (+0.2 dex) one, with a smaller number of pulsators displaying a brighter relation.

The results of the presented analysis confirm the global predictive capability of the adopted theoretical scenario. We plan to investigate residual uncertainties more thoroughly, for example the ones connected to the reddening or p-factor determination, by extending both the observational sample and parameter space of the adopted model set and making use of the data that will be provided by upcoming facilities devoted to variability studies, such as the Rubin LSST.

Data availability

Appendix D is available on Zenodo.

Acknowledgments

The Authors thank the Referee, for his detailed analysis of this work, that impressively improved its contents and readability.

References

Appendix A: Initial Fourier model

After folding the input data with the observed pulsational period, the fitting routine constructs a Fourier model for each input time series of photometric and radial velocity data. The routine uses some useful information provided by this initial Fourier model. In particular, the phase difference between the maximum of light (minimum of velocity) of the observed time series and that of the theoretical curves, gives a first rough guess of the phase shift to be applied to the model in order to fit observations. Furthermore, the Fourier peak-to-peak amplitude is used by the fitting routine to compute the χ_Amp² term of Eqs. 1 and 2.

Appendix B: The PTLM and ML equations

The PTLM relation adopted in this work is from De Somma et al. (2020):

$$\begin{array}{c}\hfill logP=a_1+a_{2}log{T}_{\mathit{eff}}+a_{3}log(M/M_{\odot })+a_4(log(L/L_{\odot })+{\mathrm{\Delta }}_{\mathit{ML}}),\end{array}$$(B.1)

where the coefficients a_i depend on the selected α value (see Table B1), while the Δ_ML takes into account the over-luminosity level with respect to the canonical ML relation (Δ_ML = 0), having the following form (Bono et al. 2000b):

$$\begin{array}{c}\hfill logL=b_1+b_{2}logM+b_{3}logZ+b_{4}logY,\end{array}$$(B.2)

where the coefficients are given in Table B1.

A generic model is defined by its elemental composition (Z, Y), the pulsational period (P), the effective temperature (T_eff), the mass (M), and the luminosity (log L/L_⊙). Fixing the elemental composition, for every mass (M) and effective temperature values (T_eff), the luminosity of the model is given by summing the canonical value, obtained from Eq. B.2, to the non-canonical correction Δ_ML given by:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathit{ML}}=\frac{1}{a_4}(logP-a_1-a_{2}log{T}_{\mathit{eff}}-(a_3+a_4{b}_2)log(M/M_{\odot })\\ \hfill -a_4{b}_{3}logY-a_4{b}_{4}logZ-a_4{b}_1)\end{array}\end{array}$$(B.3)

such that the computed theoretical period matches the value obtained from the observations.

We note here that, we did not consider the PTLM coefficient errors in our computation because they affect the inferred parameters by amounts that are significantly smaller than the adopted steps in the corresponding model grid.

Appendix C: Minimization of the χ²

Requiring that the first derivatives of the eq. 1 and eq. 2 with respect to the parameters of the fit, are equal to zero, it is possible to find the following equations:

$$\begin{array}{c}\hfill -{\displaystyle \sum _{i=1}^{N_{\mathit{pts}}^j}}\left(\frac{m_i^j({\varphi }_i^j){\dot{M}}^j({\varphi }_i^j+\delta {\mathrm{\Phi }}^j)}{{\left({\sigma }_i^j\right)}^2}\right)+{\displaystyle \sum _{i=1}^{N_{\mathit{pts}}^j}}\left(\frac{M^j({\varphi }_i^j+\delta {\mathrm{\Phi }}^j){\dot{M}}^j({\varphi }_i^j+\delta {\mathrm{\Phi }}^j)}{{\left({\sigma }_i^j\right)}^2}\right)+\frac{1}{B}({\displaystyle \sum _{i=1}^{N_{\mathit{pts}}^j}}\left(\frac{M^j({\varphi }_i^j+\delta {\mathrm{\Phi }}^j)}{{\left({\sigma }_i^j\right)}^2}\right)-A){\displaystyle \sum _{i=1}^{N_{\mathit{pts}}^j}}\left(\frac{{\dot{M}}^j({\varphi }_i^j+\delta {\mathrm{\Phi }}^j)}{{\left({\sigma }_i^j\right)}^2}\right)=0\end{array}$$(C.1)

$$\begin{array}{cc}\hfill C{\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}{\left(\frac{V_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{\sigma }_l^{\mathit{RV}}}\right)}^2{\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\left(\frac{v_l{\dot{V}}_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{\left({\sigma }_l^{\mathit{RV}}\right)}^2}\right)-{\left({\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\left(\frac{V_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{\left({\sigma }_l^{\mathit{RV}}\right)}^2}\right)\right)}^2{\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\left(\frac{v_l{\dot{V}}_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{\left({\sigma }_l^{\mathit{RV}}\right)}^2}\right)& \hfill \\ \hfill -C{\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\left(\frac{v_l{V}_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{\left({\sigma }_l^{\mathit{RV}}\right)}^2}\right){\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\left(\frac{{\dot{V}}_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})V_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{\left({\sigma }_l^{\mathit{RV}}\right)}^2}\right)+& \hfill \\ \hfill +A^{\prime }{\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\left(\frac{V_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{\left({\sigma }_l^{\mathit{RV}}\right)}^2}\right){\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\left(\frac{{\dot{V}}_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})V_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{\left({\sigma }_l^{\mathit{RV}}\right)}^2}\right)+& \hfill \\ \hfill -A^{\prime }{\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}{\left(\frac{V_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{\sigma }_l^{\mathit{RV}}}\right)}^2{\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\left(\frac{{\dot{V}}_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{\left({\sigma }_l\right)}^2}\right)+& \hfill \\ \hfill +{\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\left(\frac{V_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{\left({\sigma }_l^{\mathit{RV}}\right)}^2}\right){\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\left(\frac{v_l{V}_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{\left({\sigma }_l^{\mathit{RV}}\right)}^2}\right){\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\left(\frac{{\dot{V}}_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{\left({\sigma }_l\right)}^2}\right)=0\end{array}$$(C.2)

Here $A={\displaystyle \sum _{i=1}^{N_{\mathit{pts}}^j}}\frac{m_i^j}{{\left({\sigma }_i^j\right)}^2}$, $B={\displaystyle \sum _{i=1}^{N_{\mathit{pts}}^j}}\frac{1}{{\left({\sigma }_i^j\right)}^2}$ and $C={\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\frac{1}{{\left({\sigma }_l^{\mathit{RV}}\right)}^2}$, $A\prime ={\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\frac{v_l}{{\left({\sigma }_l^{\mathit{RV}}\right)}^2}$, while the unknown parameters are the two phase shifts, δΦ^j and δΦ_RV.

After solving the equations for the two unknown phase shifts, the apparent distance modulus in the j-th band is given by

$$\begin{array}{c}\hfill {\mu }^j=\frac{1}{B}·(\frac{{\displaystyle {\sum }_{i=1}^{N_{\mathit{pts}}^j}}{M}^j({\mathrm{\Phi }}_i^j+\delta {\mathrm{\Phi }}^j)}{{\sigma }_i^2}-A),\end{array}$$(C.3)

while the barycentric

γ velocity and the p-factor can be obtained from the following equations:

$$\begin{array}{cc}& \gamma =\frac{A^{\prime }·(E-F)+H·B^{\prime }-D·B^{\prime }}{B^{\prime }·G-B^{\prime 2}-C·F+C·E},\hfill \end{array}$$(C.4)

$$\begin{array}{cc}& 1/p=\frac{A^{\prime }-\gamma ·C}{B},\hfill \end{array}$$(C.5)

where $B\prime ={\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\frac{V_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{({\sigma }_l^{\mathit{RV}})}^2}$, $D={\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\left(\frac{v_l·V_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{({\sigma }_l^{\mathit{RV}})}\right)$, $E={\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\frac{V_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{\sigma }_l^{\mathit{RV}}}$, $F={\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\frac{V_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})·\dot{V}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{\sigma }_l^2}$, $G={\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\frac{{\dot{V}}_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{\sigma }_l^2}$, and $H={\displaystyle \sum _{l=1}^{M_{\mathit{pts}}}}\frac{v_l·{\dot{V}}_{\mathit{puls}}({\varphi }_l+\delta {\mathrm{\Phi }}_{\mathit{RV}})}{{\sigma }_l^2}$. Despite their long and complicated expressions, Eqs. C.3, C.4, and C.5 consist of summations and are very fast to compute once the δΦs are obtained from the solution of Eq. C.2.

All Tables

All Figures

	Fig. 1. Distribution of the selected classical Cepheids period values obtained from the DR3.
In the text

	Fig. 2. Distribution of the selected classical Cepheids [Fe/H] values obtained through high-resolution spectroscopy (T24).
In the text

	Fig. 3. Distribution of the Z values of all the selected variables computed by converting the spectroscopic [Fe/H] values from Trentin et al. (2024).
In the text

	Fig. 4. Top panel: Modelled period values (vertical axis) versus the DR3 periods (horizontal axis). The 1:1 line is also plotted to facilitate the comparison. Bottom panel: Percentage difference between the theoretical periods and the DR3 values as a function of the observational period.
In the text

	Fig. 5. p-factor as a function of the best-fit pulsational period. Three kinds of linear fits are also reported: the weighted linear fit (blue line), the unweighted linear fit (green line), and the robust linear fit (magenta line). The dashed black line indicates the median p-factor computed by considering the whole sample.
In the text

	Fig. 6. Best-fit α parameter plotted against the best values of the Teff. Different colours indicate different metallicity bins.
In the text

	Fig. 7. Top panel: Comparison between the Gaia DR3 parallaxes, corrected according to the L21 recipe, and those obtained from our pulsation model fitting technique, shown together with 1:1 line (grey solid). Bottom panel: Residuals, $\mathrm{\Delta }={\varpi }_{\mathrm{models}}-{\varpi }_{\mathrm{DR}3}$, around the 1:1 line plotted against the DR3+L21 parallaxes.
In the text

Fig. 8. Absolute Wesenheit magnitude plotted against the pulsation period. The fundamental CCs from T24 are also plotted and coloured according to their [FeH] values, together with their best PWZ fitted relation (dashed line). The same quantities for the sample studied in this work are also plotted, together with their errors, with magenta circles. Their residuals around the plotted fit are also shown in the bottom panel.

In the text

	Fig. 9. Mass and luminosity best values obtained by our pulsation model fit technique for all the selected sources (violet points), together with the expected canonical ML relation (solid black line), the mild (+0.2 dex) non-canonical ML relation (dashed red line), and the full (+0.4 dex) non-canonical ML relation (dashed blue line).
In the text