Cepheid distances from the SpectroPhotoInterferometry of Pulsating Stars (SPIPS)
Application to the prototypes δ Cephei and η Aquilae^{⋆}
^{1}
European Southern Observatory,
Alonso de Córdova 3107, Casilla 19001,
19,
Santiago,
Chile
email:
amerand@eso.org
^{2}
LESIA (UMR 8109), Observatoire de Paris, PSL, CNRS, UPMC, Univ. ParisDiderot,
5 place
Jules Janssen, 92195
Meudon,
France
^{3} Unidad Mixta Internacional FrancoChilena de Astronomía,
CNRS/INSU (UMI 3386), France
^{4}
Departamento de Astronomía, Universidad de Chile, Camino El
Observatorio 1515, Las Condes,
1058
Santiago,
Chile
^{5}
Universidad de Concepción, Departamento de
Astronomía, Casilla
160C, Concepción,
Chile
^{6}
Center for High Angular Resolution Astronomy, Georgia State
University, PO Box
3965, Atlanta,
Georgia
303023965,
USA
^{7}
National Optical Astronomical Observatory,
Tucson,
AZ
85719,
USA
Received: 24 February 2015
Accepted: 24 August 2015
Context. The parallax of pulsation, and its implementations such as the BaadeWesselink method and the infrared surface brightness technique, is an elegant method to determine distances of pulsating stars in a quasigeometrical way. However, these classical implementations in general only use a subset of the available observational data.
Aims. Freedman & Madore (2010, ApJ, 719, 335) suggested a more physical approach in the implementation of the parallax of pulsation in order to treat all available data. We present a global and modelbased parallaxofpulsation method that enables including any type of observational data in a consistent model fit, the SpectroPhotoInterferometric modeling of Pulsating Stars (SPIPS).
Methods. We implemented a simple model consisting of a pulsating sphere with a varying effective temperature and a combination of atmospheric model grids to globally fit radial velocities, spectroscopic data, and interferometric angular diameters. We also parametrized (and adjusted) the reddening and the contribution of the circumstellar envelopes in the nearinfrared photometric and interferometric measurements.
Results. We show the successful application of the method to two stars: δ Cep and η Aql. The agreement of all data fitted by a single model confirms the validity of the method. Derived parameters are compatible with publish values, but with a higher level of confidence.
Conclusions. The SPIPS algorithm combines all the available observables (radial velocimetry, interferometry, and photometry) to estimate the physical parameters of the star (ratio distance/pfactor, T_{eff}, presence of infrared excess, color excess, etc). The statistical precision is improved (compared to other methods) thanks to the large number of data taken into account, the accuracy is improved by using consistent physical modeling and the reliability of the derived parameters is strengthened thanks to the redundancy in the data.
Key words: techniques: interferometric / stars: variables: Cepheids / stars: distances / circumstellar matter / methods: observational
FITS data files for each star are only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/584/A80
© ESO, 2015
1. Introduction
Cepheids are the backbone of the extragalactic distance ladder because their pulsation periods, which are easily measured observationally, correlate directly with their luminosities through Leavitt’s law (the periodluminosity relation, 38; 39). Thanks to their very high intrinsic brightness, they are visible in distant galaxies, as demonstrated for instance by Freedman et al. 22 or Riess et al. 62. They overlap with secondary, farreaching distance indicators, such as type Ia supernovae (SN Ia) or the TullyFischer relation, whose scales are anchored to Cepheid luminosities. Direct distance estimation of nearby Cepheids plays a crucial role in the calibration of Leavitt’s law and, as a consequence, of the extragalactic distance ladder used to observationally estimate the Hubble constant H_{0}e.g. 62. This importance has recently been reaffirmed by Suyu et al. 65: to the question “Are there compelling scientific reasons to obtain more precise and more accurate measurements of H_{0} than currently available?”, the authors answered “A measurement of the local value of H_{0} to one percent precision (i.e. random errors) and accuracy (i.e. systematic errors) would provide key new insights into fundamental physics questions and lead to potentially revolutionary discoveries”. These authors also recognized the role of the Cepheids and the problem of controlling the systematics in their distance determinations. An elegant and powerful method of directly measuring distances to Cepheids is the parallax of pulsation, also known as the BaadeWesselink (BW) method 2; 68, although Lindemann 40 suggested the same method eight years earlier, but has never been credited for it. In the BW technique, the variation of the angular diameter θ is compared to the variation of the linear radius (from the integration of the pulsation velocity V_{puls}). The distance d of the Cepheids is then obtained as the ratio between the linear and angular amplitudes, (1)The BW method uses in practice a combination of two quantities: (1) diskintegrated radial velocities, estimated from the changing Doppler shift of photospheric absorption lines; and (2) angular diameters, either derived from multicolor photometric measurements and surface brightness relations, or from interferometric measurements. One common property of these quantities is that they are derived from observations using models or some physical assumptions, therefore breaking the geometric nature of the parallax of pulsation. The BW method has demonstrated its capability to reach the onepercent statistical precision regime e.g., 43, and its true current limitation lies in the systematic uncertainties, which are probably between five and ten percent. Two problems directly contribute to these systematics: the projection factor p and the presence of circumstellar envelopes (CSEs). The projection factor is a multiplicative correction factor applied to the radial velocity derived from a spectroscopic absorptionline Doppler shift. This factor is used to unbias the spectroscopic measurement and estimate the true pulsation velocity. To first order, the radial velocity can be seen as the projection of the pulsation velocity, integrated over the surface of the star. Since the pulsation of Cepheids is radial, the limb of the star does not have a Doppler shift, whereas the point at the center of the apparent stellar disk has a maximum projected velocity toward the observer. Assuming a pulsation velocity of 1 km s^{1}, the measured diskintegrated radial velocity would be 1 /p = 1 / 1.5 = 0.67 km s^{1} for a uniformly bright sphere. p is lower than 1.5 for a limbdarkened star and more than 1.5 for a limbbrightened star. The pfactor is important because it biases the derived distance linearly: d/p is the unbiased measurement in the parallax of pulsation equation (Eq. (1)). For a long time, the adopted values of p were based on the linear periodpfactor relation established by Hindsley & Bell 29; 30: p = 1.39−0.03log P. This gives a value of p ≈ 1.36 for a typical tendayperiod Cepheid, which was the most commonly used value in the literature (see, e.g., 12). But with the first direct determination of the pfactor of 1.27 ± 0.06 for the star δ Cep 43, there has been a renewed interest in estimating the value of p. This work was based on the availability of a geometrical distance measurement, using the Fine Guidance Sensor (FGS) of the Hubble Space Telescope (HST). Since then, a dozen Cepheids have had their parallax measured directly in the same fashion 8. This allows us to estimate more values of p, and even calibrate it as a function of the pulsation period, using the infrared surface brightness (IRSB) version of the parallaxofpulsation method 64. Stars are limbdarkened in the spectral continuum and more darkened at shorter wavelength. However, it should be noted that stellar surfaces are slightly limbbrightened inside absorption lines. This leads to an apparent paradox: one would expect the pfactor to be 1.5 or higher, even though direct measurements instead lead to values of around 1.3.
To avoid the need of calibrating the projection factor, another approach is to include its contribution in the pulsation model. In their recent work, Gray & Stevenson 26 attempted to directly extract the pulsation velocity by using a simple geometric model of an absorption line deformed by the pulsation: the resulting pfactors they found for the radial velocity published using different measurement techniques vary from 1.30 to 1.38 for given star, leading to a systematic error of 6% on the parallax of pulsation distances. Again, this value is for a given star and results from the various datareduction techniques (e.g., bisector, crosscorrelation) used to extract the radial velocity from spectra 49. Another potential source of bias is the presence of circumstellar envelopes, which have been discovered and studied in the infrared by Kervella et al. 33, Mérand et al. 44; 45, Kervella et al. 34, and Gallenne et al. 23; 24; 25. In the context of the parallax of pulsation, these envelopes affect the infrared apparent brightness of the star from the Kband (2 μm) and longward of this. They also bias the angular diameters measured by infrared longbaseline interferometry. The geometry of the CSE seems to be almost universal (33; 44; 45) and to vary only in intensity. Even in the Gaia era, when a few hundred Galactic Cepheids will have their distance measured accurately, the parallax of pulsation will still be a invaluable tool for distance investigation. One might think, for instance, of studying the Large Magellanic Cloud Cepheids using this technique. In addition, it should be noted that the parallax of pulsation will remain an important tool for studying the physics of Cepheids: Gaia providing the distances, the BW studies of Galactic Cepheids will investigate the physics which it relies on.
2. Integrated method
2.1. Motivations
This work is the natural evolution of the method suggested by Barnes & Evans 3 to estimate the angular diameter from photometry. The generalization of the idea was proposed by Freedman & Madore 21 to provide a better physical basis for the parallax of pulsation and to call for taking into account all possible observables. They proposed to use a universal surface brightness to compute magnitudes, based on the following formula (for example, for band B): (2)where θ is the Rosseland angular diameter, T_{eff} the effective temperature, E(B − V) the color excess, B_{0} and C_{B} a set of parameters describing the surface brightness relation, and A_{B} the bandpassdependent reddening coefficient. This method has the disadvantage of requiring a calibration of B_{0} and C_{B}, and, more important, assumes a dependency of the surface brightness (here, a linear relation in effective temperature). These relations were recently calibrated by Pejcha & Kochanek 55 by analyzing thousands of measurements for dozens of Cepheids. We propose to use a different method that is unique thanks to a combination of two things:

