Introduction

Cepheids are the best-established standard candle. They link the distance scale in the Local Group with type Ia supernova host galaxies. A thorough understanding of the pulsation of these stars is required to obtain the best accuracy on the Hubble constant H 0 (Breuval et al. 2020;Riess et al. 2021).

Obtaining accurate distances to Cepheid stars is still a nontrivial issue. Cepheid distances may be derived through main-sequence fitting for Cepheids in clusters or through the measurement of their parallax. Recently, very precise geometric parallaxes for about 9500 Cepheids were measured by the Gaia satellite (Gaia Collaboration 2020), which is the first com-petitive alternative to Hubble Space Telescope (HST) parallaxes (Benedict et al. 2007;Riess et al. 2018).

In addition, distances to classical Cepheids (CCs) can be obtained from the parallax-of-pulsation method (PoP). In this approach, the variation in the angular diameter of a Cepheid is compared with the variation of its linear diameter, derived from the integration of its pulsation velocity. The true pulsational velocity of a star is derived by multiplying the disk-integrated radial velocities (measured by spectroscopy) by a projection factor (hereafter p-factor). In the absence of interferometric measurements, angular diameters can be derived from surfacebrightness-color relations (SBCR): this particular implementation of the PoP technique is known as the Baade-Wesselink (BW) method (Baade 1926;Wesselink 1946). The PoP method is the most geometrical way, except for measuring the direct parallax, to estimate the distance of Cepheids. This method is therefore valuable in calibrating the period-luminosity (P-L) relation, also called the Leavitt law (Leavitt & Pickering 1912). However, the assumptions behind the PoP method may introduce strong sources of error on the derived distances. Especially the current uncertainty on the p-factor value is still the main reason for recent determinations of the Hubble constant based on Cepheid distances to avoid relying on the PoP technique (Riess et al. 2009). Mérand et al. (2015) developed the code called spectrophoto-interferometry of pulsating stars (SPIPS). This is a variant implementation of the PoP method that uses atmospheric models and combines all types of available data in order to bypass the limitations of the traditional BW method that affect the accuracy and precision of the derived parameters of a pulsating star. Unfortunately, previous studies using this method (Mérand et al. 2015;Breitfelder et al. 2016;Kervella et al. 2017;Gallenne et al. 2017;Trahin 2019) or alternatives (Ngeow et al. 2012;Storm et al. 2011;Pilecki et al. 2018) did not converge to a consistent dependence of the p-factor because the few available HST parallaxes were not very precise, because of the Gaia DR2 zeropoint uncertainty, or because the datasets were incomplete.

In this paper, we present the application of the SPIPS method to a sample of CCs for which we used the best and most complete data, in combination with the new Gaia EDR3 parallaxes, and we derive various precise and consistent parameters and investigate their dependences. This paper is similar to the study by Gallenne et al. (2017), who performed a SPIPS analysis of Large Magellanic Cloud (LMC) and Small Magellanic Cloud (SMC) Cepheids for which they disposed of light curves in order to derive the period-p-factor relation. The difference is that our work is based on Milky Way Cepheids and uses a larger set of data (effective temperatures, more complete photometry, and radial velocities and diameters).

In Sect. 2 we introduce our sample of 63 Galactic Cepheids with their data and present the SPIPS method. In Sect. 3 we adopt Gaia EDR3 parallaxes as an input into the SPIPS algorithm and apply this method to our sample of Cepheids. Our calculations converge to a robust estimate of their parameters such as radius, reddening, mean multiband magnitudes, effective temperature, and p-factor. Finally, in Sect. 4 we test the accuracy of the parameters derived from the SPIPS modeling by calibrating the P-L and period-radius (P-R) relations, and we investigate the dependences of the projection factor.

Cepheid data and fitting method

Cepheid sample and data

