The Baade-Wesselink projection factor of RR Lyrae stars -- Calibration from OHP/SOPHIE spectroscopy and Gaia DR3 parallaxes

The application of the parallax of pulsation (PoP) technique to determine distances of pulsating stars implies the use of a scaling parameter, the projection factor (p-factor), required to transform disc-integrated radial velocities (RVs) into photospheric expansion velocities. The value of the p-factor is poorly known and debated. Most PoP applications assume a constant p-factor. However, it may actually depend on the physical parameters of each star. We aim to calibrate p-factors for RR Lyrae stars (RRLs) and compare them with classical Cepheids (CCs). Due to their higher surface gravity, RRLs have more compact atmospheres, and provide a valuable comparison with their supergiant siblings. We determined the p-factor of 17 RRLs using the SPIPS code, constrained by Gaia DR3 parallaxes, photometry, and new RVs from the OHP/SOPHIE spectrograph. We carefully examine the different steps of the PoP technique, particularly the method to determine RV from spectra using the classical cross-correlation function (CCF) approach. The method employed for RV extraction from the CCF has a strong impact on the p-factor, of up to 10%. However, this choice of method results in a global scaling of the p-factor, marginally affecting the scatter within the sample for a given method. Over our RRL sample, we find a mean value of $p = 1.248 \pm 0.022$ for RVs derived using a Gaussian fit of the CCF. There is no evidence for a different value of the p-factor of RRLs, although its distribution for RRLs appears significantly less scattered than that for CCs. The p-factor does not appear to depend in a simple way on fundamental stellar parameters. We argue that large-amplitude dynamical phenomena occurring in the atmospheres of RRLs and CCs during their pulsation affect the relative velocity of the spectral line-forming regions compared to the velocity of the photosphere.


Introduction
Pulsating stars, and in particular Cepheids and RR Lyrae, are essential indicators for the determination of distances in the local Universe.Precise calibration of the relations between their apparent luminosity, period of pulsation (Leavitt & Pickering 1912;Breuval et al. 2022), and metallicity is therefore particularly important.These calibrations are based on geometrical determinations of the distance using, for instance, parallaxes (see Groenewegen 2018, for example) or detached eclipsing binaries (Riess et al. 2019;Pietrzyński et al. 2019).Such empirical relations play a central role in the extragalactic distance ladder.Period-luminosity relations of Cepheids in the Milky Way, the Magellanic Clouds, and NGC 4258 are used to calibrate the lu-mated from the integration of the pulsation velocity curve of the photosphere of the star over its pulsation cycle.The resulting radius curve is then compared to the angular diameter curve, which can be either measured directly by interferometry for nearby stars, or predicted using the surface-brightness-colour relation (SBCR; see Kervella et al. 2004;Nardetto et al. 2023, for instance).The SBCR is generally calibrated using interferometric measurements.Unfortunately, it is not possible to directly measure the pulsational velocity by spectroscopy, but only the radial velocity (RV), which is an averaged quantity over the stellar disc.The ratio p = v puls /v rad between the pulsation velocity v puls and the RV v rad is called the projection factor, or p-factor.Considering a uniformly bright pulsating sphere, the projection on the line of sight of the spherical distribution of the pulsation velocity vector corresponds to p = 1.5.The p-factor also incorporates the effect of limb darkening, which results in a higher weight attributed to the RV in the centre of the star compared to the limb.This results in a lower value of p than the purely geometrical value, and typically around p ∼ 1.3.In addition, p is also affected by velocity gradients that may exist in the atmosphere of the star between the line-forming regions and the continuumforming photosphere.This aspect was studied theroretically in Nardetto et al. (2004) and with observations in Nardetto et al. (2007).Finally, the p-factor also depends on the method used to estimate the RV from the spectra.Such measurements need to be taken carefully, and well documented, as demonstrated by Anderson (2020) and Borgniet et al. (2019), but this is difficult because of the asymmetry of the lines in the spectrum of a pulsating star.
The intrinsic astrophysical complexity of the p-factor makes its calibration challenging.Using independent distance measurements (e.g. using Gaia parallaxes), several authors determined the p-factor of pulsating stars empirically (Storm et al. 2011;Ngeow et al. 2012;Breitfelder et al. 2016;Kervella et al. 2017;Nardetto et al. 2017;Trahin et al. 2021;Nardetto et al. 2023).The resulting large scatter of the determined p-factors led to a withdrawal of the PoP technique as a way to determine the distances of pulsating stars at the percent accuracy level.However, studying the p-factor is a powerful way to progress in our understanding of the dynamics of the atmosphere of pulsating stars.Among classical pulsators, Cepheids are known to exhibit relatively complex atmospheric dynamics (Fokin et al. 1996;Hocdé et al. 2020a), including mass loss (Neilson & Lester 2008) and circumstellar envelopes (Hocdé et al. 2020b;Gallenne et al. 2021).Conversely, RR Lyrae stars (RRLs), which present intrinsic differences in their atmospheres when compared to Cepheids (e.g. they are physically smaller than Cepheids, with a higher surface gravity, a more compact atmosphere, and a higher temperature), make for a valuable comparison in order to characterise their p-factor.This is the motivation behind the present study, which is intended to expand the investigation into classical Cepheids presented by Trahin et al. (2021) towards the lowermass RRLs.
In Sect.2, we present the new RV data that we collected with the OHP/SOPHIE spectrograph, as well as the sources of ancillary literature data (photometry and effective temperature).In Sect.3, we detail the methodology we used to extract RV measurements from our SOPHIE spectra.Section 4 is dedicated to the analysis of the pulsation of 17 RRLs using the SPIPS modelling tool.In Sect.5, we summarise our results on the p-factor, as well as those on the physical properties of the modelled stars.Finally, we conclude this study in Sect.6.