We propose a “fit all at once” method (for a given star), which takes into account all the observables and fit all the parameters. This has the advantage of offering the best statistical accuracy and confidence in the result. Usually, BW methods are implemented by steps: first a radial velocity function is fitted analytically, then it is integrated, and finally compared to the angular diameter measurements to derive the distance. Unless treated properly (using a bootstrapping method, for example), this leads to an underestimation of the uncertainty of the final distance, unless the uncertainties on prior steps of the methods are propagated properly (e.g., the uncertainty on the radial velocity Fourier fit).

We try, as much as possible, to physically model the observables. For example, we propose synthesizing photometry based on atmospheric models and using calibrated bandpass filters, instead of using analytical surface brightness relations linear in color (such as V − K), which we know are not observationally linear, see for example Kervella et al. 31.
This approach also offers the potential of investigating, for example, why, in the case of δ Cep, the interferometric angular diameters of Mérand et al. 43 and the angular diameters derived by IRSB by Ngeow et al. 54 seem to systematically disagree by about 4%. A global method should be able to provide an answer to this contradiction. Another advantage of such a method is also to relax the constraint of uniform phase coverage to a certain extent; this was previously recognized by Freedman & Madore 21.
It is remarkable that global methods using physicsbased models are quite widespread in the field of determining fundamental parameters of eclipsing binaries. Implementations such as PHOEBE^{1}61 or ROCHE 59 use the same philosophy as we mentioned above. As a first path to implement such a method for Cepheids (this work), we developed a global approach for deriving fundamental parameters of the eclipsing binary δ Vel 46, which we successfully checked against the ROCHE model of the same system 60.
2.2. Description of the model
We assumed that Cepheids are radially pulsating spheres, with perfect cycletocycle repetition of their physical properties. The pulsation velocity and the effective temperature as a function of phase are described by periodic functions of the pulsation phase φ, interpolated using splines or Fourier series. Mérand et al. 43 showed that periodic spline functions often offer a better description of the pulsation of Cepheids than do Fourier series, since Cepheids often exhibit pulsation velocity variations that are very different from a simple sinusoidal wave. This requires many Fourier harmonics to describe the pulsation profile properly. Additionally, Fourier series fits are very sensitive to poor phase coverage and tend to introduce nonphysical oscillations. This means that Fourier decomposition requires a very uniform and dense phase coverage, which is not always available. However, Fourier series offer a good numerical stability, which is not always the case for a spline with freefloating nodes. In practice, we implemented both methods to allow for more flexibility. By default, Fourier series are used because they allow quicker computation and certain numerical convergence. We then switched to splines and kept this option if the goodness of fit was improved. Another important assumption was that Cepheid photospheres can be approximated by hydrostatic models in terms of energy distribution and centertolimb darkening. We used the set of astrophysical constants recently recommended by Harmanec & Prša 28.
Atmospheric models:
to compute synthetic photometry, we used ATLAS9 atmospheric models^{2}, with solar metallicity and a standard turbulent velocity of 2 km s^{1}. The effect of metallicity on the magnitudes is very weak, as noted by Casagrande & VandenBerg 14. We used a grid of models spaced by 250 K in effective temperatures and by 0.5 in log g. In practice, for each photometric bandpass, we reduced the models to a grid of magnitudes computed for an angular diameter of 1 mas. We then modeled the photometry by using the formula (here in B band) (3)This equation is similar to Eq. (2), except that the linear surface brightness relation is replaced by a grid of interpolated values B_{θ = 1 mas}, which is a function of the model: T_{eff} and log g. T_{eff}(φ) is fitted to the data (using either splines or Fourier series). On the other hand, log g is deduced from the parameters of the model: the mass of the star is assumed using the periodradiusmass relation of Bono et al. 11, and the linear radius is known internally in the model. The sensitivity of the M_{θ = 1 mas} to the gravity is, in any case, very low: this means that the choice of mass for the model is quite unimportant. As noted by Casagrande & VandenBerg 14, atmospheric models are poorly suited for reproducing synthetic photometry bluer than the B band, hence we limit our modeling to a range of 0.4 μm (B band) to about 2.5 μm (K band): the data presented here used the Johnson system in the visible (B and V bands), as well as the Walraven system (B and V band) and the CTIO system in the nearinfrared (J, H, and K bands).
Photometric bandpasses and zeropoints:
the photometric magnitudes were computed for each model of the grid, using bandpasses and zeropoints from the Spanish Virtual Observatory (SVO) database^{3} and the Asiago Database on Photometric System^{4}48 for the Walraven systems. Note that in the case of Walraven, we multiplied all the magnitudes by −2.5 since this unusual system expresses magnitude as the logarithm of the flux, without using the conventional −2.5 multiplicative factor. This allows for a uniform numerical treatment of all the photometric measurements. For the zero points, we chose the filters in the SVO that were recently calibrated by Mann & von Braun (see Table 1).
Adopted filters and zero points.
Reddening:
we parametrized the interstellar reddening using the B − V color excess, E(B − V), and the reddening law from Fitzpatrick 19, taken for Rv = 3.1. Because the correction depends on the spectrum of the observed object, we computed all our reddening corrections using a template spectra for actual effective temperature at the phase at which the photometric observations were made. Reddening values for T_{eff} = 4500 K, 5500 K, and 6500 K are listed in Table 2 for the various photometric systems we used. This is significantly different from traditional BW implementation. Reddening correction factors R_{λ} are usually computed for Vega, a star much hotter than the Cepheids. For example, Fouqué et al. 20 quotes R_{V} (i.e., for the V band) values between 3.10 and 3.30 and adopted a value of 3.23. As seen in our Table 2, our value for V_{GCPD} (Johnson) ranges from 3.00 to 3.05 between T_{eff} = 4500 K to T_{eff} = 6500 K (it would be 3.1 for T_{eff} = 10 000 K). We note that the effect of our choice of computation of the reddening is most notably different for blue filters and makes the least difference for the nearinfrared Kband. Our choice of Rv = 3.1 is mostly based on consensus and does not play a important role in the result: as far as we are concerned, the degeneracy is onetoone between the reddening law Rv and the color excess E(B − V). In other words, changing the fixed value of Rv changes the fitted value of E(B − V) and maintains the other parameters of the fit within their fitted values.
Subsets of magnitudes for θ = 1 mas and reddening law (for Rv = 3.1) for 3 values T_{eff} and logg = 1.5.
Centertolimb darkening:
the effect of the centertolimb darkening (CLD) needs to be taken into account to properly interpret interferometric angular diameters. Interferometers do not measure diameters directly, they measure visibilities, which need to be modeled in order to estimate an angular diameter. This is easiest to do using a uniform disk (UD) model. However, the derived diameter is not the true stellar diameter. Many authors have published tables of diameter corrections UD/LD, but we found that none are satisfactory, for the simple reason that the UD/LD correction depends on the spatial frequency at which the observations were made, because of the slight difference between UD and LD visibility profiles. For this reason we computed our own θ_{UD}/θ_{Ross.} corrections.
The truly interesting radius in our case is the bolometric radius, which almost matches the Rosseland value (where the average optical depth is 1). The Rosseland radius is the one that enters in the identity 6. In the context of this work, we used a grid of photospheric models tabulated in effective temperature: this is why the apparent Rosseland diameter (θ_{Ross}) is the one that allows to compute accurate synthetic photometry.
We did not use ATLAS models for our own CLD correction because these models are planeparallel and cannot produce accurate CLD profiles. Instead, we used grids of SATLAS models in the Cepheid range (51). The actual CLD profiles are available in the Vizier database (52, via FTP^{5}). We extracted the radial intensity profile I(r), which was converted to a visibility profile using a Haenkel transform, for various spatial frequencies (expressed as x = πBθ/λ, where B is the baseline in meters, θ the angular diameter in radian and λ the wavelength in meters). For each spatial frequency, we scaled the spatial frequency of a uniform disk visibility profile to match the synthetic profile: the scaling factor was the ratio θ_{UD}/θ_{Ross}. An example is shown in Fig. 1. We note that spherical models, tabulated as I(μ) (where ), do not have their limb for r = 1, in contrast to planeparallel models. This is because for spherical models, r = 1 is the outer boundary of the model (defined as the optical depth in the case of SATLAS, 50) and does not correspond to the Rosseland radius. We used a separate tabulation of R_{Rosseland}/R_{outer} extracted from the grid of SATLAS models (H. Neilson, priv. comm.). The mathematical justification of the equivalence of the scaling in r in the intensity profile and scaling the visibility curve to estimate the unbiased Rosseland angular diameter is a fundamental property of the Fourier transform: V [ I(a × r),Bθ_{LD} ] = V [ I(r),Bθ_{LD}/a ] = V [ I(r),Bθ_{Ross.} ], where B is the baseline and a = 1 /r_{Ross.}
We note that our results notably depart from those of Neilson et al. 53 for two reasons: 1) we took the radius of the star as the Rosseland radius, not the outer layer of the SATLAS model (defined as θ_{LD} by 53); and 2) our θ_{UD}/θ_{Ross} is a function of angular diameter and baseline. Overall, we found our values of θ_{UD}/θ_{Ross} to be higher than those published in Neilson et al. 53.
A limitation of our approach is that we used hydrostatic atmospheric models to compute our UD/Rosseland correction. This is not the latest way, since Marengo et al. 42 have used updated models to take into account nonhydrostatic effects. These authors found that the UD/Rosseland correction is, on average, comparable with the hydrostatic values and that the variation of the correction, due to the pulsation, is very small: about 0.3% in the nearinfrared and up to 1.5% in the visible. This translates more or less into the same bias in d/p bias. Since we mostly used nearinfrared optical interferometric data, the bias from our choice of using hydrostatic models is, to the best of our knowledge, only about 0.3%, at most. Moreover, there are no published grids of hydrodynamic models.
Fig. 1 Example of deriving the interferometric correction factor θ_{UD}/θ_{Ross.} for SATLAS model T_{eff} = 6000 K, log g = 1.5 and M = 10 M_{⊙}. Left: radial intensity profile, close to the limb (± 1%), for various bands; upper right: corresponding visibility functions as a function of the dimensionless spatial frequency x = πBθ/λ; lower right: corresponding factors θ_{UD}/θ_{Ross.} for each band as a function of x. 

