A thin shell of ionized gas explaining the IR excess of classical Cepheids
V. Hocdé N. Nardetto E. Lagadec G. Niccolini A. Domiciano de Souza A. Mérand P. Kervella A. Gallenne M. Marengo B. Trahin W. Gieren G. Pietrzynski S. Borgniet L. Breuval B. Javanmardi


A thin shell of ionized gas explaining the IR excess of Cepheids
Hocdé et al. 


Université Cot̂e d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, France,

email : 
European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching, Germany
LESIA (UMR 8109), Observatoire de Paris, PSL, CNRS, UPMC, Univ. Paris-Diderot, 5 place Jules Janssen, 92195 Meudon, France  
European Southern Observatory, Alonso de Córdova 3107, Casilla 19001, Santiago 19, Chile 
Department of Physics and Astronomy, Iowa State University, Ames, IA 50011, USA
Nicolaus Copernicus Astronomical Centre, Polish Academy of Sciences, Bartycka 18, PL-00-716 Warszawa, Poland
Universidad de Concepción, Departamento de Astronomía, Casilla 160-C, Concepción, Chile



Despite observational evidences, InfraRed (IR) excess of classical Cepheids are seldom studied and poorly understood, but probably induces systematics on the Period-Luminosity (PL) relation used in the calibration of the extragalactic distance scale.

 This study aims to understand the physical origin of the IR excess found in the spectral energy distribution (SED) of 5 Cepheids : RS Pup (P=41.46d), ζ Gem (P=10.15d), η Aql (P=7.18d), V Cen (P=5.49d) and SU Cyg (P=3.85d).

 A time series of atmospheric models along the pulsation cycle are fitted to a compilation of data, including  optical and near-IR photometry, Spitzer spectra (secured at a specific phase), interferometric angular diameters, effective temperature estimates, and radial velocity measurements. Herschel images in two bands are also analyzed qualitatively. In this fitting process, based on the SPIPS algorithm, a residual is found in the SED, whatever the pulsation phase, and for wavelengths larger than about 1.2μm, which corresponds to the so-determined infrared excess of Cepheids. This IR excess is then corrected from interstellar medium absorption in order to infer or not the presence of dust shells, and is finally used in order to fit a model of a shell of ionized gas. For all Cepheids, we find a continuum IR excess increasing up to ≈-0.1 magnitudes at 30μm, which cannot be explained by a hot or cold dust model of CircumStellar Environment (CSE). 
However, a weak but significant dust emission at 9.7 μm is found for ζ Gem, η Aql and RS Pup, while clear interstellar clouds are seen in the Herschel images for V Cen and RS Pup. We show, for the first time, that the IR excess of Cepheids can be explained by free-free emission from a thin shell of ionized gas, with a thickness of  ≃15% of the star radius, a mass of 10^-9-10^-7and a temperature ranging from 3500 to 4500K.The presence of a thin shell of ionized gas around Cepheids has to be tested with interferometers operating in visible, in the mid-IR or in the radio domain. The impact of such CSEs of ionized gas on the PL relation of Cepheids needs also more investigations.


A thin shell of ionized gas explaining the IR excess of classical Cepheids
    V. Hocdé 1 
N. Nardetto 1
E. Lagadec 1
G. Niccolini1
A. Domiciano de Souza1
A. Mérand 2
P. Kervella 3
A. Gallenne1,4
M. Marengo5
B. Trahin3
W. Gieren7 
G. Pietrzyński6
S. Borgniet3
L. Breuval3
B. Javanmardi3

    Received ... ; accepted ...
====================================================================================================================================================================================================================






§ INTRODUCTION




Cepheids have been the keystone of distance scales determination in the Universe for a century, because their pulsation period correlates directly with their luminosity, through the Leavitt law <cit.> also called the period-luminosity relation (hereafter PL relation). 

The recent one percent precision on the Large Magellanic Cloud distance <cit.> has led to a new determination of the Hubble constant H_0 <cit.>. Moreover, upcoming space and ground-based telescopes such as the  James Webb Space Telescope (JWST) and the Extremely Large Telescope (ELT) will make it possible to obtain light curves of extragalactic Cepheids up to one hundred  megaparsecs.
However, this distance ladder is still largely based on Cepheids PL relation whose uncertainties on both zero point and slope are today one of the largest contributors to the error on H_0 <cit.>. One possible bias could be due to IR excesses from CSEs such as the ones discovered using near- and mid-infrared interferometry around nearby Cepheids <cit.>. Indeed, if the brightness of CSEs is found to be significantly different in the Milky Way, the SMC, the LMC and in galaxies hosting SNIa due to metallicity effects for instance, then the use of an universal PL relation could introduce a bias in the distance scale calibration.

Envelopes around Cepheids have been discovered by long-baseline interferometry in the K-Band with VLTI and CHARA <cit.>, and four Cepheids CSEs have been observed in the N band with VISIR and MIDI <cit.> and one with NACO in the near-IR <cit.>. The presence of a motionless Hα absorption component using high-resolution spectroscopy around l Car <cit.> has also been attributed to a CSE, and, recently, a CSE was detected in the visible domain with the VEGA/CHARA facility around δ Cep <cit.>. These various studies determined a  CSE radius of around 3 stellar radii and a flux contribution in the K band, ranging from 2% to 10% of the continuum, for medium- and long-period Cepheids respectively, while it is around 10% or more in the N band. However, we still do not know how these CSEs are produced, neither their nature, nor their characteristics (density and temperature profiles, chemical composition...).

