A thin shell of ionized gas as the explanation for infrared excess among classical Cepheids

Context. The infrared (IR) excess of classical Cepheids is seldom studied and poorly understood despite observational evidence and the potential for its contribution to induce systematics on the period-luminosity (PL) relation used in the calibration of the extragalactic distance scale. Aims. 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). Methods. A time series of atmospheric models along the pulsation cycle were ﬁtted to a compilation of data, including optical and near-IR photometry, Spitzer spectra (secured at a speciﬁc phase), interferometric angular diameters, effective temperature estimates, and radial velocity measurements. Herschel images in two bands were also analyzed qualitatively. In this ﬁtting process, based on the SPIPS algorithm, a residual was 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 was then corrected from interstellar medium absorption in order to infer the presence (or absence) of dust shells and was, ultimately, used in order to ﬁt a model for a shell of ionized gas. Results. For all Cepheids, we ﬁnd a continuum IR excess increasing up to approximately − 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 signiﬁcant 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 ﬁrst 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 (cid:39) 15% of the star radius, a mass of 10 − 9 − 10 − 7 M (cid:12) and a temperature ranging between 3500 and 4500K. Conclusions. The presence of a thin shell of ionized gas around Cepheids must be tested with interferometers operating in the visible or mid-IR, or using radio telescopes. The impact of such CSEs of ionized gas on the PL relation of Cepheids also calls for further investigation.


Introduction
Cepheids have been the keystone of distance scale determination in the Universe for a century because of the direct correlation between their pulsation period and their luminosity as per the Leavitt law (Leavitt 1908;Leavitt & Pickering 1912), also referred to as the period-luminosity relation (hereafter PL relation).
The recent determination of one percent precision for the Large Magellanic Cloud distance (Pietrzyński et al. 2019) has led to a new value for the Hubble constant H 0 (Riess et al. 2019). 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 a PL relation for Cepheids whose uncertainties on both zero point and slope today stand as one of the largest contributors to the error on H 0 (Riess et al. 2019). One possible bias could be due to IR excesses from CSEs, such as those discovered using near-and mid-infrared interferometry around nearby Cepheids Mérand et al. 2006). Indeed, if the brightness of CSEs is found to be significantly different in the Milky Way, SMC, LMC, as well as 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 longbaseline interferometry in the K-Band with VLTI and CHARA Mérand et al. 2006), and four Cepheids CSEs have been observed in the N band with VISIR and MIDI (Kervella et al. 2009;Gallenne et al. 2013), as well as one with NACO in the near-IR (Gallenne et al. 2011). The presence of a motionless Hα absorption component using high-resolution spectroscopy around l Car (Nardetto et al. 2008) has also been attributed to a CSE, and, recently, a CSE was detected in the visible domain with the VEGA/CHARA facility around δ Cep (Nardetto et al. 2016). These various studies determine 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 mediumand 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, nor do we understand their nature or their characteristics (density and temperature profiles, chemical composition...). This paper is aimed at building a phase-dependent spectral energy distribution (SED) of a sample of Cepheids from visible to mid-IR wavelengths and comparing 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. 2 using photometric and Spitzer observations in various bands, and we study qualitatively far-infrared images from Herschel. In Sect. 3 we correct the spectra from interstellar foreground absorption along the line-of-sight and seek residuals at 9.7 µm that could be attributed to a dusty CSE. In Sect. 4, we use the radiative transfer code DUSTY (Ivezic et al. 1999) to model the IR excess continuum and to show that such continuum cannot be attributed to a dusty CSE. In Sect. 5, we propose, instead, a model of CSE composed of ionized gas in order to reproduce the observed IR excess. Our results are discussed in Sect. 6 and we present our conclusions in Sect. 7.

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 to 38 µm, makes it the best tool for studying dust features and IR excess continuum emission. We selected the five Cepheids with Spitzer low resolution spectra and the highest signal-tonoise 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 and 160 µm to study larger cold environments of dust. These space-based observations have the advantage of avoiding biases caused by the Earth atmosphere, meaning the entire flux is preserved and not hidden by the ozone band, which is essential for the study of 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 (Mérand et al. 2015). In order to compute synthetic photometry to match what is drawn from the dataset, SPIPS uses a grid of ATLAS9 atmospheric models 1 (Castelli & Kurucz 2003) with solar metallicity and a standard turbulent velocity of 2 km s −1 . SPIPS models are available for η Aql (Mérand et al. 2015), ζ Gem (Breitfelder et al. 2016) and RS Pup (Kervella et al. 2017). We updated them with new datasets, including Gaia photometry 1 http://wwwuser.oats.inaf.it/castelli/grids.html (Gaia Collaboration 2018), and we also applied the SPIPS algorithm to SU Cyg and V Cen.  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 observational data. While the fit is satisfactory from the visible domain to the near infrared, a disagreement results for wavelengths larger than about 1.2 µm. For λ > 1.2 µm, while 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, for λ < 1.2 µm α(λ − 1.2) β , for λ 1.2 µm (1) 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. (1) was adopted originally in SPIPS modeling because it provides a reasonable description of the IR excess without impacting the results of the fit. In Fig. 1, 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.05 mag) from Fig. 1, we conclude that the IR excess of Cepheids is not time-dependent or only time-dependent in the slightest.
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 (Scowcroft et al. 2016) and this is, indeed, probably what is seen in IRAC I2 Spitzer bands for ζ Gem and V Cen (see Figs. A.1 and A.2,respectively). In these bands, the relation 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 dynamic 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. 1), and is considered in the following part of the analysis.
In addition, the photometric bands longward of 5 µm are problematic for two reasons. First, their cycle coverage is poor, meaning there is only one data point in each band with an undetermined phase, as is the case for η Aql, ζ Gem and RS Pup at λ > 5 µm (see, for instance, the A MSX band of RS Pup in Fig. A.3). 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 have not been considered are indicated by orange bars in Fig. 1.