Open with DEXTER 
Circumstellar envelopes:
the CSEs have two observational effects. The first one is on the near infrared photometric measurements, which are potentially biased for wavelengths in the K band (2.2 μm) and redder. The second effect is on the interferometric angular diameters. Kervella et al. 33 and Mérand et al. 44 showed that the fringe visibility as a function of the baseline length departs from the classical function of a limbdarkened star. In the case of the CSE, the bias on the measurements depends on the baselines and angular diameter. The approach we adopted was to use a grid of models using the parametrization reported by Perrin et al. 58, allowing the tabulation of the angular diameter bias as a function of infrared excess. Biases (θ_{observed}/θ_{real}) for different strengths of CSEs are shown in Fig. 2. We also allowed for an excess in H band, since these two bands are relatively close in wavelengths and it is hard to imagine that the CSEs produce Kband excess and no Hband excess. If no H excess is given as an input parameter, we chose to consider an Hband excess twice as low as the K band excess. The numerical process is very similar to the one we described for the limbdarkening correction: we synthesized the visibilities of a limbdarkened disk surrounded by the CSE, with the relevant observational parameters, and we fitted a uniform disk model to estimate the bias. This is numerically costly, but it is the only accurate way to estimate the bias.
Fig. 2 Kband interferometric angular diameter bias (observed/real) due to the CSE as a function of the dimensionless spatial frequency. 