Spectroscopic data
We focus our study on a selection of 23 fundamental RRLs observable from the Northern hemisphere; these are listed in Table 1.We obtained new spectra for these RRLs using the SO-PHIE spectrograph installed at the 193cm telescope of the Observatoire de Haute-Provence in southern France.SOPHIE is a cross-dispersed échelle spectrograph (Perruchot et al. 2011) with excellent stability.This instrument typically provides a signal-tonoise ratio (S/N) of 40 within 10 minutes on the brightest RRLs in our sample.This relatively short exposure time is important, as 10 minutes already corresponds to 1-2% of the pulsation cycle of a RRL with a period of 10 hours.Longer exposure times would potentially result in smearing of the RV curve and underestimation of its amplitude.Our observations were divided into four runs of four or five nights to survey our RRL sample.We used the high resolution mode of SOPHIE, with the slow readout mode to limit the noise.During the night, we organised the sequence of stars according to the previous observations and the predicted phase of the star to obtain regular and complete coverage of the pulsation cycle.We were particularly careful to properly sample the phases between 0.8 and 1 when the stars rebound at their minimal radius, as the velocity changes quickly.For these phases, we obtained several consecutive exposures when possible.The automatic pipeline at OHP reduces the raw spectra and computes a cross-correlation function (CCF) with a specific mask.The RV is then determined as the mean value of a Gaussian profile adjusted to the CCF.For our RRLs, we initially chose an F0 type mask suitable for their typical effective temperature.In the present work, we recomputed the CCF using various methods that are detailed in Sect.3. Within the list of 26 stars originally selected, we obtained a complete phase coverage for 11 of them (the number of observations is greater than 15, the maximum gap in phase between two points is lower than 0.2, and the number of observations at the rebound is more than four).Six additional stars have a satisfactory phase coverage, although some non-critical phases are incompletely covered.For four stars, part of the RV curve is missing, and two stars have too few RV measurements for an independent analysis.Three stars were not observed at all and are not presented in the table.Finally, the two stars RS Boo and VX Her showed a Blazhko effect during the observations (as confirmed by the literature) and were removed from our sample.We also note the binarity of TU UMa (Liška et al. 2016), which is remarkable when comparing different RV datasets.For the stars with insufficient phase coverage, we completed our RV measurements with data from the literature.A limitation of this approach is that these additional RV measurements are potentially not determined from the spectra using the same method.

Photometric data
Photometric data were taken exclusively from the literature.To perform the modelling as described below, a large number of observations is required, and in different bands and photometric filters, in order to constrain the pulsating star model.Photometric light curves from Gaia (Gaia Collaboration et al. 2016, 2023) and Hipparcos (Detre & Szeidl 1973) were employed, as well as 2MASS magnitudes (Cutri et al. 2003).We also used the collection of photometric data assembled by Monson et al. (2017), and we added photometric measurements from the following sources: Barnes et al. (1988), Cacciari et al. (1987), Notes.The maximum gap and N rebound are two criteria that we used to select RV curves that are sufficiently well covered.The three last columns list the references for the RV, effective temperature, and photometric data.Star symbols ⋆ indicate the data sets that were effectively used in the model fitting.‡ refers to complete RV curves, † to satisfactorily covered RV curves, and * to incomplete RV curves.B indicates the stars that exhibit the Blazhko effect.The absence of any symbol indicates that there are too few observations. References.
(1) Barnes et al. (1988) Cacciari et al. (1989), Carney & Latham (1984), Fernley et al. (1993), Fernley et al. (1989), Fernley et al. (1990), Jones et al. (1987a), Jones et al. (1987b), Jones et al. (1988), Liu & Janes (1989), Liu & Janes (1990), Meylan et al. (1986), Skillen et al. (1989), Skillen et al. (1993a), andSkillen et al. (1993b).Photometry in the bluer filters (U, B, and Gaia B P bands) is generally not included in the model fit, as the stellar atmosphere models that we are using are insufficiently accurate for comparison with observations in this wavelength range because of some effects of the microturbulence (Casagrande & VandenBerg 2014).The Gaia G band (G_GAIA_GAIA3 in the models) covers a very broad wavelength range, which leads to difficulties for atmosphere models and for reddening corrections.For this reason, the G band light curves are not systematically used to constrain the atmosphere models along the pulsation cycle.Generally speaking, light curves in visible bands are much more common than in infrared filters in the literature.However, as it provides a strong constraint on the effective temperature and is also mildly sensitive to interstellar reddening, infrared photometry provides a particularly valuable constraint.This encouraged us to only consider part of the visible photometric measurements available in the literature in order to give them comparable weight in the fit to infrared measurements.The choice of data sets was based on phase coverage, agreement with other data sets, and dispersion or uncertainties.The data that were effectively used for the modelling of each RRL of our sample are marked with a star symbol ⋆ in Table 1.