We built a database including most of the observations collected in the past 50 years for more than 300 Cepheids (including our own observations) in order to identify the stars with the best dataset. The realization of the resulting database was made possible using the McMasterfoot_0 , Vizier (Ochsenbein et al. 2000), Simbad (Wenger et al. 2000), AAVSOfoot_1 , and ADSfoot_2 databases. For the application of the PoP technique, we only used a subset of this database for which the data were ideal. We assumed that a suitable dataset corresponds to a full phase-coverage, which is associated with a good accuracy and a minimum dispersion of the data. Moreover, we required that all data had the corresponding epoch of observation. The Modified Julian Date (MJD) of an observation allowed us to determine the period and the period changes of the star with the best precision, which is not possible with the indication of the phase alone. This preliminary selection led to a sample of 63 Cepheids that covers a broad range of periods from 3 to 68 days. This sample constitutes one of the most complete, precise, and homogeneous samples of Galactic Cepheids that are available for the application of the PoP method. The references of all the data we used are provided in Table A.1.

Photometry. Each Cepheid of our sample has at least photometric data in the optical B and V bands (which contain the information about the temperature and reddening), and in the near-infrared (NIR) J, H, K bands (which are less sensitive to interstellar reddening and are more sensitive to the variation in radius). We also used photometric data from spatial observatories and surveys such as Hipparcos (H p band) and Tycho (B and V bands), Spitzer (I 1 and I 2 bands), 2MASS (J, H, and K S bands), and Gaia (G, BP, and RP bands). As recommended in Breitfelder et al. (2016), we did not use R-and I-band photometry because the effective bandpasses are poorly defined. In the SPIPS algorithm, all photometric observations are modeled with the dedicated filters available in the Spanish Virtual Observatory databasefoot_3 (SVO, Rodrigo & Solano 2020). Most of the data used in this study are originally in the California Institute of Technology (CIT) system. However, only the South African Astronomical Observatory (SAAO) filters are not available in the SVO database. We converted the infrared photometry from the SAAO system into the CIT system in order to include these data and to obtain a better phase coverage of the NIR photometry. We used the following equations from Carter (1990):

In Sect. 4.1 we perform a second transformation of mean magnitudes from the CIT to the 2MASS system in order to compare our PL relations in infrared bands with other calibrations from the literature. In the SPIPS adjustments, we decided to keep the data (except for the SAAO data) in their original system as far as possible in order to avoid introducing potential systematics in the derived parameters. For safety, we introduced a conservative systematic uncertainty of 0.01 mag in order to take the different instrumental calibrations and photometric zeropoints into account. This value is consistent with the average offset that is generally observed when data from different instruments and magnitude systems are combined (see, e.g., Barnes et al. 1997;Breitfelder et al. 2016).

Radial velocities. The p-factor depends on the method that is used to extract the radial velocity (such as cross-correlation or broadening functions) because the velocity curves that are obtained with different techniques can have a difference of up to 5% in amplitude (Nardetto et al. 2009). This must be taken into account for studies that use the p-factor, in particular regarding its dependence on other parameters such as the period. In this work, we only used radial velocities determined from crosscorrelation techniques. As the p-factor directly depends on the integrated radial velocity curve, we took care to use only precise observations with full phase coverage and with a well-defined amplitude. As observed in Kervella et al. (2019a), at least 80% of the Cepheids belong to a multiple system. For most stars of our sample, binary Cepheids are not excluded, but the effect on the radial velocities and photometry is considered to be negligible. For some Cepheids, radial velocities are clearly affected by a spectroscopic companion and were corrected for the Keplerian motion using the orbital parameters from the Konkoly databasefoot_4 in order to retain only the pulsation component. They are indicated by a star in Table A.1.

A conservative uncertainty of 0.5 km s -1 was quadratically included as a systematic error in order to take all the systematic effects due to the combination of different datasets into account.

Effective temperatures. For some stars, we disposed of effective temperature measurements, which are mostly provided by the series of papers by Luck & Andrievsky (2004), Andrievsky et al. (2005), Kovtyukh et al. (2005), Luck et al. (2008), Luck (2018), andProxauf et al. (2018). In these papers, the authors estimated the depth ratio of about 50 spectral lines (described in Kovtyukh & Andrievsky 1999) in order to derive the effective temperature of the star. These observations allowed us to constrain the SPIPS models better and to evaluate the consistency of the atmospheric models. We included an error of 50 K as a systematic error for the effective temperatures (Breitfelder et al. 2016).