Open with DEXTER 
Fig. 3 δ Cep data fit. Various panels show pulsation and radial velocities with spline model and residuals (panel a)); angular diameters and residuals, with the baseline colorcoded for the data and CSEbiased model – as a dash line, based on the model shown in Fig. 2 – (panel b)); effective temperatures (panel c)); photometric measurements and models (panels d) to k)) for different photometric bands or colors. Typical error bars are shown on the right side of the plot, below the reduced χ^{2} values. 

Open with DEXTER 
2.3. Fitting strategy
We used a standard χ^{2} minimization, (4)where O_{i} is the ith observations, e_{i} its associated error, and M_{i} the prediction from the model. The strategy to compute the overall χ^{2}, for all observations, necessitates some care. A normal χ^{2} would weight each measurement by its error bar. However, when we mix various observables, those that are present in large numbers are favored compared to scarce ones. A more general approach is to compute a χ^{2} by computing the final χ^{2} as the average of χ^{2} computed for each observable: (5)This is to ensure that each group G_{j} of observables contributes equally to the final likelihood estimation: for example, there are usually many more photometric observations than radial velocity or interferometric diameters. We used a LevenbergMarquardt (LM) leastsquares fit based on SciPy ^{6} scipy.optimize.leastsq . Using the total χ^{2} would have given more importance to data in highest numbers. Contrary to the approach taken by Pejcha & Kochanek 55, we did not fit the zero points of photometric systems, so we do not suffer degeneracy. After we found the best fit, we estimated the uncertainties in the derived parameters by using the covariance matrix around the bestfit solution.
Another aspect of the fitting process is the phasing of the data. It is known that Cepheids are not perfectly stable pulsators. For example, the slow (compared to the pulsation time) evolution of the star’s interior leads to a firstorder period change. The amount of linear change is an indicator of the evolutionary stage of the Cepheids and can be computed theoretically (see, for example, 18). We allowed the period to change linearly in our model.
3. Prototypical stars
Note that the observational data, and best fit model are available as FITS tables at the CDS.
3.1. δ Cep
δ Cep is the prototypical Cepheid and has been observed extensively, in particular by optical interferometer. We took the photometry from Moffett & Barnes 47, Barnes et al. 4, Kiss 35, Berdnikov 9 and Engle et al. 16. We also added photometric observations from Tycho and Hipparcos from van Leeuwen et al. 66 and ESA 17. We took the crosscorrelation radial velocities from Bersier et al. 10 and Storm et al. 63. The angular diameters are the ones published in Mérand et al. 43 and Mérand et al. 44. In addition, to properly interpolate the photospheric models, we adopted a metallically of [Fe/H] = 0.06, based on Andrievsky et al. 1. We note that the metallicity has a very weak effect on surface brightness values and is undetectable with our data set.
For the χ^{2} averaging, we used four groups of observables: radial velocities (91 measurements) angular diameters (67 measurements), photometric magnitudes (483 measurements), and colors (421 measurements). Error bars for each of these groups were multiplied by ~0.59, ~0.50, ~1.26, and ~1.35, respectively. We show the fit in Fig. 3, and the most important parameters are listed in Table 3.
Parameters of the δ Cep fit.
It is interesting to compare the result we obtain here with that of our previous study, which did not include photometry (43). The value of the pfactor is very similar: Using only the radial velocities and angular diameters reported by Bersier et al. 10, we found p = 1.27 ± 0.01. The uncertainty was smaller since we took into account correlations in interferometric error bars (using the formalism of 57), which we do not yet have implemented in our current SPIPS fitting algorithm. The actual pfactor uncertainty should take into account the distance uncertainty (0.050), however, which is much larger that the statistical uncertainty (0.020).
The CSE is noticeable in the interferometric data as a bias affecting the angular diameter measured at the shortest baselines. Mérand et al. 44 did not fit the excess, but rather compared the fit using a simple star model to a fit using the model we fitted on another Cepheid (Polaris), for which we had extended baseline coverage. At the time, we used a 1.5% excess (0.016 mag). In the case of SPIPS, we have photometric data that allow anchoring the model and allow using the CSE contribution as a free parameter. Thanks to this, we confirmed the infrared excess and estimated it to be 0.025 ± 0.002 mag in K band. We also let the H excess free to vary to fit the photometry and found it to be 0.018 ± 0.004. This latter is solely based on the photometric measurements.
The good agreement with all the observables is remarkable and increases our confidence in the method. In particular, our SPIPS modeling is able to combine all data and does not show apparent discrepancies between optical interferometry and IRSB, such as noted by Ngeow et al. 54. Admittedly, we added the complexity of having an infrared excess, which probably explains the discrepancy (which 54 did not take into account). One could argue that the Kband magnitudes do not agree the best agree in our fit (Fig. 3, panel “h”). We also performed a fit using only photometric measurements (omitting our interferometric measurements) and found the pfactor to be 1.29 ± 0.06, which, apart from the poorer statistical uncertainty, agrees perfectly well with our fit using optical interferometry. The K excess was also let free in the photometric fit, and its value was found to be 0.010 ± 0.004 magnitude.
Additionally, the period change (−0.07 ± 0.03 s/yr) is found to agree well with the recent estimate by Engle et al. 16, even though these authors have a much greater accuracy (−0.1006 ± 0.0002 s/yr).
3.2. η Aql
η Aql is another important prototypical Cepheid because of its proximity (and hence large apparent size), which makes it accessible to optical interferometry. We observed η Aql in July 2006, using the FLUOR instrument (15) at the CHARA Array. We used the same data reduction approach as in previous works, in particular the δ Cep data used in the previous section.
We took photometry from Welch et al. 67, Moffett & Barnes 47, Barnes et al. 4, Kiss 35, Berdnikov 9. Photometric measurements in the Walraven system were taken from Pel 56. We also added photometric observations from Tycho and Hipparcos from van Leeuwen et al. 66 and ESA 17. Radial velocities were taken from Barnes et al. 5 and Kiss 35. Finally, we also took additional angular diameter measurements: H band longbaseline measurements from Lane et al. 36 and shortbaseline Kband measurements from Kervella et al. 32. We adopted a metallically of [Fe/H] = 0.05 , based on Andrievsky et al. 1.
Fig. 4 η Aql fit. Various panels show pulsation and radial velocities with spline model and residuals (panel a)); angular diameters and residuals, with the baseline colorcoded for the data and CSEbiased model – as a dash line, based on the model shown in Fig. 2 – (panel b)); effective temperatures (panel c)); photometric measurements and model (panels d) to m)) for different photometric bands or colors. Typical error bars are shown on the right side of the plot, below the reduced χ^{2} values. 