Effective temperature measurements
Spectroscopic effective temperature estimates offer a valuable constraint on the atmosphere models of the star.This is due to the fact that these estimates are independent of the interstellar reddening of the star, and therefore help to mitigate the effects of the strong correlation that exists between the effective temperature and the colour excess E(B − V) of the star.We took such measurements from Andrievsky et al. (2018Andrievsky et al. ( , 2020Andrievsky et al. ( , 2021)), who used the Spectroscopy Made Easy tool, which compares a given spectrum with synthetic ones with different atmospheric parameters.Measurements from Takeda et al. (2006) and Takeda (2022) using the same approach were also used.APOGEE (Jönsson et al. 2020) and LAMOST (Cui et al. 2012) catalogues are of little use for this kind of work as they do not provide timeseries of effective temperature.APOGEE temperatures are measured from a mean spectrum and LAMOST do not provide precise epochs of measurements.However, we were able to check that these measurements are compatible with the obtained models for effective temperature.The number of collected effective temperature estimates is small for a single star.In order to balance the different observables in the fit of the model, we set a similar weight to each kind of observable (RVs, effective temperatures, magnitudes in the visible, and magnitudes in the infrared) in the global χ 2 .Even if photometric measurements are much more numerous and precise than RVs, the weight of the latter is properly balanced in the fit.A drawback of this choice is that the normalisation can potentially give disproportionate weight to the temperature measurements compared to the other observables.When temperature measurements were found to hamper the convergence of the model, leading to phase mismatching or unlikely light-curve shapes, we omitted them from the fitting process.However, models are still consistent with the measured temperatures.

Methods to determine RV and its impact on the p-factor
One key ingredient in the procedure to compute the p-factor of a pulsating star is the RV measurement.From the spectra acquired to the extraction of an RV value, some steps need to be carried out carefully.The most common way to do this is by fitting a Gaussian profile to the CCF.However, this implies a few choices.Firstly, the choice to compute the CCF with a specific mask or template.Secondly, there exists different definitions of the RV given the CCF; it can be defined as the centroid of the line profile as well as the mean of a Gaussian fit or as a parameter of another fitted function.Such differences in the measurements have a large impact on the p-factor as the angular diameter curve is compared to the integrated RV curve.One consequence is that RVs from different methods can lead to very different p-factors, which need to be quantified to identify the origin of the scatter of the p-factors in the pulsating-star population.In this section, we aim to quantitatively describe the impact of each of the aforementioned choices on the p-factor and to motivate one or the other choice, considering physical aspects but also the diversity observed in literature measurements.
To compare the various methods, we interpolate the RV curve using simple linear functions in order to avoid an impact of the choice of function used for the interpolation.This analysis is therefore only possible for really good phase coverage of the curve.We can then integrate ∆V = V − V mean with time on a pulsation cycle.This results in the integrated RV curve, of which we can measure the amplitude.This amplitude is directly correlated to the p-factor as it corresponds to the scaling of the integrated RV curve to the angular diameter curve.For each aspect studied in the following part, we present an example using a star with good phase coverage (especially at the rebound), along with its RV and integrated RV curves.A comparison with a reference is provided for the integrated RV curve, which is the method used afterwards.In this work, we used two different tools in addition to the SOPHIE pipeline.Firstly, iSpec (Blanco-Cuaresma et al. 2014;Blanco-Cuaresma 2019) was used to compute the CCF with different masks and templates, either standard or newly defined.In employing the CCF calculation method of this code, several strategies were tried, detailed in the subsequent section of the present paper.In parallel, the RaveSpan code (Pilecki et al. 2017) was also used because it allows a comparison of the CCF and broadening function (BF) approaches.Both functions show similar results but the BF has the advantage that it is easily computable from a poor S/N spectrum.Comparisons of those different tools and of the BF and CCF approaches are presented in Figs. 1 and 2 using a Gaussian profile fitting.The cross-correlation was done with the same template as that in Coelho et al. (2005), which is commonly used in RaveSpan, with T eff = 6000K, log g = 3.0, and [Fe/H] = −1.0.We compare the CCF and BF approaches using RaveSpan with the CCF computed by iSpec.The difference observed in the integrated RV curve is characterised by a difference of 1% and 2% in the amplitude respectively for the CCF and BF.

Masks and templates
We studied the choice between binary masks and synthetic templates for the CCF computation.We compared the standard SO-PHIE masks provided in the iSpec tools, and different templates derived from synthetic spectra.The comparison is not provided for spectral masks that are very far from the expected spectral type of RRLs, nor for templates with very different atmospheric properties.One can see the corresponding RV and integrated RV curves in Figs. 3 and 4, between two masks and two templates.Radial velocities were computed from a Gaussian fitting of the cross-correlation with an F0 and an A0 mask and with synthetic templates from Coelho et al. (2005) with properties similar to the masks.The corresponding error in the amplitude of the integrated RV curve is of 1%, and we conclude that it does not have a major impact on the p-factor.
Consequently, we have to keep in mind that the use of generic binary masks or templates, which contain a large number of spectral lines as they are usually used for SOPHIE observations, will probably average the gradient in RV in the atmosphere of the star, leading to a smoothed RV curve.
The choice between masks and templates leads us to the question of the line selection.The CCF is a profile computed from a large number of different spectral lines.However, it has been shown that not all lines are pertinent for the Baade-Wesselink method (Petterson et al. 2005).Indeed, different lines from different elements are forming at different levels on the atmosphere and do not offer a good comparison to angular diameter observations.We therefore studied the impact of the selection of the lines on the p-factor.

Line selection
The CCF is created by analyzing a set of specific lines from a spectrum.These lines are determined by the mask or template used to do the cross-correlation.Nardetto et al. (2007) reported that the amplitude of the RV curve increases with the depth of the line used to measure the RV.Borgniet et al. (2019) observed the effect of cross-correlating lines that are forming at different levels near the photosphere of Cepheids.To do so, these latter authors built masks that select either weak, medium, or deep lines, considering the depths of lines to be a proxy for the optical depth.They built their masks by selecting lines in a synthetic spectrum at 5250K and with solar metallicity.We followed the same approach for our RRLs, with new masks using a synthetic spectrum with a metallicity of -1 dex, a surface gravity of log g = 3.0, and a temperature of 6200K in order to be closer to RRLs.From such a normalised synthetic spectrum computed with the PHOENIX code (Husser et al. 2013), we applied some masks to extract the position of lines with specific ranges of depths.We defined three ranges: the deep lines for which the line depth is between 0.65 and 0.95; medium lines, with the depth between 0.45 and 0.65; and weak lines with a depth of higher than 0.2 and lower than 0.45.We tried to extend the limits from the definition in Borgniet et al. (2019) for the weak lines mask to include more lines in order to improve the S/N of the CCF.In addition, a mask and a template including all these lines are also constructed.Contrary to the masks, the template is not binary, as it is weighted using the depth of each line in the synthetic spectrum.For this reason it is quite similar to the mask for deep lines as these lines have a major impact on the CCF for the general mask.All these masks are constructed in the wavelength range from 4500Å to 6800Å1 .We note that the difference between the mask and the template for all lines is so small that it is not really pertinent to present both for each measurement.For each step in the analysis, we specify whether we use the binary mask or the template.
Figure 5 shows part of the synthetic spectrum as an example, with the delimitation of the different ranges of depths used and the resulting masks.It is important to note that in the line selection, we only keep unblended lines in order to have a shape that is as accurate as possible for the CCF, and not smoothed by some blended lines.The properties of the masks are summarised in Table 2. Nardetto et al. (2007) proposed that the CCF of the weak lines is the most pertinent for the PoP technique as it measures the velocity of the lines that were formed near the photosphere (those probed by interferometric and photometric measurements).Unfortunately, the S/Ns of some of our observations were not sufficient to have a well-defined peak in the CCF, some of them being too noisy.A satisfactory comparison for RV curves for different depths is therefore only provided for deep and medium lines in Fig. 7, as for weak lines, the phase coverage was too poor to obtain a good interpolation on the RV curve.
The RV curves for the three depth-restricted masks are shown in Fig. 6.All RVs are measured by fitting a Gaussian profile to the CCF.
We do not see a significant difference between RV measurements from masks on different line depths (1.3% and 0.3% in the amplitude of the integrated RV curve for deep lines and medium lines, respectively, with respect to all lines), which is similar to  (2007) showed that, for Cepheids, the difference in RV for different individual lines at different layers is stronger for long-period stars.We can extrapolate these observations and consider this effect to be negligible for RRLs.To further investigate the different CCFs using these different masks, we chose one star with high-quality spectra (S/N typically greater than 35 at 550nm) and compared the CCF using the masks for the three ranges of depth.Among our 23 stars, DX Del provides good spectra for this analysis at different phases, and so we can monitor how the shape evolves during the pulsation.For this star, the asymmetry of the CCF, which reverses between the phases of higher and lower velocity, is visible to the naked eye (Fig. 8).This observation is compatible with the conclusions of Nardetto et al. (2006) based on single lines, who pointed out that the asymmetry of the line variation curve has, in the case of Cepheids and at first order, a shape similar to that of the RV curve itself.Some characterisation of the CCF is given in Table 3.We computed the bisector of the line, spanning over 3σ (with σ being the standard deviation of the Gaussian) on both sides of the centre of the line to avoid the noise from the wings.
From this bisector, we obtain ⟨b⟩, the mean value of the bisector in this part of the line, and the amplitude in velocity ∆b.The CCF is also fitted by a bi-Gaussian profile with widths of σ 1 and σ 2 on the bluer and redder side (positive RV) of the CCF, respectively.This table shows the correlation between the RV measurement from a Gaussian fit and the mean of the bisector b.The absolute value ∆b and the signed value σ 2g 2 −σ 2g 1 are markers of the asymmetry of the CCF.The sign of the latter indicates which of the CCF wings is more extended: if this parameter is positive or negative, the larger wing extension is on the red or blue part of the CCF, respectively.Both evolve in the same way with the phase or the depth of the mask.We note that the asymmetry is systematically higher for medium and weak lines than for deep lines, which is consistent with previous observations (Anderson 2016).However, it is either of the same order or higher for medium lines with respect to weak lines.Borgniet et al. (2019) observed that weak lines are more sensitive to the asymmetry.

Measurement method on the CCF
The RV is then measured through the Doppler shift of the CCF.However, spectral lines in pulsating stars show a high level of asymmetry, which can bias the various possible estimators of this Doppler shift to a varying extent.For instance, a Gaussian fit of the CCF may be biased due to the asymmetry.An alternate approach is to fit a bi-Gaussian profile, defined as a Gaussian with different widths for the two sides of the profile.The bi-Gaussian profile reproduces the CCF of most of the data well (and therefore the line asymmetry derived from it will be properly determined), but does not produce a suitable RV for the projection factor (besides the additional information of the asymmetry that it provides, its mean value is almost consistent with the velocity Table 3. RV and asymmetry measurements on three spectra of DX Del.All velocities are in km.s −1 .Q is the quality factor as defined by Borgniet et al. (2019) associated to the minimum of the line, albeit less noisy) because it is highly sensitive to the line asymmetry generated, for example, by rotation.Instead, the centroid method, which is insensitive to rotation and width variation, is often preferred (Nardetto et al. 2006).
Here we present the effect of the choice of method on the RV curve (plotted in Fig. 9), and the resulting difference in the amplitude of the integrated RV curve (plotted in Fig. 10).Comparisons of RV curves using different techniques have been carried out for single lines (e.g.Nardetto et al. 2006), and on CCF for Cepheids (Anderson 2016) with the conclusion that there is a non-negligible difference in the RV curves when using one or the other method (i.e.Gaussian, centroid, bi-Gaussian etc.).We studied five definitions: the mean of a Gaussian fit, the mean of a bi-Gaussian fit, the mean of a Gaussian fit of the core of the CCF (defined as containing approximately 68% of the CCF absorption), the centroid, and the minimum using a fourth-order polynomial fit.For AB UMa, the difference in the amplitude of the curve ranges from 0.2% (centroid method) to 11% (bi-Gaussian fit) of the amplitude compared to the classic Gaussian fit method.Based on these differences, we estimate that the choice of the definition of the RV given a CCF can lead to a change in the p-factor of up to 10%.
As a summary of this analysis, Table 4 provides the major sources of variation in the amplitude of the integrated RV curve for each aspect studied, and for two different stars.We conclude that the main source of variation of the p-factor is caused by the choice of the CCF measurement method -rather than the choice of mask or template, the code used to compute the CCF, or even the line selection.It therefore appears that the choice of CCF measurement method can significantly affect the RV, and therefore also the p-factor.
For the remainder of our analysis and the determination of the p-factor, we used the mean of the Gaussian fit of the CCF, as this is the most common method but is also the most precise with a poor S/N (Anderson 2016).This allows us to directly compare our values of the p-factor with previous determinations from the literature.
All RV curves are plotted in Fig. 11  in grey (see details in Sect.4).For stars without a model, a dashed curve is plotted, which is the weighted mean of models of the most similar stars in terms of period (of those for which a model has been determined).This is interesting as it reveals the continuity of the shape of RV curves with changing period.

The SPIPS code
The Spectro-Photo-Interferometry for Pulsating Stars (SPIPS, Mérand et al. 2015) is a modelling code that simultaneously fits a broad range of observables of a pulsating star to determine its physical parameters (pulsation period, bolometric luminos- ity, distance or projection-factor, effective temperature, radius, reddening, etc.).SPIPS allows us to simultaneously fit different types of observations for a single star: RVs, spectroscopic effective temperatures, angular diameters from interferometry and photometric magnitudes in all bands and colours.Firstly, an interpolated curve is fitted to the RV and another for temperature data.Then, the other observables are predicted as a function of pulsation using the fitted parameters, and the photometry is modelled thanks to a grid of static ATLAS9 models (Castelli & Kurucz 2003) for the given temperature, metallicity, and surface gravity.The latter is estimated from the star's radius and the period-mass-radius relation from Bono et al. (2001).Such relations are calibrated for Cepheids but surface gravity does not have a large influence on the large-band photometry, with an impact on the p-factor of ≲ 1%.Filter transmissions and photometric zero points are used to compare with observations in a given set of filters.

The modelling strategy
To determine the p-factor for our RRLs, we adopted the distances from Bailer-Jones et al.Table 5. Measurements of RV from spectra with the SOPHIE spectrograph.We provide values from a Gaussian (RV g ) and a bi-Gaussian (RV 2g ) fit of the CCF of spectra with the computed binary mask containing all unblended lines and from a Gaussian fit of the CCF of spectra with the computed binary masks containing deep lines (RV deep ) or medium lines (RV medium ).The errors come from the fit of the Gaussian or bi-Gaussian fit and may be underestimated.This  (2021).The distances of the stars in our sample range between 536 and 1258 pc.For such nearby stars, errors on the parallax are quite small, and therefore one could also choose the distance as the inverse of the parallax.We checked that this choice does not affect the result.The impact is around 2% on the p-factor, which is quite similar to the uncertainty on this value.As we do not have spectroscopic temperature measurements for all stars, we adopted the colour excess E(B − V) predicted by the Stilism 3D maps using the same distances (Capi-tanio et al. 2017).Metallicities were set to the homogeneous values reported in Monson et al. (2017) (which include original values from Fernley et al. 1998a;Fernley & Barnes 1997;Fernley et al. 1998b), and we considered a typical uncertainty of ±0.15 dex.Because of the large number of parameters in the global fit, we individually checked the convergence of the optimisation to convincing solutions.As stars in our sample are relatively close and bright, they have been well observed and realistic initial parameters are easy to find.We chose to initialise the period, the  9).The bottom plot shows the difference from the more commonly used method, that is, the fit of a Gaussian profile to the CCF.gamma velocity, and the amplitude of the RV curve to the Gaia values found by Clementini et al. (2023).The starting value of the star's radius is defined according to the initial period thanks to the period-radius relation of Marconi et al. (2015).
It is important to note that the various data used are separated by up to several decades, as many photometric measurements were collected in the 1980s.Such an interval corresponds to tens of thousands of pulsation cycles.Therefore, the fit is highly sensitive to the period.This is why we only fit the period in the last step, allowing only very small variations.We also set the reference time to be in the average of all the data in order to avoid adding the errors as many times as the number of cycles between the older and the more recent observation.Finally, we occasionally considered the temporal change in the pulsation period as an additional parameter, depending on the goodness of the fit.
Here is the strategy we adopted to model the pulsation cycle of a typical star in our sample.First, we only fit the effective temperature curve with temperature measurements (if any) and photometry.The curve is fitted using Fourier series with as many harmonics as necessary to properly reproduce the data (between 6 and 8).Photometric colours are not used for the general fit as we favoured single-band photometric measurements.In a second step, we only fit the RV model curve to RV measurements.Because the phase coverage is generally not very dense, Fourier series are not suitable for RV curves as they introduce small oscillations.We therefore used cubic splines, with between six and eight nodes.The choice of the number of Fourier series or nodes of splines to use was made through a visual inspection of the model and of the individual and global χ 2 .Given these first guesses of both RV and temperature curves, we can then fit all the parameters but the period, before a finer and final tuning including the period determination.Metallicities, colour excesses, and distances are fixed during the whole procedure, but the impact of each of these parameters has been considered in the final p-factors and the associated uncertainties.
Because we were not able to obtain a good model with this procedure for some stars, we then adapted the code individually for each star to reach a good result in the fitting curves, that is, for RS Boo, AV Peg, and TU UMa.For these three stars, because of an observed phase shift between the model and some datasets that was not reduced by simply adding a period change, we performed a first model based on contemporaneous data and used the fitted parameters as initial parameters for a second run, including the other datasets.This helped the fitting process to converge.For RS Boo, the initial parameters were computed on data exclusively from Jones et al. (1988).For AV Peg we used the Monson et al. (2017) dataset together with the OHP/SOPHIE RV curve and Andrievsky et al. (2018) temperature for the first step.For TU UMa, the initial parameters were obtained using Monson et al. (2017) and Gaia photometry, the OHP/SOPHIE RV curve, and the temperature provided by Andrievsky et al. (2018) .In the particular case of RS Boo, its Blazhko effect prevented us from fitting one observable with various datasets, and so the RV curve from OHP/SOPHIE is not fitted.This leads to a higher uncertainty in the fit, especially for the period determination.
We found a good model for eight RRLs with new measurements taken at OHP using a Gaussian fit of the cross-correlation with a binary mask that includes all the identified unblended lines in a synthetic spectrum.We completed this sample with six additional stars from our list, using literature RVs that were offset to match our SOPHIE RV measurements (such an additive offset does not affect the p-factor).We also added three more southern stars using exclusively data from the literature.
Figures 12 and 13 present the result of the SPIPS modelling of two RRLs with new SOPHIE RV measurements.Others are provided in Appendix A 2 .The plot shows velocity, temperature, angular diameter, and photometric curves, both observed (points) and modelled (solid lines).The model curves shown in grey are computed based on the green curves whose parameters are optimised during the fitting process to fit the various observables simultaneously.For the photometric curves, the filter is indicated in blue.