This paper aims at building a phase-dependent Spectral Energy Distribution (SED) of a sample of Cepheids from visible to mid-IR wavelengths and compare it with dedicated atmospheric models in order to quantify and study their IR excess. We present the IR excess of the stars in the sample in Sect. <ref> using photometric and Spitzer   observations in various bands, and we study qualitatively far-infrared images from Herschel. In Sect. <ref> we correct the spectra from interstellar foreground absorption along the line-of-sight and seek for residuals at 9.7μm that could be due to a dusty CSE. In Sect. <ref>, we use the radiative transfer code  <cit.> to model the IR excess continuum, and show that such continuum cannot be due to a dusty CSE. In Sect. <ref>, we instead propose a model of CSE composed of ionized gas in order to reproduce the observed IR excess. Our results are discussed in Sect. <ref> and we conclude in Sect. <ref>



§ OBSERVATIONS AND DATA REDUCTION OF RS PUP, Ζ GEM, Η AQL, V CEN, AND SU CYG




 §.§ Star sample

In order to study the IR excess of Cepheids we selected a sample of Galactic Cepheids with Spitzer observations. The combination of its sensitivity and wide IR wavelength coverage from 5μ m to 38μ m, makes it the best tool to study dust features and IR excess continuum emission. We selected the five Cepheids having Spitzer low resolution spectra with the highest signal-to-noise ratio 

:  RS Pup (P=41.46d), ζ Gem (P=10.15d), η Aql (P=7.18d) , V Cen (P=5.49d) and SU Cyg (P=3.85d). In addition, we retrieved Herschel data at 70μ m and 160μ m to study larger and cold environment of dust. These space based observations have the advantage of not being biased by the Earth atmosphere, e.g. the entire flux is preserved and it is not hidden by ozone band which is essential to study dust spectral features. 



 §.§ Building the infrared excess using the SPIPS algorithm