Open with DEXTER 
Parameters of the η Aql fit.
The results of the fit are presented in Fig. 4 and Table 4. As for δ Cep, we applied a correction factor to the error bars to equally weight the four following groups: radial velocities (57 measurements, 0.5 factor), angular diameter (70 measurements, 0.55 factor), photometric magnitudes (377 measurements, 1.3 factor), and photometric colors (432 measurements, 1.35 factor).
We detect a slight H and Kband infrared excess (0.016 ± 0.003 and 0.018 ± 0.002, respectively). Like δ Cep, this is allowed by the combination of infrared photometry and infrared interferometric angular diameters.
Regarding the accuracy of E(B − V), 37 reported 0.126 and also quoted an older value of 0.143 (13), as well as 0.138 (metallicity corrected, computed by the software “BELRED”). Groenewegen 27 quoted 0.130 ± 0.009. Storm et al. 64 used 0.129. Our estimate is in this range, at 0.161 ± 0.005, on the redder side. The statistical uncertainty we obtain, ± 0.005, is underestimated because we did not take into account the fact that all photometric measurements in a same band and from a same source share a common error, namely the zero point and the photometric calibrators. If we perform a Jackknife resampling, removing one set of photometric measurements every time, the uncertainty on E(B − V) increases by a factor of 3, to ± 0.015.
Regarding the distance, η Aql appears in Table 5 of Groenewegen 27 with a distance of 261 ± 6 ± 7 pc for p = 1.321 (d/p = 198 ± 5 ± 4 pc). Storm et al. 64 determined a distance of 255 ± 5 pc using IRSB method, for p = 1.39 (d/p = 183 ± 4 pc). Using a subset of data we used, Lane et al. 36 obtained d = 320 ± 32 pc with p = 1.43 (d/p = 223 ± 22 pc). Our method gives a distance of 296 ± 5 pc (d/p = 228 ± 4 pc), which is not consistent with Storm et al. 64. We note that our uncertainty is on the same order as that of Storm et al. 64, and surprisingly, they used only radial velocity and twoband photometry. If we restrict ourselves to IRSB data (radial velocities and V, K photometry), our fit leads to ± 15 pc. Since we cannot fit E(B − V) (because of the degeneracy with T_{eff}), we should estimate the sensitivity of the distance estimate to change in E(B − V). We computed that decreasing E(B − V) by 0.05 leads to a distance 4 pc smaller. In other words, restricting our data set to the IRSB method leads to similar distances. The reason why we find an uncertainty in the estimated distance three times larger than Storm et al. 64 is the following: we suspect that since we fitted all parameters at once (radial velocity profile, T_{eff} profile, distance, etc.), our uncertainties are more realistic. If we keep our η Aql model and only use the IRSB dataset, and if we assume that we know everything in the model except for the distance and only adjust for this parameter, the uncertainty decreases to ± 5 pc, which is the claim of Storm et al. 64. In other words, our analysis of η Aql is a perfect example of why fitting all parameters at the same time provides more realistic uncertainties.
4. Conclusions
Our model makes many simplistic assumptions about Cepheids, most of which are known to be incorrect at a certain level. However, in the context of the parallaxofpulsation distance estimation, our approach is more complete than most (if not all) implementations that are variations of the BaadeWesselink method (BWM): 1) we include all possible observables, including redundant ones; and 2) we use observation modeling based on a physical model (as opposed to adhoc parameters, such as the surface brightness relations). Our implementation includes the traditional BWM, if one restricts the input data set. Using our modeling, we address some shortcomings of the BWM:

We adopted an approach of modeling the observables rather than using illdefined corrective factors. For example, we used modeled interferometric visibility profiles to compute the interferometric bias θ_{UD}/θ_{Ross.} whereas it is traditionally derived for brightness profile fits to analytical functions. We still make use of the projection factor, but we are working on a spectral synthesis modeling to allow us to use a consistent pulsation velocity estimation.

We used atmospheric models (ATLAS9 in our case) to compute synthetic photometry. This works very well, as proven by the agreement with interferometric angular diameters on our two prototypical stars. We note that the resulting surface brightness relation cannot be approximated by a linear function of the effective temperature (or color), as is done with a traditional implementation of the BWM. Because the BWM lacks redundancy in the dataset it uses, this shortcoming cannot be detected and propagates as a color bias on the distances.

Circumstellar envelopes (CSE) are consistently taken into account in the nearinfrared photometry and optical interferometric diameters.

Reddening is fitted from the data in a selfconsistent way. Conversely, BWM uses an E(B − V) that was determined for a certain reddening law and often applies it using another reddening law. Our method does not suffer from this bias.

Our approach permits very good phasing of data, even taken at different epochs. Not only does it improve the accuracy of the distance determination (because poorly phased data often have underestimated amplitude), it also allows us to study the period change of Cepheids.