Results
In total, 17 RRLs have been modelled using the SPIPS code.The properties of the modelled RRLs are summarised in Table 6.

p-factor
We compare our determinations of p-factors (listed in Table 6) with those of several previous Baade-Wesselink analysis -in terms of the ratio between the distance and the p-factor d/p.All these studies adopted a constant p-factor of between 1.3 and 1.36 to determine the distances of RRLs and their absolute magnitudes.A detailed comparison is provided in Table 7, and we summarise our main conclusions below.
-We note that our results for the ratio d/p are in good agreement with those of Fernley et al. (1989), Skillen et al. (1989), 2 Models are all available at the CDS Fernley et al. (1990), Skillen et al. (1993a) for X Ari, DX Del, SS Leo, VY Ser, and RV Oct, but not for BB Pup and WY Ant, and are also in agreement with the results of Liu & Janes (1990) for RR Cet, RX Eri, AV Peg, and TU UMa, but not for TT Lyn and SW Dra.-For the stars present in the studies by Jones et al. (1987aJones et al. ( ,b, 1988) (X Ari, SW Dra, VY Ser), our results do not agree in terms of the d/p ratio but the results for X Ari and VY Ser in the mentioned papers are also not consistent with the ones in Fernley et al. (1989Fernley et al. ( , 1990)), closer to the results of the present paper.-Similarly, our results for RR Cet and DX Del are comparable with those in Skillen et al. (1989) and Liu & Janes (1990), but differ from the results of Burki & Meylan (1986b,a).-For SW Dra, we find a different ratio d/p from all other authors (Jones et al. 1987b;Cacciari et al. 1989).-Finally, it is worth noting that, even if they agree within the error bars, ratios d/p derived from the SPIPS algorithm are almost systematically higher than those in the literature, leading to smaller p-factors.
The resulting p-factors of the 17 RRLs (defined for a Gaussian fit of the CCF) are plotted against the pulsation period in Fig. 14, encoded according to the origin of the RV measurements, and together with a linear fit and a constant fit.We find a mean value of p = 1.248 ± 0.022 In Fig. 15, we compare the dispersion of the RRL p-factors with that of the Cepheids studied by Trahin et al. (2021), who also used Gaia EDR3 parallaxes.The scatter for RRLs is smaller (7%) than for Cepheids (12%).It would be interesting to compare this to other classes of pulsating stars like δ Scuti (Nardetto et al. (2014) found higher p-factors with quite low dispersion of 2% for a sample of four stars) or type II Cepheids for instance (work by P.Wielgórski, in preparation).However, the scatter in p for RRLs remains too high to use the PoP technique to measure their distances with an accuracy of 1%.
Our p-factor distribution is compatible with a constant pfactor for all stars (RRLs and Cepheids) of p = 1.245 ± 0.015, which is uncorrelated with period.This constant value is compatible with previous results obtained with the SPIPS algorithm for Cepheids in the Magellanic Clouds and in the Milky Way (Gallenne et al. 2017), and is relatively close to the result of Breitfelder et al. (2016).
Comparison of the p-factor of the RRLs in our sample with some of their parameters (see Fig. 16) did not reveal any correlation, that is, with the effective temperature, the radius, the metallicity, or the amplitude of the RV curve.A more complex dependence is not excluded but further analysis would be necessary to investigate this properly.

Period-luminosity relation
From the SPIPS models, we also have access to the main parameters of the star, which allows us to determine period-radius and period-luminosity relations.While period-luminosity relations are the principal tool for distance determinations, they are not always defined for RRLs (Bhardwaj et al. 2023).Figure 17 shows this latter relation for the dereddened absolute magnitude in the Wesenheit index W JK = K S − 0.756(J − K S ) (using the reddening law from Fitzpatrick (1999) with R V = 3.1), taking into account a metallicity effect.We find the relation M W JK = −2.838log P + 0.126[Fe/H] − 1.161.The scatter around this relation is about 0.05 mag, which is smaller than for Cepheids (0.086 mag; see Breuval et al. 2022) Article number, page 10 of 26    TT_Lyn (P~ 0.597d) p=1.205 d=679.4pcE(B-V)=0.025Fig. 13.SPIPS model of TT Lyn with RV from SOPHIE spectra (Gaussian fit of the CCF using the binary mask built for RRLs with all unblended lines).Photometry measurements are from Monson et al. (2017) and Liu & Janes (1989) (29% of available photometric measurements).Table 6.Main properties of the best-fit models for 17 field RRLs.The distance (d), metallicity ([Fe/H]), and colour excess E(B − V) are fixed.We fit the period (P with uncertainty σ P ), if needed a change in period ( Ṗ), and the p-factor.From the RV and photometry curves, we determine the mean radius R, the mean effective temperature T eff , and the gamma velocity V γ .The three groups correspond to stars with SOPHIE RV data exclusively, stars with some additional literature RV data, and stars with literature RV data only, respectively.servables of a pulsating star, but is easy to reach thanks to SPIPS modelling.This relation is usually tighter than the periodluminosity relation.We show the relation in Fig. 18, the equation for which is log(R/R ⊙ ) = 0.770 ±0.003 log(P) + 0.9189 ±0.0002 .We compare with the relation found from non-linear convective models, where log(R/R ⊙ ) = 0.866 ±0.003 + 0.55 ±0.02 log(P) (Marconi et al. 2015).This relation is within the confidence interval for our sample.These coefficients are in agreement with those that we find within the range of period containing the majority of our sample.It is truly remarkable that the period-radius relation is relatively tight.This is a clue that the conversion from photometry to a radius curve is realistic.Hence, the dispersion of the p-factors we obtain may come from a spectroscopic effect related to shocks or other specific behaviours of the atmosphere of pulsating stars.blue and red edges of the instability strip for a metallicity similar to our RRLs, as predicted by Marconi et al. (2015).We note that the RRLs do not usually cool down beyond the red edge.On the contrary, they tend to leave the strip on the blue edge during their pulsation cycle.Nevertheless, their average positions are well confined within the instability strip.

Applicability to a Gaia sample
In this paragraph, we aim to test the possibility of performing a precise calibration of the p-factor, taking advantage of the Gaia photometry and RV curves available in the DR3.As also reported by Clementini et al. (2023), we find some significant differences between the RV curves from ground-based spectrographs and from Gaia.For some stars, the Gaia Radial Velocity Spectrograph (RVS) curves exhibit a bump before the maximum velocity.We show an example in Fig. 20.This difference is likely due to the fact that this measurement is done by cross-correlating a template on a narrow spectrum centred on the calcium triplet.These spectra can be affected by shock waves in a specific way.Evidence of such shocks mainly comes from studies of Balmer lines (Gillet & Crowe 1988;Fokin & Gillet 1997;Gillet et al. 2017), but our knowledge of its effect on the entirety of the spectrum is still incomplete.Hocdé et al. (2020a) already showed that a desynchronisation between Hα and the calcium triplet lines may occur for long-period Cepheids, for which the surface gravity is lower.The presence of such a bump will increase the amplitude of the integrated RV curve and therefore poses a problem for the PoP method.For the star SW Dra, we estimate that the amplitude of the radius curve will differ by 15%, and the p-factor will decrease to p = 1.04.
All RV curves from the 23 RRLs in our sample with their Gaia RVS and other RV measurements from the literature are plotted in Fig. 21.To enhance readability, only one additional literature dataset is included in the plot.This dataset is either the one used in SPIPS (if available) or the one offering the best phase coverage.At least 12 (V0341 Aql, RR Cet, W CVn, UY Cyg, SW Dra, RX Eri, WZ Hya, RR Leo, SS Leo, TT Lyn, UU Vir and BN Vul) of these stars quite clearly exhibit the same feature as SW Dra.For 6 other stars, we are not able to carry out the comparison, as either Gaia or ground-based measurements are missing.Only DX Del appears to really show a similar shape between all curves.14).Thick lines are the edges of the instability strip from models (Marconi et al. 2015) for a metallicity globally corresponding to RRLs.These predictions have been extended linearly to lower luminosities in semi-dotted lines.
It would be interesting to study the systematic uncertainties linked to the presence of this bump (if it is related to atmospheric parameters of the star or to some range of p-factors for instance), how it affects the curve, and why it appears for such spectral lines.This would require Gaia data release 4 (DR4), which will provide individual RVS spectra for each observing epoch.