Angular diameters. In the past ten years, improvements in interferometry enabled the direct determination of the angular diameter for some Cepheids. Several stars of this sample were regularly observed with the CHARA and VLTI interferometers in order to obtain direct measurements of their angular diameter variations. These observations associated with the SPIPS method already allowed us to obtain a better precision on the projection factor (Breitfelder et al. 2016). The new raw data that we obtained with the PIONIER instrument of the VLTI were reduced using the pndrs data reduction software (Le Bouquin et al. 2011). We then adjusted the calibrated squared visibilities with a uniform disk (UD) model to obtain the UD angular diameters. A conservative uncertainty corresponding to 2% of the angular diameter values was quadratically added as a systematic (Kervella et al. 2004).

Distances. As input in the SPIPS code, we adopted the parallaxes from Gaia EDR3 (Gaia Collaboration 2020) and inverted them to obtain the Cepheid distances. We note that using the Bailer-Jones et al. (2021) approach to derive distances of 9000 Cepheids, the geometric distance (based on the parallax and on the direction on the sky) and the photo-geometric distance (which also includes the color and apparent magnitude of the star) differ by 9 pc from the inverted Gaia EDR3 parallaxes on average, with a largest difference of 230 pc. This comparison confirms that the inversion of Gaia EDR3 parallaxes in order to obtain Cepheid distances does not add biases to the computed parameters. We corrected each parallax for their individual zeropoint offset by using the dedicated Python codefoot_5 described by Lindegren et al. (2021a). Alternative procedures to determine the zeropoint offsets were realized (Groenewegen 2021, e.g., suggested that the Gaia EDR3 parallaxes may be underestimated by about 5%), but we limit this study to the Lindegren et al. (2021a) corrections. The new Gaia EDR3 catalog also provides the renormalized unit weight error (RUWE) indicator, which represents the quality of a star's parallax compared with other stars of the same type. Lindegren et al. (2021a) recommended to avoid the use of parallaxes with a RUWE indicator higher than 1.4. We find 16 stars in this case in our sample of 63 Cepheids: we performed the SPIPS modeling successfully for these stars, but we did not use them to calibrate the P-p, P-L, and P-R relations, which depend on the distance. However, we made an exception for δ Cep, for which we had one of the best available datasets, with a full coverage of the interferometric angular diameters and spectroscopic effective temperature curves. For this star, Kervella et al. (2019b) found a close companion that has a very precise Gaia EDR3 parallax with a RUWE of 1.415, which is only slightly higher than the threshold for the other stars of our sample. Only one other star (RS Pup) has a similar dataset, which permits constraining the different parameters better.

The range of magnitudes G = [10.8 -11.2] corresponds to a transition of window classes (see Fig. 1 in Lindegren et al. 2021b) that might affect the accuracy of the zeropoint offset, but none of our stars falls in this range. Finally, we followed the conservative recommendation by Riess et al. (2021) and increased each parallax error by 10% to account for potential additional excess uncertainty.

SPIPS fitting method

We used the SPIPS modeling toolfoot_6 from Mérand et al. (2015) to reproduce our observational dataset. This algorithm is inspired by the classical BW technique. We here present the general idea of the SPIPS method and refer the reader to Mérand et al. (2015) for more details.

The motivation behind the SPIPS method is to bypass the limitations of the traditional BW implementation, which affect the accuracy and precision of the derived parameters. A main limitation of the BW method results from the determination of angular diameters through surface brightness-color relations using only two photometric bands (generally V and K). In this case, the effective temperature and the angular diameter of the star are adjusted from only two photometric measurements. Finally, a poor phase coverage or a low-order interpolation of the different quantities can prevent the precise determination of the parameters.