SpectroPhoto-Interferometric modeling of Pulsating Stars (SPIPS) is a model-based parallax-of-pulsation code which includes photometric, interferometric, effective temperature and radial velocity measurements in a robust model fit <cit.>. To compute synthetic photometry to match those from the dataset SPIPS uses a grid of ATLAS9 atmospheric models[<http://wwwuser.oats.inaf.it/castelli/grids.html>] <cit.> with solar metallicity and a standard turbulent velocity of 2km/s.

SPIPS models are available for η Aql <cit.>, ζ Gem <cit.> and RS Pup <cit.>. We updated them with new datasets including Gaia photometry <cit.> and we also applied the SPIPS algorithm to SU Cyg and V Cen. Fig. 1 to 5 in Appendix <ref> show the results of the fitting process. For each star in the sample, SPIPS provides a fit of the photometry along the pulsation cycle, which is in agreement with the observations. While the fit is satisfactory from the visible domain to the near infrared, a disagreement is obtained for wavelengths larger than about 1.2μm. For λ > 1.2μm, the observed photometries (m_obs) are indeed significantly brighter than the synthetic ones (m_kurucz), which corresponds to an infrared excess. SPIPS does not physically model this IR excess, but instead takes it into account using an ad-hoc analytic power law IR_ex, which is defined as:







   



    IR_ex=Δmag = m_obs - m_kurucz= 0,    λ < 1.2μm
    α (λ - 1.2)^β,    λ⩾ 1.2 μm
 
with two parameters, α and β. 

These ad-hoc laws take the visible domain as a reference, and thus assume that there is no excess nor deficit in this wavelength domain due to the CSE (Δm=0 for λ < 1.2 μ m). The choice of the ad-hoc analytic power law in Eq. (<ref>) was adopted originally in SPIPS modeling because it provides a reasonable description of the IR excess without changing the results of the fit.
 
In  Fig. <ref>, we show the IR excess analytic law obtained for the star sample together with the measurements. These measurements correspond to the cycle-averaged magnitude difference Δmag = m_obs - m_kurucz (in a specific band), while the uncertainties are the corresponding standard deviations. Since these deviations over the cycle are quite small (≈ 0.05mag) from Fig. <ref>, we conclude that the IR excess of Cepheids is not time dependent or, at most, slightly time dependent. 


In the fitting process described above, the photometric band corresponding to the carbon-monoxyde (CO) band-head at 4.6μ m has to be considered with a special care. Indeed,   CO can form at low temperature <cit.> and this is indeed probably seen in IRAC I2 Spitzer bands for ζ Gem and V Cen (see Figs. <ref> and <ref>, respectively). In these bands, the agreement between the ATLAS9 model and the observations is in disagreement at certain phases. Even if ATLAS9 includes the modeling of CO molecules, these models are static and do not reproduce the dynamical structure of the atmosphere of Cepheids satisfactorily. However, the cycle-averaged infrared excess we obtain at this specific wavelength of 4.6μ m is consistent with the general trend of the IR excess law (see Fig. <ref>), and is considered in the following of the analysis. 

Besides, the photometric bands longward 5 microns are problematic for two reasons. First their cycle coverage is poor i.e. one data point in each band with undetermined phase, this is the case for η Aql, ζ Gem and RS Pup at λ>5μm (see for instance the A MSX band of RS Pup in Fig. <ref>). Hence the reliability of the fitted ad-hoc IR excess model depending on pulsation phase is questionable. Secondly these bands are more sensitive to interstellar dust emission around the star than Spitzer because of their lower spatial resolution. Thus we have decided to discard these photometric bands from the analysis. The measurements which are not considered are indicated by orange bars in Fig. <ref>.






 




















  





   


 §.§ Spitzer data

The spectroscopic observations were made with the InfraRed Spectrograph IRS <cit.> onboard the Spitzer telescope <cit.> and the full spectra were retrieved from the CASSIS atlas <cit.>. CASSIS identified the five stars as point-like sources therefore we retrieved the best flux calibrated spectrum obtained from optimal extraction method <cit.>.

The Spitzer dataset of this paper is presented in Table <ref>. Short-Low (SL) and Long-Low (LL) are IRS modules placed in the focal plane instrument providing low spectroscopic resolution (R=60 - 128) from 5.2 μ m to 38 μ m. 
In this study we use low-resolution IRS spectra obtained with the SL and LL modules, for which aperture sizes together with covered spectral ranges are described in Table <ref>. 




We derive the IR excess of each star in the sample at the specific phase of Spitzer using:

    Δmag = m_Spitzer - m_kurucz[ϕ_Spitzer]

where  m_Spitzer is the magnitude of the Spitzer observation and m_kurucz[ϕ_Spitzer] is the magnitude of the ATLAS9 atmospheric model interpolated at the phase of Spitzer observations (ϕ_Spitzer). The T_eff and log g values of the star at the phase of Spitzer are provided by the SPIPS algorithm, while the interpolation is then done in a ATLAS9 grid of models with steps of  250K in effective temperature and 0.5 in log g, respectively. The angular diameter derived by SPIPS is then used to calculate m_kurucz[ϕ_Spitzer]. The corresponding stellar parameters are summarized in Table <ref> and the Δmag values we finally obtain are shown by green dots and red uncertainties in Fig. <ref>.   The same is done for all the photometric bands of observations, i.e. we do not use the average IR excess values shown in Fig. 1, but instead we recalculate the actual values corresponding to the phase of Spitzer observations (red dot with blue uncertainties in Fig. <ref>). Note that in Fig. <ref>, the cycle averaged ad-hoc analytic IR excess law is shown only for comparison.


First, we correct several discontinuities in the Spitzer spectra. Even if absolute uncertainties <cit.> are not excluded, we find a discontinuity at 14μ m between SL and LL spectra, indicated by a dashed line in Fig. <ref>. According to CASSIS atlas a residual offset appears in the case of slightly extended sources <cit.>, because more flux enters the larger LL aperture. Indeed the two Spitzer detectors have indeed different Fields Of View (FOV, see Table <ref>) with different orientations on the sky (see Fig. <ref>), which explains the different level of flux measured. As an indication, the FOV of both detectors are overlaid on Herschel images (see Fig. <ref> in Sect. 2.4). 

For RS Pup, we obtain a jump of 7% in flux between the SL and LL detectors, while for other stars we find a difference of around 2% (1.3% in the case of ζ Gem). Extended dust emission has been discovered around RS Pup <cit.>, so that this jump is not surprising. 



We correct these discontinuities at 14μ m in the spectra by scaling up SL flux (with the latter flux ratio) so that it corresponds to LL flux at 14 μ m. This correction assumes that the LL flux calibrated by CASSIS pipeline is reliable and the emission of the environment is homogeneous between SL and LL apertures.



Also, orders mismatches can appear because of telescope pointing accuracy and causing discontinuity in the spectra between 2 adjacent orders. We corrected this feature in SU Cyg and ζ Gem spectra at 7.5 μ m between SL orders 1 and 2 by scaling up the lower part. Note also that the Spitzer data at larger wavelengths than 30μm were not considered due to their extremely large uncertainties. 

Second, for the five stars in the sample, we find an agreement between the level of IR excess at 5μm in all photometric bands  (λ < 5μm) and in the Spitzer spectroscopic observation (λ > 5μm). This is verified at the 0.05 magnitude level, and before the correction of different Spitzer discontinuities (see Fig. <ref> for a detailled analysis of corrected IR excess). 
The agreement between IR excesses from SPIPS models and Spitzer data shows the consistency of this approach. It  reveals that Cepheids exhibit a continuum IR excess from about 2μm up to 30μm in the star sample. This result is different from <cit.> who found no IR excess for short-period classical Cepheids (see Discussion <ref>). We find however that the ad-hoc analytic laws of SPIPS overestimate the IR excess in the Spitzer wavelength domain (see Fig. <ref>).

Third, RS Pup, SU Cyg and V Cen present unambiguous silicates absorption bands at 9.7μ m and 20μ m. The latter absorption is known to peak usually at 18μ m however it can be shifted to longer wavelength depending on various conditions such as temperature or minerals proportions into silicates <cit.>. These silicate absorptions are likely due to the presence of interstellar clouds on the line of sight. 
 








 §.§ Herschel images 

We retrieved Herschel images from the PACS instrument at 70 and 160μ m to study the presence of extended environment around cepheids. No data have been found for SU Cyg. Herschel observation products are ordered according to the level of the data processing ranging from raw data (level 0) to highly processed scientific data (level 3). All the images in the sample are highly processed data with a level of 2.5. The FOV of Herschel is of 8' and its spatial resolution is of  7” at 100μ m. As an indication, we performed aperture photometry measurements of the flux within 1' in order to include most of the emission due to the environment of the cepheids. Then we derived the IR excess at 70 and 160μ m by comparing with atmospheric model in the same way than Eq. <ref>. 

Observations are presented in Fig. <ref> and the derived photometry in Table <ref>. 
From these data we observe important IR excesses at 70μ m and 160μ m, which are associated to cold environments with temperatures of about 40K and 20K respectively according to Wien's displacement law (λ T = 2900μm.K). Cold dusty material around RS Pup  was already observed, while the large cloud observed around V Cen (see Fig. <ref>) at 160μ m could be related to a star forming region, since V Cen is likely a member of the open cluster NGC 5662 <cit.>. Herschel is sensitive not only to emission from the star itself and a potential CSE, but also from cold extended emission from the interstellar medium. We thus do not consider the Herschel photometry in the following of the paper.

Interestingly, we find a qualitative correlation between silicate absorption features in Spitzer data and the presence of extended emission in the Herschel images. For instance, both RS Pup and V Cen show at the same time strong silicate absorptions and extended emission in the Herschel data (see <ref>), while η Aql and ζ Gem present no obvious silicate absorption and weaker interstellar environment (<ref>). 


    




§ CORRECTING THE INTERSTELLAR SILICATE ABSORPTION IN SPITZER DATA

Absorption of silicates at 9.7μ m is easily identified in the Spitzer spectra of RS Pup, V Cen and SU Cyg whereas it is not observed in η Aql and ζ Gem data (see Fig. <ref>). Such absorption has in principle two components, an emission due to a CSE close to the star (if present), which is then absorbed by interstellar environment. Thus, the silicate absorption observed by Spitzer A^Spitzer_9.7, if it is corrected from the expected silicate absorption from the ISM (A^ISM_9.7), can indicate, in case of residual emission, whether there is a silicate emission from the CSE E^CSE_9.7 or not. We thus have the following relation:




    A^Spitzer_9.7 - A^ISM_9.7  = E^CSE_9.7

 All values are magnitudes  at 9.7μ m. We will now determine A^Spitzer_9.7 (Sect. <ref>) and A^ISM_9.7 (Sect. <ref>) in order to estimate the CSE emission at 9.7μ m (Sect. <ref>) and then, we will extend our correction to the whole wavelength range of  Spitzer (Sect. <ref>).





 §.§ Quantifying the observed silicate absorption from Spitzer

In order to estimate the apparent silicate absorption at 9.7μ m (A^Spitzer_9.7), we have to define the IR excess continuum. For that, we fit the IR excess continuum on each side of the absorption feature over the wavelength range [6,7]∪ [12,13] (dashed line in Fig. <ref>), and we interpolate the excess continuum value at 9.7μ m (blue dot in Fig. <ref>). We then subtract to this value the IR excess corresponding to the core of the silicate absorption at 9.7μ m (red dot). The uncertainty on A^Spitzer_9.7 is derived by adding in quadrature both the error on the continuum and the error on the Spitzer observation at 9.7μ m. For ζ Gem we do not find any evidence for absorption.





 §.§ Estimation of the interstellar medium (ISM) silicate absorption.

 The ISM extinction theory can be used to estimate the silicate absorption directly from the color excess E(B-V). In order to select the most reliable E(B-V) values we compared for each star (1) the E(B-V) fitted by the SPIPS algorithm, (2) the E(B-V) from the David Dunlap Observatory Database of Galactic Classical Cepheids[<http://www.astro.utoronto.ca/DDO/research/cepheids/>]<cit.> (hereafter DDOD). This database computes the mean extinction and standard error from various values obtained in the literature (3) the E(B-V) obtained with the Stilism 3D map <cit.> using the bayesian inversion method to a wide color excess dataset and parallaxes (4) the 2D map from <cit.> (hereafter SFD98) using 100μ m data from COBE/DIRBE and IRAS. In the latter 2D map we have corrected the extinction by taking into account the star's location using a three-dimensional model of the Milky Way from <cit.>. We present these E(B-V) values in Table <ref>. The best agreement is obtained when comparing SPIPS and DDOD values, two totally independent methods, which reinforces the reliability of the SPIPS fitting. 
 

For the Stilism and SFD98 2D maps, the E(B-V) values are based on models of dust distribution within the Milky Way, and thus cannot consider local over- or under-densities in the vicinity of the Cepheids, which probably explains the large discrepancies found in particular for RS Pup, SU Cyg, and also V Cen. For these stars the uncertainty on E(B-V) is also particularly large in the case of Stilism. The E(B-V) values provided by SPIPS are the most precise by an order of magnitude but uncertainties are underestimated <cit.>. Therefore we decided to adopt an independent and more conservative approach using the DDOD values.





 
 
 

 
 As a second step, we derive the visible absorption A_v assuming an extinction law A_v=R_vE(B-V) with a ratio of total-to-selective extinction of R_v=3.1 which corresponds to a diffuse ISM along the line of sight <cit.>. Then we used the relation A_v/τ_9.7=18.5 <cit.> which is suited to diffuse ISM in the solar vicinity, in order to derive τ_9.7. This was done for each star except RS Pup. Indeed, RS Pup is known for being embedded in a large environment which is most probably the remnant of the molecular clouds in which the star formed <cit.>. Hence we treated RS Pup as a special case since the extinction towards dense clouds is different from diffuse ISM <cit.>. We assumed most of the extinction is due to a dense cloud and we used R_v=5 and A_v/τ_9.7=11.46 following the empirical work in star forming regions  by <cit.>. 
Finally, the extinction A_λ is given by the intensity absorption along the line of sight considering an optical depth τ_λ. 
Combining both the Beer-Lambert law and absorption definitions for any wavelength λ we have A_λ=1.086τ_λ, thus A_9.7=1.086τ_9.7 at the specific wavelength 9.7μm.

This approach can be finally summarize using using two equations :
[ht!]

    A^ISM_9.7 = 1.0863.1/18.5=0.182 E(B-V),

Assuming silicate reddening law for diffuse ISM.



[ht!]

    A^ISM_9.7 =1.0865/11.46= 0.474 E(B-V),

Assuming silicate reddening law for dense clouds.

for stars in diffuse ISM (4)  and for RS Pup (5).

The final values of A^ISM_9.7 we consider are listed in Table <ref>, with their corresponding uncertainties.







 §.§ Residual silicate CSE emission at 9.7μ m 

Using values of A^Spitzer_9.7 (Sect. 3.1) and A^ISM_9.7 (Sect. 3.2), we now calculate E^CSE_9.7 using Eq. (<ref>). 
Results are illustrated in Fig. <ref>.

Using this method, silicate emission from CSE exists only when E^CSE_9.7 < 0. We find no residual CSE emission for SU Cyg and V Cen, and significant but weak emissions for ζ Gem (-0.008 mag± 0.004), η Aql (-0.019± 0.004 mag) and RS Pup (-0.110 ± 0.057 mag). In the case of V Cen, the

emission is likely underestimated when using Eq. (<ref>). If we assume instead the presence of dense cloud on the line of sight using Eq. (<ref>), we obtain a silicate emission of -0.091±0.013 mag (see Fig. <ref>). 

From this figure (and Table <ref>), it is difficult to conclude whether there is a dusty CSE of silicate around the Cepheids in the sample or not, in particular because our method is sensitive to the Equation used (<ref> or <ref>) when correcting the silicate ISM absorption (case of V Cen for instance), and to various sources of uncertainties, in particular in the E(B-V) estimate. We also assume that the IR excess is not varying with time, as suggested by the SPIPS analysis of Sect. <ref>. But still, with this method, we find residual dust emission at 9.7μm for the long-period Cepheids in the sample (η Aql, ζ Gem and RS Pup).










 §.§ Extending the correction of the silicate absorption from 9.7μm to the whole wavelength range of Spitzer.

We correct the entire Spitzer observations by subtracting a synthetic interstellar medium composed of silicates. Since we assumed an averaged ISM temperature of 20K, the dust emission is negligible in the Spitzer wavelength range according to Wien's law. Thus, we simply derive the absorption A^ISM_λ analytically using Mie theory. In the calculation we take into account only the effective absorption cross-section C^abs_λ and we neglect the scattering effects since the radius of grain a is small compared to the mid-IR wavelength. Hence we adopted the following expression for λ between 5 and 30 μ m

    A^ISM_λ∝κ_λ=∫ C^abs_λ(a)π a^2 n(a)da

We first derived C^abs_λ using complex refractive index for silicates from <cit.> (hereafter DL84) assuming an uniform distribution of ellipsoidal shapes given by <cit.>. Then we derived the absorption coefficient κ_λ by taking into account a standard grain size distribution n(a)∝ a^-3.5 <cit.>. Finally we normalize A^ISM_λ using its specific value A^ISM_9.7 at 9.7μ m we already derived in Sect. <ref>.




































For all the stars in the sample, we plot the IR excess (found to be constant with the pulsation phase, see Sect. <ref>) corrected from the ISM silicate absorption in Fig. <ref>. In the following sections we first show that a dust envelope cannot reproduce the IR continuum excess, then we explore the possibility of a thin envelope of partially ionized gas in order to model the IR excess continuum.




§ INCOMPATIBILITY OF DUST CSE MODEL TO EXPLAIN IR EXCESS CONTINUUM.

Cepheids are oxygen-rich stars, and most of the carbon in their envelopes is thus locked in CO molecules. In the condensation sequence described by <cit.> corundum (Al_2O_3) is expected to form first because of its high condensation temperature of about 1400K in typical pressures encountered in circumstellar shells. For lower temperatures, Al_2O_3 is depleted and silicates such as gehlenite (Ca_2Al_2SiO_7) and forsterite (Mg_2SiO_4) can form. All these components present emission components mostly observed in the N band (8-13 μm, see Fig. <ref>), with a width of at most a few microns. The IR excess we observe in Spitzer and SPIPS data does not present any clear spectral feature and is broader than 20 microns.  Fig. <ref> shows the best fit we obtain with silicate dust. Silicate dust features are clearly unlikely to explain the observed IR excess continuum from near- to mid-IR. This conclusion is in agreement with the work made by <cit.> who has found that CSE made of silicate dust cannot explains IR excess of 132 classical and type-II Cepheids. However, since the opacity of iron exhibit no particular feature but rather a continuum (see Fig. <ref>) we investigate whether a warm dust envelope of iron could explain the IR excess continuum observed with SPIPS and Spitzer dataset.

The CSEs were modeled using  <cit.>, which solves the radiative transfer equations in a dusty environment. The method is based on a self-consistent equation for the spectral energy density, including dust scattering, absorption, and emission. We present 2 CSEs models with different inner shell temperatures. Iron is a component in oxygen-rich star mineralogy with a condensation temperature of ≈ 1000K for circumstellar pressures. Hence, the cold model takes into account the typical condensation temperature of iron (≈1000K) as the inner shell temperature which is thought to be realistic. A reduced χ^2 fitting is applied to adjust the optical opacity τ_0.55. 
 The hot model lets both the inner shell temperature and the optical opacity as free parameters during χ^2 fitting. In both hot and cold model we used a standard MRN size distribution <cit.>, and we computed the density distribution in the case of an envelope expansion which is driven by radiation pressure on the dust, the wind structure is derived taking into account the dust drift and the star's gravitational attraction. The outer shell radius is set to 500 times the inner shell radius. We test the consistency of these models on η Aql using interpolated ATLAS9 atmospheric models (see Table <ref> in Sect. <ref>) as a central source's radiation in the  computation.
The results of the computation are summarized in Table <ref> and the computed IR excess is presented Fig. <ref>. 
Both models fails to reproduce the IR excess continuum modeled with SPIPS and observed with Spitzer. The cold model is well below the observed IR excess and the temperature required for iron in the hot model is much higher than iron condensation temperature for circumstellar shell pressures. In that case  solid iron would not form or would be sublimated. Last, we can also argue that a CSE made of iron exclusively for each Cepheid is very unlikely. In conclusion we do not expect a warm dust envelope to cause the IR excess continuum of Cepheids  in the 1 - 30 μ m range. 













§ THE IR EXCESS FROM A THIN SHELL OF IONIZED GAS

We investigate hereafter if an ionized, spherical gas shell can explain at the same time the near- (∼1 μm < λ < 5 μm) and mid (5 μm < λ < 30 μm) IR excess of Cepheids. As a first step, this study is made only at the specific pulsation phase of Spitzer data. We discuss the possible time-variability of such shell in Sect. <ref>



We consider the emission of a thin gas shell around the star with constant density and temperature for the sake of simplicity. The shape of the mid-IR excess, saturating to a constant flux ratio at large wavelengths (see Fig. <ref>), suggests an opacity source increasing with wavelength. We used the free-free and bound-free opacities for a pure H shell presenting such behaviour. The combined absorption coefficient (in m^-1; SI(MKS) unit system) for these two opacities sources is given by <cit.>







    κ_λ=3.692×10^-2[1-^-hc/λ k ] ^-1/2×(λ/c)^3
    
    	(γρ/)^2[(λ,)+(λ,)],


where h, c, and k have they usual meanings, γ is the degree of ionization (between 0 and 1),  the temperature of the shell,  the hydrogen mass and  and  are the free-free and bound-free Gaunt factors respectively. These factors were computed mainly from approximation formulas given by <cit.>, <cit.>, and references therein. As an example we present a typical shape of the absorption coefficient κ_λ computed for V Cen in Fig. <ref> using the shell parameters from Table <ref>.


We computed the SED of the star plus the gas shell as described in Appendix <ref>, taking into account the latter absorption coefficient κ_λ. In order to match the SPIPS photometries plus corrected Spitzer spectra presented in Sect <ref> we performed a χ^2 fitting using the Levenberg-Marquardt method. We used filters bandpass to convert the flux of the physical model into the corresponding SPIPS photometries for the fitting procedure. In addition, since the SPIPS fitting assumes that there is no excess in the visible domain i.e Δm=0 for λ < 1.2 μ m (see Sect. <ref>) it is necessary to relax this assumption to allow the data to present deficit or excess in the visible depending on the physical behaviour of the ionized shell. Soften this assumption is equivalent to suggest an improvement in the SPIPS physical treatment of the circumstellar environment. As a first step, we consider a simple thin, spherical, and partially ionized gas shell in order to reproduce the IR excess of the Cepheids in the sample. Hence, we fitted 4 parameters: three parameters from the gas shell i.e the ionized shell mass γ, its temperature  and radius  plus one parameter corresponding to the IR excess offset corresponding to Δm≠0 for λ < 1.2 μ m.
Since Spitzer data have a higher statistical weight than SPIPS data we first find the different parameters by hand to use it as first guesses. The results are presented in Table <ref> and in Fig. <ref>.

In Fig. <ref> we observe discontinuities at short wavelengths, from visible to near-IR, which are due to bound-free opacities (see Fig. <ref>). Bound-free opacity decreases for longer wavelength whereas free-free opacity increases. Beyond 5 μm the bound-free contribution is negligible compared to the free-free contribution which explains the observed smooth continuum. For the Cepheids in the sample, we obtained shell temperatures between 3500 and 4500K and an envelope thickness of ≃ 15% of the radius of the star, while γ is ranging from 10^-9-10^-7. The IR excess correction offsets found have positive values which means Δm⩾0 for λ < 1.2 μ m for all stars except for SU Cyg which presents a slight negative value. Indeed the several computed models present absorption in the visible domain (see Fig. <ref>), thus correction offsets have positive values to allow a deficit in the visible due to ionized shell absorption.



In the prospect of JWST observations, we provide the synthetic IR excess anticipated from the shell of ionized gas in several bands (Table <ref>). We used filter transmissions of wide filters: F070W and F200W from Near-Infrared Camera (NIRCam)[<https://jwst-docs.stsci.edu/display/JTI/NIRCam+Filters>] and F560W and F1000W from Mid-Infrared Instrument (MIRI)[<https://jwst-docs.stsci.edu/display/JTI/MIRI+Filters+and+Dispersers>]. These filters, centered on 0.70μ m, 2μ m, 5.6μ m and 10μ m respectively, are suitable to measure the IR excess continuum from the ionized shell model.



§ DISCUSSION



 §.§ The IR excess and dust environment 

We estimated the infrared excess of Cepheids using SPIPS algorithm (λ <5μm) and Spitzer data (5 < λ < 30μm). The observed IR excess presented in Fig <ref> is thus calculated at the specific phases of Spitzer and makes the assumption that there is no excess and/or deficit in the visible domain due to the CSE (Δm=0 for λ < 1.2 μ m). This leads to the following conclusions. 

First, these IR excess emissions show a continuum which is consistent across the wavelength range for all stars. It ranges in the data from about 2μm to 30μm and corresponds to differences of magnitudes of up to -0.1 and even -0.2 magnitudes for RS Pup and η Aql, between visible and far-infrared. Importantly, this IR excess rises in the near-infrared (around 2μm) from about 0 magnitude of differences (assuming Δm=0 for λ < 1.2 μ m) to -0.1 magnitude around 5μm, for each star, which brings strong constrains on the models and in particular invalidate a pure CSE of dust to explain the IR excess continuum. 

Second, in order to unveil CSE emission in the N-band we performed an independent correction of the ISM silicate absorption (Sect. <ref>) which seems to fill almost perfectly the silicate absorption seen in the Spitzer data. We have determined the excess at 9.7μ m and we have found that there is no emission within the uncertainty for SU Cyg and V Cen, while weak emissions are found for ζ Gem, η Aql and RS Pup. 

The slight residuals that we found from our spectroscopic analysis at 9.7μm could be attributed to faint dusty CSE around Cepheids with a flux contribution (compared to the stellar flux) of few percents, which is in average ten times less than what was predicted by other studies based on interferometry <cit.>, and consistent with <cit.> who found no silicate emission for a large sample of Cepheids.


Third, using photometric bands and basically the same approach but without considering the pulsation of the Cepheids, <cit.> found that 21 over 132 classical and type 2 Cepheids in his sample has a clear or weak IR infrared excess. Moreover, these 21 Cepheids have periods larger than 11 days. Conversely, we have 4 short periods Cepheids in the sample (except RS Pup) and all of them show a clear IR excess. However, there is no star in common between the two studies to go deeper in the analysis. 





Four, we find cold, large and inhomogeneous circumstellar environment around the Cepheids RS Pup, V Cen, which is seen in the Herschel images. As Cepheids are relatively young stars, they are still likely close to the cloud where they formed and such environment can contribute to the overall IR excess, either in absorption or emission. This could also affect the mid-infrared photometry of Cepheids, and thus potentially estimations of the PL relation when using instruments such as the JWST and forthcoming mid-infrared instruments on the future 40 meter-class telescopes.




 §.§ The IR excess explained by a thin shell of ionized gas

We show that a thin shell of ionized gas with a temperature ranging from 3500 to 4500K depending on the Cepheid considered, and with a width of typically ≃ 15% of the radius of the star can reproduce the IR excess. Up to now, the only attempt to detect ionized material were carried out by the Very Large Array (VLA) at 5GHz on η Aql and 4 others classical Cepheids <cit.>. Since no 3σ detection has been reported, only upper limits on flux density were derived. From the ionized shell model presented in this paper, we derived a flux density between 0.01 and 0.1 μJy at 5GHz (≈20% above the star continuum) which is below the upper limit on flux density of ≈100μJy estimated by <cit.>.


On the other hand interferometric observations have resolved CSEs around Cepheids. The first detection was reported around l Car <cit.> followed by δ Cep and Polaris <cit.>. These CSEs were modeled with a ring at a distance of 2 to 3R_⋆, i.e. close to the star, in a region where the temperature is high enough (>2000K) to prevent dust condensations <cit.>. Thus, these observations are more likely explained by a shell of partially ionized gas. This occurence in detection by interferometry plays in favor of a widespread phenomenon amongs classical Cepheids.



Also, extensive studies of Hα profiles in the atmosphere of short-, mid- and long-period Cepheids have shown that strong increases of turbulence occurs when the atmosphere is compressed during its infalling movement, or because of shock waves dynamics <cit.>. In the case of long period Cepheids, several shock waves can be observed and P Cygni profiles show that there is an expanding shell of Hα emission which is detached from the photosphere <cit.>. Besides, analytical work of <cit.> have shown that mass loss is enhanced by pulsations and shocks in the atmosphere. The present paper suggests this mass loss could be in the form of partially ionized gas.



The model of the shell of ionized gas could also be linked to the chromospheric activity of Cepheids. <cit.> report HeIλ10830 observation on seven Cepheids providing the evidence of a high temperature plasma and steady material outflow in the highest part of the atmosphere. In addition, high resolution profiles of Mg II h and k lines (2900A) of five Cepheids with International Ultraviolet Explorer instrument <cit.> revealed extended dynamical chromospheres up to tenth of stellar radii composed of rising and falling material. Such extensions could be compatible with the shell thickness presented in this paper. Outflowing materials, with velocity of 50 to 100 km/s, can in turn eject material to several stellar radii if the velocity lasts through the pulsation cycle. These strong mass transfers could explain the enigmatic X-rays detections with Chandra and XMM Newton <cit.>.



In parallel, it is interesting to compare Cepheids to very long period Mira stars for which a radiosphere near 2R_⋆ due to free-free emission has been reported <cit.>.
 


 §.§ Limit of the model and perspectives 

We have found relatively low temperatures for the hydrogen ionization (between 3500 and 4500K see Table <ref>) and only a small fraction of the gas should be ionized according to Saha equation. In particular under 4000K the number of free electrons should be mostly provided by metals such as iron or aluminum which have low ionization potential. This effect could in turn produce a fainter ionized shell for low-metallicity Cepheids in the Magellanic clouds. Thus, it would be interesting to quantify theoretically the impact of metallicity on the shells of ionized gas, and the PL relation. Moreover, our model does not take into account temperature nor density gradients in the star's atmosphere, and in particular compression and/or shock waves which could also heat up the shell and ionize the gas. 

There are indications from SPIPS analysis that the IR excess of Cepheids is not or only slightly time-dependent. Evidences of slight cycled variations exist due to the opacity change at 4.5μ m caused by periodic formation and destruction of CO molecules in the atmosphere <cit.>. As we have Spitzer data only at a specific phase of pulsation for each Cepheid, we cannot firmly conclude for the time-dependency of the mid-IR excess. Nevertheless, if we assume that the IR excess is constant, then, as our thin shell of ionized gas is close to the star (≃ 15% of the radius of the star), it is supposed to be dynamic and its parameters should vary with time. Doing the test for η Aql we find a rather stable relative size of the shell (≃ 15 ±2 %) and a temperature variation from about 3500 to 4500K, which is similar to the values we obtained from one star to the other in the sample.







Last, our simple model suggests that the thin gas shell absorbs the light coming from the star in the visible domain (from 0.01 up to 0.13 mag), which invalidates the initial assumptions of Δm=0 for λ < 1.2 μ m in the SPIPS algorithm. We obtain satisfactory results only if we apply a correction offset Δm>0 for all wavelengths (except for SU Cyg with a slight negative offset). To explain this shift, we suggest that the distance found by the SPIPS algorithm might be too large by few percents (factor 10^Δm/2.5), all other parameters being unchanged. In other words, if obscured by a shell of ionized gas, Cepheids could be slightly closer than expected by SPIPS.
On the interferometric side, the angular diameter of the star would be also lower by few percents, but this can be compensated by the effective size of shell, and would have little impact in the SPIPS fitting. Thus, a spatial and chromatic analysis of the shell including interferometric constrains in all available bands with in particular VEGA/CHARA (visible), PIONIER/VLTI (infrared) and MATISSE/VLTI (L, M, N bands) is still necessary to better understand the environment of Cepheids, and eventually, check the impact on the PL relation.



§ CONCLUSION



  * For the five Cepheids, we report a continuum IR excess increasing up to ≈-0.1/-0.2 magnitudes at 30μm, which cannot be explained by a hot or cold dust model of CSE.

  * Within the limits of our assumptions, we do not firmly conclude about the presence of CSE emission in the N-band but it is likely weak (> -0.1mag) according to our results. 

  * We show for the first time that IR excess of Cepheids can be explained by a free-free emission of a thin shell of ionized gas with a thickness of ≃8-17% star radius, an ionized mass of 10^-9-10^-7and a temperature of 3500-4500K. In this simple model, density and temperature have a constant radial distribution.
 
  * The presence of a thin shell of partially ionized gas around Cepheids has to be tested with interferometers operating in visible, in the mid-IR or in the radio domain. The impact of such CSEs of ionized gas on the PL relation of Cepheids needs also more investigations.
























 






The research leading to these results has received funding
from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme under grant agreement No
695099 (project CepBin). The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-15-CE31-0012- 01 (project UnlockCepheids). W.G. and G.P. gratefully acknowledge financial support for this work from the BASAL Centro de Astrofisica y Tecnologias Afines (CATA) AFB-170002.We acknowledge financial support from “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU, France. This project was partially supported by the Polish Ministry of Science grant Ideas Plus. This research made use of the SIMBAD and VIZIER[Available at < http://cdsweb.u- strasbg.fr/>] databases at CDS, Strasbourg (France) and the electronic bibliography maintained by the NASA/ADS system. This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. PACS  has  been  developed  by  a  consortium  of  institutes  led  by  MPE  (Ger-many) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA,LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA(Italy); IAC (Spain). This development has been supported by the funding agen-cies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR(Germany),  ASI/INAF  (Italy),  and  CICYT/MCYT  (Spain). 
This research also made use of SAOImage DS9,  an  astronomical  imaging  and  data  visualization  application <cit.> and  Astropy, a  community-developed corePython package for Astronomy <cit.>.




aa  
 

 


§ SPIPS DATA SET AND FITTED PULSATIONAL MODEL OF THE STAR SAMPLE.

The plots are organized as follows: the pulsational velocity,  the effective temperature and the angular diameter curves according to the pulsation phase are shown on the left panels while  the right panels display photometric data in various bands. Above each figure is indicated the projection factor set to p=1.270, the fitted distance d using parallax-of-pulsation method, the fitted color excess E(B-V), and the ad-hoc IR excess law. In the photometric panels, the gray dashed line corresponds to the magnitude of the SPIPS model without CSE. It actually corresponds to the magnitude of a Kurucz atmosphere model, m_kurucz, obtained with the ATLAS9 simulation code from <cit.> with solar metallicity and a standard turbulent velocity of 2km/s. The gray line corresponds to the best SPIPS model, which is composed of the latter model without CSE plus an IR excess model. Note that for WISE, MSX and IRAS filters observations above 5 microns, only one data point is obtained without information on the phase. Hence it is represented by a horizontal gray strip for which the vertical width is the uncertainty of the measurement. In the angular diameter panels the gray curve corresponds to limb-darkened (LD) angular diameters.
For metallic stars, when effective temperature is low enough, CO molecules can form in the photosphere and absorb light in the CO band-head at 4.6 μ m <cit.>. This effect is observed in the Spitzer I2 IRAC dataset of ζ Gem and V Cen (see Fig. <ref> and <ref>). In this case, these data were ignored during the fitting of SPIPS. When no effective temperatures and no angular diameters are included in the SPIPS model, there is a degeneracy between the mean temperature and E(B-V). We estimate that SPIPS can make an error of +/- 200K on the effective temperature and +/- 0.05 on E(B-V). Only V Cen has no data for both effective temperature and angular diameter, nevertheless the fitting of SPIPS is thought to be reliable (for example see E(B-V) value compared with literature in Table <ref>).



























§ THE IR EXCESS OF A THIN GAS SHELL AT CONSTANT TEMPERATURE AND DENSITY 



The shell emission is obtained integrating the radiative transfer equation along rays defined by their impact parameter p (see Fig. <ref>) following the method described in <cit.>.









We assume a constant shell temperature  and density . The gas opacities (see Eq. (<ref>)) can be written under the form κ(λ,)=^2 χ_λ().

According to Fig. <ref> we have to take into account two cases corresponding to the impact parameter p respectively larger and smaller than the stellar radius .

For p⩾, taking into account the symmetry of the problem the optical depth along the ray is given by

    τ_λ(p)=2∫_0^√(R^2_out-p^2)κ(λ,)ds  ,

which, for a constant  and density  gives

    τ_λ(p)=2^2χ_λ() √(R^2_out-p^2) .


Similarly, for p⩽ we have

    τ_λ(p)=^2χ_λ() [√(R^2_out-p^2)-√(^2-p^2)]  .



For p⩾, the specific intensity is given by 










    I_λ(p) = B_λ()(1-e^-τ_λ(p))  ,

since  is assumed to be constant along the ray.

For p⩽, both the shell and the stellar photosphere contribute to the specific intensity 

    I_λ(p) = B_λ()(1-e^-τ_λ(p))+ I^∗_λ e^-τ_λ(p) ,

where I^∗_λ is the stellar specific intensity.






The observed total emerging flux at a distance d is then computed numerically by quadrature with the following integral

    F_λ=2π/d^2∫^R_out_0pI_λ(p)dp  .





























From the result of the integral in Eq. (<ref>) we can deduce the magnitude excess defined by

    Δmag=-2.5log(F_λ/F^∗_λ) ,

with F_λ^∗=π(/d)^2 I^∗_λ.

Note that in the particular case where we have I^∗_λ = B_λ() with =, Eq. (<ref>) can be integrated analytically to give

    F_λ=π(/d)^2 B_λ() ×
    [1 + [(R_out/)^2-1] [1+2/τ^∗ 2_λ[(1+τ^∗_λ)e^-τ^∗_λ-1]]]  ,

with τ^∗_λ defined as 2 ^2 χ_λ()   √((R_out/)^2-1).




The corresponding magnitude excess is given by

    Δmag=2.5log[1+[(R_out/)^2-1]
    [1+2/τ^∗ 2_λ[(1+τ^∗_λ)e^-τ^∗_λ-1]]]

 