Conclusions
In this paper, we present a calibration of the Baade-Wesselink projection factor for a selection of RRLs thanks to new measurements of RVs collected at OHP with the SOPHIE spectrograph.We computed the p-factor for 17 RRLs using Gaia DR3 parallaxes.We examined the influence of the RV measurement method on p by studying the cross-correlation of our spectra with various masks and templates.We find that the mask or template for CCF computation does not have a significant impact on p compared to the method employed to derive the RV from the CCF.Our study points out the importance of the precise documentation of how RVs are computed in the literature.However, for consistent RV series, the difference stemming from the methodology used to derive velocities should only shift the projection factor and not affect the scatter.We calibrated the pfactor using the SPIPS modelling tool.As for Cepheids (Trahin  Gaia -Ground Fig. 20.Comparison between the RV curves from the ground (Jones et al. 1987b, SOPHIE) and the Gaia series of SW Dra.The bottom plot is the difference between the Gaia RV series and the total ground RV curve.
et al. 2021), we see an important scatter in the results, which is greater than the uncertainties in Gaia parallaxes.This result has been observed by various authors who calibrated the projection factor, and does not seem to be uniquely correlated with any of the main physical parameters of the stars (period, temperature, metallicity, etc.).Interestingly, as shown in Fig. 15, the scatter of p seems to be the highest for Cepheids with short periods (below 10 days).However, whether it is a selection effect due to the size of the sample or a real phenomenon is uncertain.We note that we do not find any infrared excess for these RRLs.Although it can have an effect on the p-factors of Cepheids (Nardetto et al. 2023), the possible presence of a circumstellar envelope cannot explain the great diversity of the p-factors of pulsating stars.Another remaining uncertainty in the SPIPS analysis is the estimation of the colour excess E(B − V), which we determined from 3D maps of the interstellar medium.Although this value does not have a drastic importance for the p-factor, it is still poorly known, because of its degeneracy with the star's effective temperature.Meanwhile, RRLs in this study are generally nearby, and so the uncertainty is low for this parameter and does not represent a limitation to the present study.
The source of the scatter in the resulting p-factors is still poorly understood, but we argue that it is linked to a misunderstanding of the pulsation of the atmosphere of the star.There appears to be a missing ingredient between the pulsational velocity and the motion of line-forming regions, which may be linked to the observations of peculiar features in RVS from Gaia DR3.
Python package for Astronomy (Astropy Collaboration et al. 2013, 2018), the Numpy library (Harris et al. 2020), the Astroquery library (Ginsburg et al. 2019) and the Matplotlib graphics environment (Hunter 2007).We used the SIMBAD and VizieR databases and catalogue access tool at the CDS, Strasbourg (France), and NASA's Astrophysics Data System Bibliographic Services.The original description of the VizieR service was published in (Ochsenbein et al. 2000).This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