The approach of the SPIPS method is first to propose a combination of all the data available in the literature for a star. This includes spectroscopic radial velocities as well as photometric measurements in any filter and optical interferometric measurements. In the current code, we use radial velocities derived from cross-correlation. A future implementation is in progress to directly reproduce the spectral lines from highresolution spectroscopy to derive RVs, effective temperatures, and other parameters. The data are then adjusted simultaneously altogether, using a standard multiparameter χ 2 minimization, in order to obtain more realistic estimates of the statistical uncertainties, as opposed to a method that would fit consecutive sets of parameters. The SPIPS code also determines the period as well as the period changes of the pulsation by phasing the data. The BW method generally makes the assumption that empirical surface-brightness relations are linear in color (e.g., V -K), which propagates a color bias on the distance. In order to bypass these uncertainties, SPIPS computes the specific surface brightness using a grid of ATLAS9 atmospheric models 8 (ignored) Fig. 1. Result of the SPIPS modeling for the Cepheid CD Cyg. The fitted observed data include radial velocities (top left), effective temperatures (middle left), and multiband photometry (right). Some main parameters derived from this modeling are listed above the plot. The distance is fixed to the Gaia EDR3 parallax with the zeropoint correction by Lindegren et al. (2021b). (Castelli & Kurucz 2004) to derive synthetic photometry from the effective temperature. The photometric magnitudes are then computed on this grid, using bandpasses and zeropoints from the SVO database. If interferometric observations of the angular diameter of a star are available, the effects of the limb darkening have to be taken into account: in the SPIPS algorithm, the uniform disk angular diameters estimated from the observed visibilities are converted into limb-darkening values using SATLAS 9 spherical atmosphere models (Neilson & Lester 2013).

The interstellar reddening is parameterized in SPIPS using the B -V color excess E(B -V) and the reddening law from Fitzpatrick (1999), adopting R V = 3.1. As explained in Mérand et al. (2015), the reddening corrections in SPIPS are computed on the basis of photometric observations of the star, whereas in classical implementations of the BW method, they are usually computed for a Vega-like star, which is much hotter (10 000 K) than Cepheids (∼5000 K). Moreover, a circumstellar envelope (CSE) is a frequent phenomenon around massive pulsating stars such as Cepheids (Hocdé et al. 2020;Gallenne et al. 2021). It introduces a bias on the interferometric angular diameters and the NIR photometric measurements. The latter are characterized by a magnitude excess and are taken into account in SPIPS by adjusting a power law for the infrared excess, assuming that there is no excess in optical wavelengths (λ < 1.2 µm). This law is defined as IR ex = α(λ -1.2) 0.4 , where α is the slope of the relation and λ is the wavelength of the photometric band in µm. Regarding interferometric measurements, the bias due to a CSE depends on the baselines and on the angular diameter. It is tabulated in SPIPS as a function of the infrared excess. The 9 SATLAS: http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/ A+A/554/A98 visibilities of a limb-darkened disk surrounded by a CSE are synthesized, and a uniform disk model is then adjusted to estimate the bias.

For each CC, the phases of the data points were calculated using their corresponding modified Julian date epoch of maximum light (MJD 0 ). A strategic approach to fitting the data was to start from a model whose general properties were close to the observed data so that the model fitting would converge faster. These starting model curves are third-order Fourier series whose amplitudes and phases agree with the data. They were built from a set of parameters found in the literature (e.g., the mean effective temperature, the period, and the MJD 0 ) and by computing mean values of the radial velocity and of the angular diameter from the available data and distance of the star. Depending on the properties of the different curves (e.g., bumps or steep variations), we then adapted the number of Fourier modes and thus of free parameters to obtain a satisfactory representation of the observed variations of the star. The Fourier series decomposition is a robust method for studying the light curves of variable stars. Fourier coefficients and parameters are commonly used today to model a Cepheid light curve (Morgan et al. 2007;Deb & Singh 2009;Bhardwaj et al. 2015). The third-order coefficients in the K band are listed for each star in Table E.1. These results can be a used in future comparative studies aiming at constraining theoretical stellar pulsation models and determining pulsation modes of Cepheids.

Results of the SPIPS modeling

The SPIPS model fitting was performed for each of the 63 Cepheids of our sample. The final SPIPS adjustment for the Cepheid CD Cyg is presented in Fig. 1. The available data for Notes. Nonfitted parameters are indicated by a star ( ). Systematic uncertainties are taken into account in the listed errors (see also Sect. 2).

Main parameters derived by the SPIPS algorithm