Fitting all parameters at once realistically estimates the statistical uncertainties, as opposed to a method that would fit consecutive sets of parameters. For example, if the analytical radial velocity function is fitted first in an implementation of the BMW and then the analytical variations of angular diameters, followed by the distance alone as the ratio between the two, the uncertainty of the distance would not account for the other uncertainties and would likely be underestimated by a factor as large as 3.
All this should come as a warning to studies using only two bands: their distance (or pfactor) determinations probably have systematic errors that are hard to estimate without using a method like the one we have presented. Even then, their statistical uncertainties might very well be underestimated by a large factor. We applied the method to δ Cep and η Aql. For δ Cep we confirm our formerly published values for the pfactor of 1.28 ± 0.06, accounting for the uncertainty of the distance by Benedict et al. 7 of 274 ± 11 pc. For η Aql, we estimated its biased distance to be d/p = 228 ± 4 pc, leading to d = 296 ± 5 pc assuming p = 1.30. In both cases, our models reproduced all the available data (about a thousand observations in each case), in a selfconsistent way. In the near future, we will continue our work by systematically studying Cepheids for which large datasets are available.
Acknowledgments
We would like to thank the referee, Hilding Neilson, for his work that led to a much improved manuscript, as well as for providing additional insights to the use of SATLAS models described in the present work. This research has made use of the Spanish Virtual Observatory supported from the Spanish MEC through grant AyA200802156. This research has made use of the VizieR catalog access tool and SIMDAD database, operated at CDS, Strasbourg, France. A.G. acknowledges support from FONDECYT grant 3130361. P.K. and J.B. acknowledge financial support from the “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU, France, and the ECOS/Conicyt grant C13U01. The CHARA Array is funded by the National Science Foundation through NSF grants AST0908253 and AST1211129, and by the Georgia State University through the College of Arts and Sciences. STR acknowledges support by NASA through grant number HSTGO12610.001A from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS 526555.
References
 Andrievsky, S. M., Kovtyukh, V. V., Luck, R. E., et al. 2002, A&A, 381, 32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baade, W. 1926, Astron. Nachr., 228, 359 [NASA ADS] [CrossRef] [Google Scholar]
 Barnes, T. G., & Evans, D. S. 1976, MNRAS, 174, 489 [NASA ADS] [CrossRef] [Google Scholar]
 Barnes, III, T. G., Fernley, J. A., Frueh, M. L., et al. 1997, PASP, 109, 645 [NASA ADS] [CrossRef] [Google Scholar]
 Barnes, III, T. G., Jeffery, E. J., Montemayor, T. J., & Skillen, I. 2005, ApJS, 156, 227 [NASA ADS] [CrossRef] [Google Scholar]
 Baschek, B., Scholz, M., & Wehrse, R. 1991, A&A, 246, 374 [NASA ADS] [Google Scholar]
 Benedict, G. F., McArthur, B. E., Fredrick, L. W., et al. 2002, AJ, 124, 1695 [NASA ADS] [CrossRef] [Google Scholar]
 Benedict, G. F., McArthur, B. E., Feast, M. W., et al. 2007, AJ, 133, 1810 [NASA ADS] [CrossRef] [Google Scholar]
 Berdnikov, L. N. 2008, VizieR Online Data Catalog: II/285 [Google Scholar]
 Bersier, D., Burki, G., Mayor, M., & Duquennoy, A. 1994, A&AS, 108, 25 [NASA ADS] [Google Scholar]
 Bono, G., Gieren, W. P., Marconi, M., Fouqué, P., & Caputo, F. 2001, ApJ, 563, 319 [NASA ADS] [CrossRef] [Google Scholar]
 Burki, G., Schmidt, E. G., Arellano Ferro, A., et al. 1986, A&A, 168, 139 [NASA ADS] [Google Scholar]
 Caldwell, J. A. R., & Coulson, I. M. 1985, MNRAS, 212, 879 [NASA ADS] [CrossRef] [Google Scholar]
 Casagrande, L., & VandenBerg, D. A. 2014, MNRAS, 444, 392 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Coudé du Foresto, V., Borde, P. J., Merand, A., et al. 2003, in Interferometry for Optical Astronomy II, ed. W. A. Traub, SPIE Conf. Ser., 4838, 280 [Google Scholar]
 Engle, S. G., Guinan, E. F., Harper, G. M., Neilson, H. R., & Remage Evans, N. 2014, ApJ, 794, 80 [NASA ADS] [CrossRef] [Google Scholar]
 ESA 1997, The Hipparcos and Tycho catalogues. Astrometric and photometric star catalogues derived from the ESA Hipparcos Space Astrometry Mission, ESA SP, 1200 [Google Scholar]
 Fadeyev, Y. A. 2014, Astron. Lett., 40, 301 [NASA ADS] [CrossRef] [Google Scholar]
 Fitzpatrick, E. L. 1999, PASP, 111, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Fouqué, P., Arriagada, P., Storm, J., et al. 2007, A&A, 476, 73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Freedman, W. L., & Madore, B. F. 2010, ApJ, 719, 335 [NASA ADS] [CrossRef] [Google Scholar]
 Freedman, W. L., Madore, B. F., Gibson, B. K., et al. 2001, ApJ, 553, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Gallenne, A., Mérand, A., Kervella, P., & Girard, J. H. V. 2011, A&A, 527, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gallenne, A., Kervella, P., & Mérand, A. 2012, A&A, 538, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gallenne, A., Mérand, A., Kervella, P., et al. 2013, A&A, 558, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gray, D. F., & Stevenson, K. B. 2007, PASP, 119, 398 [NASA ADS] [CrossRef] [Google Scholar]
 Groenewegen, M. A. T. 2008, A&A, 488, 25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Harmanec, P., & Prša, A. 2011, PASP, 123, 976 [NASA ADS] [CrossRef] [Google Scholar]
 Hindsley, R., & Bell, R. A. 1986, PASP, 98, 881 [NASA ADS] [CrossRef] [Google Scholar]
 Hindsley, R. B., & Bell, R. A. 1989, ApJ, 341, 1004 [NASA ADS] [CrossRef] [Google Scholar]
 Kervella, P., Bersier, D., Mourard, D., et al. 2004, A&A, 428, 587 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kervella, P., Nardetto, N., Bersier, D., Mourard, D., & Coudé du Foresto, V. 2004, A&A, 416, 941 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kervella, P., Mérand, A., Perrin, G., & Coudé du Foresto, V. 2006, A&A, 448, 623 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kervella, P., Mérand, A., & Gallenne, A. 2009, A&A, 498, 425 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kiss, L. L. 1998, J. Astron. Data, 4, 3 [Google Scholar]
 Lane, B. F., CreechEakman, M. J., & Nordgren, T. E. 2002, ApJ, 573, 330 [NASA ADS] [CrossRef] [Google Scholar]
 Laney, C. D., & Caldwell, J. A. R. 2007, MNRAS, 377, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Leavitt, H. S. 1908, Annals of Harvard College Observatory, 60, 87 [NASA ADS] [Google Scholar]
 Leavitt, H. S., & Pickering, E. C. 1912, Harvard College Observatory Circular, 173, 1 [NASA ADS] [Google Scholar]
 Lindemann, F. A. 1918, MNRAS, 78, 639 [NASA ADS] [CrossRef] [Google Scholar]
 Mann, A. W., & von Braun, K. 2015, PASP, 127, 948 [NASA ADS] [CrossRef] [Google Scholar]
 Marengo, M., Karovska, M., Sasselov, D. D., et al. 2003, ApJ, 589, 968 [NASA ADS] [CrossRef] [Google Scholar]
 Mérand, A., Kervella, P., Coudé du Foresto, V., et al. 2005, A&A, 438, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mérand, A., Kervella, P., Coudé du Foresto, V., et al. 2006, A&A, 453, 155 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mérand, A., Aufdenberg, J. P., Kervella, P., et al. 2007, ApJ, 664, 1093 [NASA ADS] [CrossRef] [Google Scholar]
 Mérand, A., Kervella, P., Pribulla, T., et al. 2011, A&A, 532, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Moffett, T. J., & Barnes, III, T. G. 1984, ApJS, 55, 389 [NASA ADS] [CrossRef] [Google Scholar]
 Moro, D., & Munari, U. 2000, A&AS, 147, 361 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nardetto, N., Gieren, W., Kervella, P., et al. 2009, A&A, 502, 951 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Neilson, H. R., & Lester, J. B. 2008, A&A, 490, 807 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Neilson, H. R., & Lester, J. B. 2013, A&A, 554, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Neilson, H. R., & Lester, J. B. 2013, VizieR Online Data Catalog, 355, 49098 [Google Scholar]
 Neilson, H. R., Nardetto, N., Ngeow, C.C., Fouqué, P., & Storm, J. 2012, A&A, 541, A134 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ngeow, C.C., Neilson, H. R., Nardetto, N., & Marengo, M. 2012, A&A, 543, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pejcha, O., & Kochanek, C. S. 2012, ApJ, 748, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Pel, J. W. 1976, A&AS, 24, 413 [NASA ADS] [Google Scholar]
 Perrin, G. 2003, A&A, 400, 1173 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perrin, G., Ridgway, S. T., Verhoelst, T., et al. 2005, A&A, 436, 317 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pribulla, T. 2012, in IAU Symp. 282, eds. M. T. Richards, & I. Hubeny, 279 [Google Scholar]
 Pribulla, T., Merand, A., Kervella, P., et al. 2011, A&A, 528, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Prša, A., & Zwitter, T. 2005, ApJ, 628, 426 [NASA ADS] [CrossRef] [Google Scholar]
 Riess, A. G., Macri, L., Casertano, S., et al. 2011, ApJ, 730, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Storm, J., Carney, B. W., Gieren, W. P., et al. 2004, A&A, 415, 521 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Storm, J., Gieren, W., Fouqué, P., et al. 2011, A&A, 534, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Suyu, S. H., Treu, T., Blandford, R. D., et al. 2012, ArXiv eprints [arXiv:1202.4459] [Google Scholar]
 van Leeuwen, F., Evans, D. W., Grenon, M., et al. 1997, A&A, 323, L61 [NASA ADS] [Google Scholar]
 Welch, D. L., Wieland, F., McAlary, C. W., et al. 1984, ApJS, 54, 547 [NASA ADS] [CrossRef] [Google Scholar]
 Wesselink, A. J. 1946, Bull. Astron. Inst. Netherlands, 10, 91 [NASA ADS] [Google Scholar]