Fig. 1 .Fig. 2 .
Fig.1.Comparison between the iSpec and RaveSpan codes, and between the CCF and the BF approaches in analysing the RV curve for AB UMa.

Fig. 3 .
Fig. 3. Comparison between masks and templates in the RV curve for TT Lyn.

Fig. 4 .
Fig. 4. Comparison between masks and templates in the integrated RV curves from Fig.3for TT Lyn.The bottom plot shows the difference between each integrated RV curve and that derived from the F0 mask.

Fig. 5 .Fig. 6 .
Fig. 5. Portion of the synthetic spectrum used for the template construction and the corresponding built binary masks.The colour encodes the three different depth ranges, both for the background colours that define the depth region of each mask and for the plot of the binary masks themselves.

Fig. 7 .
Fig. 7. Comparison between medium, deep, and all lines.This figure shows the integrated RV curve (showed in Fig.6) from Fig.6for AV Peg.The bottom plot shows the difference with the mask keeping all lines with a weight associated to its depth.

Fig. 8 .
Fig.8.Evolution of the asymmetry of the CCF during the pulsation cycle of DX Del.The Gaia RV curve for this star is shown in black for global visualisation.The red markers represent the measurements from SOPHIE and the red stars locate three spectra with a high S/N where we can compute the CCF using the three masks with different depths, shown in the inserts.The dashed line shows the mean RV extracted from the Gaia catalogue.