The SPIPS algorithm returns various parameters from the modeling of each of the 63 Cepheids, such as E(B -V) values, dereddened apparent and absolute magnitudes, mean radius of the star, infrared excess, and the projection factor. These parameters are provided in Table 1. The values agree well with those of Gallenne et al. (2021), who derived parameters of 45 Galactic Cepheids using the SPIPS algorithm, with a fixed p-factor and a different dataset than the one adopted here, showing the robustness of the method. In particular, they reported that the IR excess of nearly 30% of the Cepheids is likely produced by a CSE. We refer to this study for a detailed analysis of this effect. For our sample of stars, the uncertainties on IR excess are rather large and do not allow us to conclude about the presence of a circumstellar envelope at this stage. Adopting a precision threshold of 30% of the IR excess value, which corresponds to the most precise values of our sample, leads to approximately the same fraction of Cepheids with detected CSEs as was observed by Gallenne et al. (2021). In addition, given the large size of our sample, we are able to exclude a correlation between the IR excess and the period that was suggested by dusty-wind models. Hocdé et al. ( 2020) have proposed free-free emission as an explanation for the formation of circumstellar envelopes.

The mean apparent magnitudes listed in Table C.1 correspond to flux-averaged mean magnitudes. The B and V magnitudes are in the Cousins and Johnson systems, respectively. The NIR J, H, and K S mean magnitudes, originally in the CIT system, were converted into the 2MASS system using the following transformation relations from Monson & Pierce (2011), with negligible transformation errors: The corresponding scatter in J, H, and K s bands is σ = 0.018, 0.014, and 0.014, respectively.

As mentioned in Sect. 2.2, the reddening E(B -V) can be computed directly by the SPIPS algorithm instead of using values from the literature. Most values from literature are imprecise or derived from inhomogeneous methods, which can have a substantial effect on the consistence of the derived values and on the accuracy of the P-L relation calibration. We represent in Fig. 2 a comparison of the reddening values computed by the SPIPS code with those from the Stilism 3D map (Lallement et al. 2018) or derived by Kovtyukh et al. (2008) and Fernie et al. (1995). The latter are used extensively in calibrations of the periodluminosity relation (Groenewegen 2018;Breuval et al. 2021).

The dispersion between all these values is significant, but we note that the SPIPS reddening values generally agree with the others and are larger by 0.05 mag on average. The largest differences may be explained by the fact that as mentioned in Sect. 2.2, the reddening derived by SPIPS is based on the photometry of each star, whereas in the literature, it is obtained for a Vegalike star, which is significantly hotter than Cepheids. Hocdé et al. (2020) also pointed out that ionized-gas envelopes are likely to obscure Cepheids in the optical bands from 0.05 to 0.15 mag, approximately. This difference of color excess might be due to a compensation in the SPIPS algorithm for the absence of a circumstellar envelope in the optical by an increase in E(B -V). This would especially be the case for the models in which angular diameter, distance, and effective temperatures are constrained by the observations.

Position in the Hertzprung-Russell diagram

The effective temperatures and luminosities derived from our SPIPS analysis allowed us to precisely determine the variation in position of the studied Cepheids in the Hertzsprung-Russell diagram during their pulsation phase (Fig. 3). Figure 3 shows the very dynamic nature of these objects, which move significantly outside of the instability strip during a pulsation. However, the mean values of the effective temperature and the luminosity show that these objects are mainly confined between the blue and red edges defined by Anderson et al. (2016a). Moreover, the Cepheids of our sample are closer to the blue edge on average. Although the stars at the center of the strip seem to have higher amplitudes, no strong correlation between the amplitude and the proximity of one of the edges is visible, as was reported by Fernie (1990).

Our data cover almost 50 years of observations, which also allowed us to derive new period change rates Ṗ (listed in Table D.1). Negative period changes arise during the second crossing of the instability strip, and positive period changes correspond to a first or third crossing. These values are consistent with the predictions by Fadeyev (2014). We can note in particular that most of our Cepheids are in their second or third crossing of the instability strip. Figure 4 shows the linear dependence of the logarithm of the period change on the period. The scatter of this relation is mainly explained in Anderson et al. (2016a) by the rotation dependence of Ṗ. Miller et al. (2020) showed that rotation is insufficient to explain this distribution of period change rates, and that other mechanisms such as mass loss are required.

Discussion