All Tables
Subsets of magnitudes for θ = 1 mas and reddening law (for Rv = 3.1) for 3 values T_{eff} and logg = 1.5.
All Figures
Fig. 1 Example of deriving the interferometric correction factor θ_{UD}/θ_{Ross.} for SATLAS model T_{eff} = 6000 K, log g = 1.5 and M = 10 M_{⊙}. Left: radial intensity profile, close to the limb (± 1%), for various bands; upper right: corresponding visibility functions as a function of the dimensionless spatial frequency x = πBθ/λ; lower right: corresponding factors θ_{UD}/θ_{Ross.} for each band as a function of x. 

Open with DEXTER  
In the text 
Fig. 2 Kband interferometric angular diameter bias (observed/real) due to the CSE as a function of the dimensionless spatial frequency. 

Open with DEXTER  
In the text 
Fig. 3 δ Cep data fit. Various panels show pulsation and radial velocities with spline model and residuals (panel a)); angular diameters and residuals, with the baseline colorcoded for the data and CSEbiased model – as a dash line, based on the model shown in Fig. 2 – (panel b)); effective temperatures (panel c)); photometric measurements and models (panels d) to k)) for different photometric bands or colors. Typical error bars are shown on the right side of the plot, below the reduced χ^{2} values. 

Open with DEXTER  
In the text 
Fig. 4 η Aql fit. Various panels show pulsation and radial velocities with spline model and residuals (panel a)); angular diameters and residuals, with the baseline colorcoded for the data and CSEbiased model – as a dash line, based on the model shown in Fig. 2 – (panel b)); effective temperatures (panel c)); photometric measurements and model (panels d) to m)) for different photometric bands or colors. Typical error bars are shown on the right side of the plot, below the reduced χ^{2} values. 

Open with DEXTER  
In the text 