Fig. 9 .Fig. 10 .
Fig. 9. Comparison between definitions of RV given a CCF.The crosscorrelation is computed with all lines weighted by depth mask.
Fig. 11.All RV curves obtained with the SOPHIE spectrograph.The curves have been normalised using Gaia parameters if an insufficient number of measurements was obtained.An offset in velocity is added to see all curves in a single plot, which is related to the period.

Fig. 12 .
Fig. 12. SPIPS model of AB UMa with RV from SOPHIE spectra (Gaussian fit of the CCF using the binary mask built for RRL with all unblended lines).Photometry measurements are fromMonson et al. (2017) (60% of available photometric measurements).

Fig. 14 .
Fig. 14. p-factors for RRLs from the SPIPS model, fitted with a linear function in black and a constant in red.Points are for stars with an RV curve from the OHP exclusively.Inverted triangles represent stars for which the RV curve from the literature has been rescaled to OHP measurements.Finally, squares are the three additionnal stars for which all data are from the literature.

Fig. 16 .Fig. 17 .
Fig. 16.Comparison between the p-factor of 17 RRLs and effective temperature, mean radius, metallicity, and amplitude of the RV curve.For each plot, the Spearman's coefficient is given to test the statistical probability of a monotonic relation between the p-factor and the other parameters.