Period-luminosity relation from SPIPS absolute magnitudes

The Cepheid period-luminosity relation is of primary importance for measuring astronomical distances. In most recent studies, this relation is the foundation of the extragalactic distance scale on which the determination of the local Hubble constant H 0 (Breuval et al. 2020;Javanmardi et al. 2021;Riess et al. 2021) is based. Using the mean apparent magnitudes and color excesses derived by the SPIPS models, we computed the absolute magnitudes and the astrometry-based luminosities (ABL, Arenou et al. 1999) from the Gaia EDR3 parallaxes with RUWE < 1.4. Apparent magnitudes were corrected for the extinction using the reddening law A λ = R λ E(B -V) with R V = 3.10, R J = 0.815 and R K S = 0.351 (Fitzpatrick 1999). A reddening-free Wesenheit magnitude W JK was also derived, defined by W JK = K S -0.756(J -K S ).

We then performed a weighted fit of the ABL function and ensure the robustness of the fit by using a Monte Carlo approach with 10 000 iterations. The absolute magnitudes were parameterized around the pivot period log P 0 = 0.9 such as M λ = b λ + a λ (log P -0.9) in order to reduce the correlation between a λ and b λ and to minimize their respective uncertainties. We accounted for the width of the instability strip by adding in quadrature an additional term of 0.07 mag and 0.22 mag in K S and V, respectively, in the magnitude errors listed in Table C.1.

The derived PL relation in the K S , V, and W JK bands are represented in Fig. 5. The best-fit solution in the K S band corresponds to M K S = -5.529 ±0.015 -3.141 ±0.050 (log P -0.9).

(1)

In V and W JK bands, the best-fits are M V = -4.025 ±0.036 -2.775 ±0.103 (log P -0.9) (2)

W JK = -5.803 ±0.009 -3.199 ±0.034 (log P -0.9), respectively. (3)