Spitzer data
The spectroscopic observations were made with the InfraRed Spectrograph (IRS; Houck et al. 2004) on board the Spitzer telescope (Werner et al. 2004) and the full spectra were retrieved  from the CASSIS atlas (Lebouteiller et al. 2011). CASSIS identified the five stars as point-like sources, therefore we retrieved the best flux calibrated spectrum obtained from the optimal extraction method (Lebouteiller et al. 2010). The Spitzer dataset of this paper is presented in Table 1. 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 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 2. We derive the IR excess of each star in the sample at the specific phase of Spitzer using:   Notes. T eff and log g are used to interpolate a ATLAS9 atmosphere model. The limb-darkened angular diameter θ is then used to derived the observed flux. Temperature and angular diameter are in Kelvin and milliarcsecond respectively. Uncertainties on T eff and θ are provided by SPIPS, whereas uncertainty on log g is arbitrarily set to 10%.
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 3 and the ∆mag values we ultimately obtain are indicated by green dots and red uncertainties in Fig. 2. The same method is applied to all the photometric bands of observations, meaning that 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. 2). Please note that in Fig. 2, 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 (Lebouteiller et al. 2011;Sloan et al. 2015) are not excluded, we find a discontinuity at 14 µm between SL and LL spectra, indicated by a dashed line in Fig. 2. According to CASSIS atlas, a residual offset appears in the case of slightly extended sources (Lebouteiller et al. 2011) because more flux enters the larger LL aperture. Indeed the two Spitzer detectors have indeed different fields of view (FOV, see Table 2) with different orientations on the sky (see Fig. 3), which explains the different level of flux measured. As an indication, the FOV of both detectors are overlaid on Herschel images (see Fig. 3 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 (Westerlund 1961;McAlary & Welch 1986) so 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 the CASSIS pipeline is reliable and the emission of the environment is homogeneous between SL and LL apertures.
In addition, order mismatches may appear because of telescope pointing accuracy and ensuing 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. We also note that the Spitzer data at larger wavelengths than 30 µm were not considered due to their extremely large scopes of 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, prior to the correction of different Spitzer discontinuities (see Fig. 6 for a detailled analysis of corrected IR excess). The agreement between IR excesses from SPIPS models and Spitzer data demonstrates 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 that obtained by Schmidt (2015), where no IR excess was found for short-period classical Cepheids (see Sect. 6.1). We, however, find that the ad hoc analytic laws of SPIPS overestimate the IR excess in the Spitzer wavelength domain (see Fig. 2).
Third, RS Pup, SU Cyg and V Cen present unambiguous silicate absorption bands at 9.7 and 20 µm. The latter absorption is known to usually peak at 18 µm, however it can be shifted to a longer wavelength depending on various conditions, such as temperature or mineral proportions indicated by silicates (Koike et al. 2006;Henning & Mutschke 1997). These silicate absorptions are likely due to the presence of interstellar clouds in 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 that is due to the environment of the Cepheids. Then we derived the IR excess at 70 and 160 µm through comparisons with the atmospheric model in the same way as in Eq.
(2). Observations are presented in Fig. 3 and the derived photometry in Table 4. From these data we observe important IR excesses at 70 and 160 µm, which are associated with cold environments of temperatures of about 40 and 20 K, respectively, according to Wien's displacement law (λT = 2900 µmK). Cold dust material around RS Pup was already observed, while the large cloud observed around V Cen (see Fig. 3d) at 160 µm could be related to a star forming region since V Cen is likely a member of the open cluster NGC 5662 (Turner 1982;Claria et al. 1991;Anderson et al. 2013). Herschel is sensitive not only to emission from the star itself and a potential CSE, but also to that of cold extended emission from the interstellar medium.
A47, page 5 of 19 A&A 633, A47 (2020) Notes. F 70 and F 160 are the flux in Jansky at 70 and 160 µm, respectively integrated within a 1 circle centered on the star. ∆mag 70 and ∆mag 160 are the corresponding IR excess derived following the definition of Eq.
(2) but for Herschel observations. Thus, we do not consider the Herschel photometry in the rest of this 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 strong silicate absorptions and extended emission in the Herschel data at the same time (see Figs. 3a and d), while η Aql and ζ Gem present no obvious silicate absorption and weaker interstellar environment (Figs. 3b and c).

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. 2). Such a level of 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

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. Thus, we come to the following relation: All values are presented at magnitudes of 9.7 µm. We then determine A Spitzer 9.7 (Sect. 3.1) and A ISM 9.7 (Sect. 3.2) in order to estimate the CSE emission at 9.7 µm (Sect. 3.3) and we extend our correction to the whole wavelength range of Spitzer (Sect. 3.4).

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. 4), and we interpolate the excess continuum value at 9.7 µm (blue dot in Fig. 4). We then subtract 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.
A47, page 6 of 19 V. Hocdé et al.: A thin shell of ionized gas as the explanation for IR excess among classical Cepheids  is set to 0.000 ± 0.002.

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 2 (hereafter DDOD; Fernie et al. 1995) which 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 (Capitanio et al. 2017) using the Bayesian inversion method to a wide color excess dataset and parallaxes; and (4) the 2D map from Schlegel et al. (1998) (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 Drimmel & Spergel (2001). We present these E(B-V) values in Table 5. The best agreement is obtained when comparing SPIPS and DDOD values, which are based on two enitrely 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, they cannot be applied to 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 (Mérand et al. 2015). 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 of A v = R v E(B-V) with a ratio of total-toselective extinction of R v = 3.1, which corresponds to a diffuse ISM along the line of sight (Savage & Mathis 1979). Then we used the relation A v /τ 9.7 = 18.5 (Roche & Aitken 1984) 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 (Kervella et al. 2009. Hence we treated RS Pup as a special case since the extinction towards dense clouds is different from diffuse ISM (Whittet et al. 1988;Indebetouw et al. 2005;Flaherty et al. 2007;Chiar et al. 2007;van Breemen et al. 2011). 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 McClure (2009). 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 ultimately be summarized on the basis of two equations, referring to stars in diffuse ISM Eq. (4) and for RS Pup Eq. (5): A ISM 9.7 = 1.086 A ISM 9.7 = 1.086 A&A 633, A47 (2020) Notes.  Notes. The silicate absorption due to interstellar extinction A ISM 9.7 is calculated (using Eq. (4) The final values of A ISM 9.7 we have considered are listed in Table 6, along 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. (3). The results are illustrated in Fig. 5. 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, along with 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. (4). If we assume, instead, the presence of dense clouds in the line of sight using Eq. (5), we obtain a silicate emission of −0.091 ± 0.013 mag (see Fig. 5). From this figure (and Table 5), it is difficult to conclude whether there is a dusty CSE of silicate around the Cepheids in the sample or not, particularly because our method is sensitive to the used Eq. (4) or (5) when correcting the silicate ISM absorption (case of V Cen for instance), as well as to various sources of uncertainties, especially in the E(B-V) estimate. We also assume that the IR excess does not vary in time, as has been suggested by the SPIPS analysis referred to in Sect. 2.2. Still, using this method, we find residual dust emission at 9.7 µm for the long-period Cepheids in the sample (η Aql, ζ Gem and RS Pup). (4), than the E CSE 9.7 value is larger (orange bar in the figure). CSE emission appears when E CSE 9.7 < 0.
3.4. Extending the correction of the silicate absorption from 9.7 µm to the whole wavelength range of Spitzer We correct the entire scope of Spitzer observations by subtracting a synthetic interstellar medium composed of silicates. Since we assumed an averaged ISM temperature of 20 K, the dust emission is negligible in the Spitzer wavelength range according to Wien's law. Thus, we must 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 adopt the following expression for λ between 5 and 30 µm A ISM λ ∝ κ λ = C abs λ (a)πa 2 n(a)da.
We first derive C abs λ using complex refractive index for silicates from Draine & Lee (1984;hereafter DL84), assuming an uniform A47, page 9 of 19 A&A 633, A47 (2020) distribution of ellipsoidal shapes given by Bohren & Huffman (1983). Then we derived the absorption coefficient κ λ by taking into account a standard grain size distribution n(a) ∝ a −3.5 (Mathis et al. 1977). Finally we normalize A ISM λ using its specific value A ISM 9.7 at 9.7 µm that we already derived in Sect. 3.2. For all the stars in the sample, we plot the IR excess (found to be constant with the pulsation phase, see Sect. 2.2) corrected from the ISM silicate absorption in Fig. 6. 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 so, most of the carbon in their envelopes is locked in CO molecules. In the condensation sequence described by Gail & Sedlmayr (1999) corundum (Al 2 O 3 ) is expected to form first because of its high condensation temperature of about 1400 K in typical pressures encountered in circumstellar shells. For lower temperatures, Al 2 O 3 is depleted and silicates such as gehlenite (Ca 2 Al 2 SiO 7 ) and forsterite (Mg 2 SiO 4 ) can form. All these components present emission components mostly observed in the N band (8-13 µm, see Fig. 7), 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 µm. Figure 8 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 carried out by Schmidt (2015) in finding that CSE made of silicate dust cannot explain IR excess of 132 classical and type-II Cepheids. However, since the opacity of iron exhibits no particular feature, but rather a continuum (see Fig. 7), we investigate whether a warm dust envelope of iron could explain the IR excess continuum observed with the SPIPS and Spitzer dataset. The CSEs were modeled using DUSTY (Ivezic et al. 1999), which solves the radiative transfer equations in a dusty environment. The method is based on a self-consistent equation for spectral energy density, including dust scattering, absorption, and emission. We present two CSE models with different inner shell temperatures. Iron is a component in oxygen-rich star mineralogy with a condensation temperature of ≈1000 K for circumstellar pressures. Hence, the cold model takes into account the typical condensation temperature of iron (≈1000 K) 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 (Mathis et al. 1977), 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 3 in Sect. 2.3) as central source radiation in the DUSTY computation. The results of the computation are summarized in Table 7 and the computed IR excess is presented Fig. 9. 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 the 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 be the cause of IR excess continuum for Cepheids in the 1-30 µm range.

The IR excess from a thin shell of ionized gas
The next step in our investigation considers whether an ionized, spherical gas shell can simultaneously explain at the near-(∼1 µm < λ < 5 µm) and mid-(5 µm < λ < 30 µm) IR excess of Cepheids. As a first step, this study is conducted only at the specific pulsation phase of Spitzer data. We discuss the possible time variability of such a shell in Sect. 6.3. We consider the emission of a thin gas shell around the star with its density and temperature constant for the sake of A47, page 10 of 19 V. Hocdé et al.: A thin shell of ionized gas as the explanation for IR excess among classical Cepheids Notes. T in and θ in are the inner temperature and diameter in milliarcseconds respectively. τ 0.55 is the opacity at 0.55 µm. Mass lossṀ derived by DUSTY is also indicated. (a) Fitted parameters. simplicity. The shape of the mid-IR excess, saturating to a constant flux ratio at large wavelengths (see Fig. 6), suggests an opacity source increasing with wavelength. We used the freefree 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 opacity sources is given by (e.g. Rybicki & Lightman 2008) where h, c, and k hold their usual values, γ is the degree of ionization (between 0 and 1), T s the temperature of the shell, m H the hydrogen mass, and g ff and g bf are the free-free and bound-free Gaunt factors, respectively. These factors were computed mainly from approximation formulas given by Brussaard & van de Hulst (1962), Hummer (1988), and references therein. As an example, we present a typical shape of the absorption coefficient κ λ computed for V Cen in Fig. 10 using the shell parameters from Table 8. We compute the SED of the star plus the gas shell as described in Appendix B, taking into account the latter absorption coefficient κ λ . In order to match the SPIPS photometries, in addition to the corrected Spitzer spectra presented in Sect 3, we perform a χ 2 fitting using the Levenberg-Marquardt method. We use bandpass filters 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, that is, ∆m = 0 for λ < 1.2 µm (see Sect. 2.2), it is necessary to relax this assumption to allow the data to present deficit or excess in the visible domain depending on the physical behaviour of the ionized shell. This assumption is often considered equivalent to the suggestion for 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 four parameters: three parameters from the gas shell (the ionized shell mass γM s , its temperature T s , and radius R s ) plus one parameter corresponding to the IR excess offset equal to ∆m 0 for λ < 1.2 µm. Since Spitzer data have a higher statistical weight than SPIPS data, we first considered the different parameters by hand as first guesses to find the main ones for our study. The results are presented in Table 8 and in Fig. 11. In Fig. 11 we observe discontinuities at short wavelengths, from visible to near-IR, which are due to bound-free opacities (see Fig. 10). 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 4500 K and an envelope thickness of 15% of the radius of the star, while γM s ranges from 10 −9 to 10 −7 M . The IR excess correction offsets are 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 models computed present absorption in the visible domain (see Fig. 11), thus correction offsets have positive values to allow a deficit in the visible due to ionized shell absorption.

The IR excess and dust environment
We estimate the infrared excess of Cepheids using the SPIPS algorithm (λ < 5 µm) and Spitzer data (5 < λ < 30 µm). The observed IR excess presented in Fig. 6 is thus calculated at the specific phases of Spitzer and makes the assumption that there is no excess 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 exhibit a continuum which is consistent across the wavelength range for all stars. It ranges in the data from about 2 to 30 µm and corresponds to differences in magnitude 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 difference (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. 3) 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 a few percents, which is, on average, ten times less than what was predicted by other studies based on interferometry (Gallenne et al. 2012(Gallenne et al. , 2013, and consistent with Schmidt (2015) who found no silicate emission for a large sample of Cepheids.
Third, using photometric bands and the same basic approach, but without considering the pulsation of the Cepheids, Schmidt (2015) found that 21 of 132 classical and type 2 Cepheids in his sample demonstrate a clear or weak IR infrared excess. Moreover, these 21 Cepheids have periods larger than 11 days. Conversely, we have 4 short-period 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 delve deeper into the analysis.
Four, we find a large, cold 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 an 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 affect estimations of the PL relation when using instruments such as the JWST, as well as forthcoming mid-infrared instruments set on 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 4500 K, 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 was carried out using the Very Large Array (VLA) at 5 GHz on η Aql and four other classical Cepheids (Welch & Duric 1988). 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 5 GHz (≈20% above the star continuum) which is below the upper limit on flux density of ≈100 µJy estimated by Welch & Duric (1988).
On the other hand, interferometric observations have resolved CSEs around Cepheids. The first detection was reported around l Car ) followed by δ Cep and Polaris (Mérand et al. 2007). These CSEs were modeled with a ring at a distance of 2-3 R , that is, close to the star, in a region where the temperature is high enough (>2000 K) to prevent dust condensation (Gail & Sedlmayr 1999). Thus, these observations are more likely explained by a shell of partially ionized gas. This finding in detection by interferometry appears to support a widespread phenomenon among 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 occur when the atmosphere is compressed during its infalling movement or because of shock waves dynamics (Breitfellner & Gillet 1993a,c,b;Fokin et al. 1996). 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 (Gillet 2014). In addition, analytical works by Neilson & Lester (2008) have shown that mass loss is enhanced by pulsations and shocks in the atmosphere. Our study suggests that this loss in mass could be due to the effects of partially ionized gas.
The model of the ionized gas shell could also be linked to the chromospheric activity of Cepheids. Sasselov & Lester (1994) 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 (2900 A) of five Cepheids using the International Ultraviolet Explorer instrument (Schmidt & Parsons 1984) revealed extended dynamical chromospheres up to a 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-100 km s −1 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 (Engle et al. 2017).
In parallel, it is interesting to compare Cepheids to very longperiod Mira stars for which a radiosphere near 2 R attributed to free-free emission has been reported (Reid & Menten 1997).

Limit of the model and perspectives
We have found relatively low temperatures for hydrogen ionization (between 3500 and 4500 K see Table 8) and only a small fraction of the gas should be ionized according to the Saha equation. In particular, for gas temperatures below 4000 K, 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 lowmetallicity Cepheids in the Magellanic clouds. Thus, it would be interesting to theoretically quantify the impact of metallicity on the shells of ionized gas, and the PL relation. Moreover, our model does not take into account temperature or density gradients in the star's atmosphere, nor, in particular, compression 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 time-dependent at all, or only slightly. Evidence 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 (Scowcroft et al. 2016). As we have Spitzer data only at a specific phase of pulsation for each Cepheid, we cannot firmly conclude a figure for the timedependency 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 in time. In conducting the test for η Aql, we find a rather stable relative size for the shell ( 15 ± 2%) and a temperature variation from about 3500 to 4500 K, 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 also be lower by few percents, but this can be compensated by the actual size of the shell, and would have little impact within the SPIPS fitting. Thus, a spatial and chromatic analysis of the shell, including interferometric constraints in all available bands, 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, ultimately, to survey the impact on the PL relation.

Conclusion
1. For the five Cepheids, we report a continuum IR excess increasing up to approximately −0.1/−0.2 magnitude at 30 µm, which cannot be explained by a hot or cold dust model of CSE. 2. Within the limits of our assumptions, we do not assert a firm conclusion regarding the presence of CSE emission in the N-band, but we assume it is likely to prove weak (> −0.1mag) according to our results. 3. We demonstrate for the first time that the IR excess of Cepheids can be explained by a free-free emission of a thin shell of ionized gas with a thickness of a 8-17% star radius, an ionized mass of 10 −9 −10 −7 M , and a temperature of 3500-4500 K. In this simple model, density and temperature have a constant radial distribution. 4. The presence of a thin shell of partially ionized gas around Cepheids must be tested with interferometers operating in the visible domain, in the mid-IR, or in the radio domain.
The impact of such CSEs of ionized gas on the PL relation of Cepheids also requires futher investigation.   The plots are organized as follows: pulsational velocity, effective temperature and 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, the projection factor set to p = 1.270 is indicated, along with 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 Castelli & Kurucz (2003) with solar metallicity and a standard turbulent velocity of 2 km s −1 . 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 µm, only one data point is obtained without information on the phase. Hence, it is represented by a hori-zontal 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 (Scowcroft et al. 2016). This effect is observed in the Spitzer I2 IRAC dataset of ζ Gem and V Cen (see Figs. A.1 and A.2). 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 +/−200 K 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 a reliable (for example see E(B-V) value compared with those provided in the literature in Table 5).

Appendix B: 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. B.1) according to the method described in Panagia & Felli (1975). We assume a constant shell temperature T s and density ρ s . The gas opacities (see Eq. (7)) can be written under the form κ(λ, T s ) = ρ 2 s χ λ (T s ). According to Fig. B.1 we have to take into account two cases corresponding to the impact parameter p respectively larger and smaller than the stellar radius R * .
For p R * , taking into account the symmetry of the problem the optical depth along the ray is given by Similarly, for p R * we have For p R * , the specific intensity is given by since T s is assumed to be constant along the ray. For p R * , both the shell and the stellar photosphere contribute to the specific intensity I λ (p) = B λ (T s ) 1 − e −τ λ (p) + I * λ e −τ λ (p) , (B.5) 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 with F * λ = π R * d 2 I * λ . Note that in the particular case where we have I * λ = B λ (T * ) with T s = T * , Eq. (B.6) can be integrated analytically to give with τ * λ defined as 2 ρ 2 s χ λ (T s ) R * R out R * 2 − 1.
The corresponding magnitude excess is given by ∆mag = 2.5 log