Figure 19 Fig. 18 .
Figure 19 presents the pulsation cycle of all RRLs in the Hertzsprung-Russell (HR) diagram, together with the predicted

Fig. 19 .
Fig. 19.Hertzsprung-Russel diagram showing the pulsation cycle of the 17 modelled RRLs of our sample.Markers indicate the type of RV measurements used in the model (see legend of Fig.14).Thick lines are the edges of the instability strip from models(Marconi et al. 2015) for a metallicity globally corresponding to RRLs.These predictions have been extended linearly to lower luminosities in semi-dotted lines.
Article number, page 14 of 26 Garance Bras et al.: The Baade-Wesselink projection factor of RR Lyrae stars

Table 1 .
Selection of 23 RRLs in the northern hemisphere observed with the SOPHIE spectrograph at OHP, with their mean magnitude in band V (m V ), exposure time (∆t), and number of spectra acquired.

Table 2 .
Characterisation of the masks.We give the number of selected lines N, the mean depth ⟨d⟩, the mean wavelength ⟨λ⟩, and the mean width ⟨σ⟩.
. Uncertainties smaller than 0.01 km.s −1 in RV are not shown.

Table 4 .
Summary of the impact on the amplitude of the integrated RV curve when using different methods to determine RV for AB UMa and TT Lyn.
table is available in its entirety in the CDS.

Table 7 .
Comparison of results for the ratio d/p from the present study (second column) and previous determinations (last column).Stars indicate the fixed parameter for determinations from other studies.