The dispersion in the V, K S , and W JK bands is 0.22, 0.18, and 0.18 mag, respectively. We note that our K S -band calibration agrees excellently with the result by Breuval et al. ( 2020

Period-radius relation from SPIPS radii

In recent years, the calibration of the period-luminosity relation has been given particular importance. However, the periodradius relation of Cepheids also plays an important role in determining the masses and various parameters of these stars. As stated by Gieren et al. (1998), this relation may also be used to derive pulsational parallaxes of Cepheids in galaxies in which radial velocity curves cannot be observed. For this purpose, we computed the radius of each star from the angular diameter curves modeled by the SPIPS algorithm using Gaia EDR3 parallaxes for Cepheids with a RUWE < 1.4. We derive the following period-radius relation of Galactic Cepheids, represented in Fig. 6: log R = 1.763 ±0.003 + 0.653 ±0.012 (log P -0.9).

This relation has a relatively low dispersion (σ = 0.04) and agrees well with the red and blue edges of the instability strip defined by Anderson et al. (2016c). We also note that it is compatible with the relation defined by Molinaro et al. (2011) at short periods (log P < 1), and with the relation by Gallenne et al. (2017) established for LMC Cepheids.

Period-p-factor relation

Many studies recently made use of the parallax-of-pulsation method with the intention of calibrating the period-luminosity relations (Fouqué et al. 2007;Storm et al. 2011;Groenewegen 2013;Breitfelder et al. 2016;Kervella et al. 2017;Gieren et al. 2018;Trahin 2019). As discussed in the introduction, the projection factor is still the main limitation of this method to derive accurate distances that are competitive with geometrical parallaxes. Although the physics behind this parameter is better understood nowadays through the various works by Nardetto (2005), Nardetto et al. (2006Nardetto et al. ( , 2007Nardetto et al. ( , 2009Nardetto et al. ( , 2011Nardetto et al. ( , 2017)), numerous effects are still blurry. The limb-darkening is more important for the most massive stars (i.e., stars with a longer period), therefore most studies tend to conclude with a linear dependence of the p-factor on the period with a negative slope. Dynamical effects in the pulsating atmosphere might play a role as well. However, Fig. 7 clearly shows the disparity of the P-p relations found in the literature.

The unprecedented precision of the recent Gaia EDR3 parallaxes is a major tool in order to break the degeneracy of the distance over p-factor ratio in the PoP method and to constrain this parameter. Using the SPIPS implementation of the PoP technique described in Sect. 2.2, we computed the value of the projection factor for each star of our sample (with RUWE < 1.4). These values are listed in the main Table 1 and are represented as a function of the logarithm of the period in Fig. 8. We point out that the radial velocities of six Cepheids with a RUWE < 1.4 are highly affected by a spectroscopic binary that can bias the results, therefore we excluded them from the fit. These stars are AW Per, VZ Cyg, V0636 Sco, X Sgr, MW Cyg, and Z Lac. Two stars (δ Cep with the Gaia EDR3 parallax of its companion and RS Pup) have a complete dataset with a full phase coverage of interferometric angular diameters, effective temperatures, radial velocities, and multiband photometry. We note the high (Gieren et al. 2005;Nardetto et al. 2009;Feast et al. 2008;Storm et al. 2011;Ngeow et al. 2012;Neilson et al. 2012;Groenewegen 2013;Gallenne et al. 2017;Kervella et al. 2017;Trahin 2019). The gray region represents the dispersion of the values around p = 1.26 derived in this study. dispersion of the p-factor values and also some unexpected values with p > 1.5 (area delimited in gray in Fig. 8), which would physically correspond to a limb brightening of the stellar disk (instead of a limb darkening) or a reverse atmospheric velocity gradient (increase in velocity amplitudes toward the upper part of the atmosphere), which are highly unlikely. The uncertainties for these p-factors are rather large compared to bestquality p-factors, which suggests that the data are not optimal. On the other hand, we cannot firmly exclude any residual bias in the parallaxes, for instance, or an effect related to the CSE of Cepheids. Values lower than p = 1 (not found in this subsample) would be physically possible if we were to consider that longperiod Cepheids (and therefore Cepheids with a large radius) have stronger dynamics and an intense atmospheric velocity gradient. Finally, no dependence on the period is clearly visible, in agreement with the conclusion of the study by Pilecki et al. (2018) using Cepheids in eclipsing binaries systems. Fitting a linear relation through the points in Fig. 8 gives the following relation between the period P and the p-factor: p = 1.251 ±0.008 -0.172 ±0.024 (log P -0.9), (4) with a high dispersion of 0.15. Considering only stars without effective temperatures and without interferometric measure-ments (blue points), we find the same dispersion of 0.14 around the same fit. Finally, stars with only an effective temperature (red points) show a scatter of 0.17 around the fit. The two stars with both effective temperature and good interferometry agree well with the slope of the fitted relation. There is no indication that one type of data is responsible for the large observed scatter. Additionally, it reinforces the robustness of the SPIPS method even for Cepheids with a limited dataset. Fitting a constant value through the points of Fig. 8 yields a projection factor of p = 1.26 ± 0.07 with a dispersion of 0.15, which is not significantly higher than the dispersion obtained for Eq. ( 4). Fixing this value to derive new distance estimates leads to a more dispersed PL relation:

M K s = -5.488 ±0.037 -3.515 ±0.120 (log P -0.9),

with σ = 0.22 mag, which is higher by ∼21% than the previous calibration (Eq. ( 1)). The quality criteria from Lindegren et al. (2021a) were verified for these stars, and we assume that biases due to a potential chromaticity effect (Breuval et al. 2020) are negligible in the EDR3, therefore this suggests that Gaia EDR3 parallaxes are sufficiently precise to let an intrinsic dispersion of the projection factor appear. Thus, the dispersion of the p-factors and the presence of values outside of the expected range suggest potential additional dependences of the P-p relation, or physics of the projection factor that is still not well understood. Pilecki et al. (2018) already suggested a dependence of the p-factor on other parameters than the period, such as the mass or radii. However, after some investigations, we did not find any correlation between the projection factor and these parameters or any other parameter, such as the mean effective temperature (Fig. From a theoretical point of view, Nardetto et al. (2011) predicted no correlation between the metallicity and the p-factor. Moreover, when we consider stars with the same period but extremely different p-factors in detail, no issue in the SPIPS modeling was highlighted. No correlation of radius and mass is clearly visible, in contrast to the suggestion by Pilecki et al. (2018).

A simplification made in the SPIPS algorithm is the parameterization of the infrared excess as a function of the wavelength with the assumption that there is no excess or deficit in optical bands. However, Hocdé et al. (2020) showed that this effect, physically understood as being due to a circumstellar envelope, can affect not only the infrared bands, but also optical ones.

Conclusions

We have presented the application of the SPIPS method to 63 Galactic Cepheids for which the most precise and complete dataset is available for the application of the PoP technique. This database covers almost 50 years of Cepheid observations, including multiband photometry, spectroscopic radial velocities, effective temperatures, and interferometric angular diameters. This modeling allowed us to derive new precise and consistent mean values of several parameters such as color excesses, period changes, angular diameters, effective temperatures, multiband mean apparent magnitudes, and the p-factor.

We established new calibrations of the period-luminosity and period-radius relations. We finally investigated the value and dependences of the projection factor: Gaia EDR3 parallaxes did not allow us to highlight a significant correlation between the p-factor and the period, but rather indicated that the p-factor is consistent with a constant value of p = 1.27 ± 0.06, with a significant dispersion of 0.15. This dispersion and the presence of unexpected p-factor values suggest that other important physical phenomena affect the PoP technique that have not yet been identified. Additionally, this study suggests that the period-pfactor relation may have an intrinsic width and/or may depend on many individual properties. However, its physical origin is still unknown and should be investigated in the future. We found no correlations between the p-factor or other parameters such as the mass, radius, effective temperature, or metallicity.

There are still several aspects to overcome before the pfactor is understood, and the very first is probably to wait for the final Gaia data release to obtain the best parallaxes possible in terms of precision and accuracy. In particular, improved Gaia distances in the next releases for Cepheids with many interferometric observations such as δ Cep, RS Pup, β Dor, or ζ Gem would permit us to obtain a better constraint on the p-factor. Another aspect to improve is the measurement of atmospheric velocity gradient using dedicated contribution functions of the line-forming regions. One of the best hopes is also related to the environment of Cepheids: recent studies appear to show that the circumstellar environment of Cepheids might not be static and may have some effects in the optical domain, and most probably in a different way, depending on the position of the Cepheids in the instability strip. This might explain the dispersion that we observe in the p-factors. This means that before we model the Cepheids and the p-factor in greater detail, we first need to understand the general scheme of the physics of the close circumstellar environnements of Cepheid through ongoing Cepheid observations in the NIR (MATISSE/VLTI) and opti-cal (CHARA/SPICA) domains. Moreover, parallel independent applications of the PoP technique would allow us to understand the physics of pulsating stars in more detail in order to conclude about the reliability of this method for the calibration of the extragalactic distance scale. ,2,3,4,5,9,10,11,12,13,14,15 7,16,17,18 1,3,4,5,9,10,11,12,13,14,30 7,12,16,17,19,36 ,2,3,4,5,9,11,13,14,15,24,30,43 7,16,17,18 ,3,4,5,9,10,11,12,13,14,39 7,12,17,19 8,25 21,22 U Sgr 1,2,3,4,5,11,13,14,15,24,30,32 7,16,17,36 8,25 21 BB Sgr 1,3,4,5,11,13,15,24,32 7,16 8 -XX Sgr 1,3,4,5,11,13 6,16 8,26 -W Sgr • 1,3,4,5,11,13,14,32,39 17,18,36 8,25 -A102, page 12 of 23 1,2,3,4,5,11,13,14,24,30,34,39,43 7,31 8,20 -Y Sgr • 1,3,4,5,11,13,14,15,32 6,16,17 8,26,33 21 X Sgr 1,3,4,5,11,13,14,15,32,42 16,17 8 -V0350 Sgr • 1,3,4,5,11,13 6,7,17, ,2,3,4,5,9,10,11,12,14 7,17,19,36 8,25 21,22 T Vul 1,3,4,5,9,10,11,12,13,14,32 12,17,19,36 8,2,5,13,14,24,30 7,17,31 8,20 21,22 SV Vul 1,2,3,4,5,9,10,11,12,13,14,24,30 7,16,17,19

Appendix A: References of the complete dataset