© ESO, 2011

1. Introduction

For almost a century, Cepheid stars have been used as distance indicators thanks to the well-known period-luminosity relation (P-L), also called Leavitt’s law. A good calibration of this relation is necessary to obtain an unbiased estimate of the distance. After the detection of circumstellar envelopes (CSEs) around many Cepheids in the near- and mid-infrared (Kervella et al. 2006; Mérand et al. 2006, 2007; Kervella et al. 2009a), it was shown that these CSEs could bias the angular diameter determination using near-infrared interferometry (Mérand et al. 2006). A circumstellar envelope could have an impact on the distance estimate if we do not take into account its flux contribution. It could lead to a bias both in interferometric and photometric measurements and therefore influence the accuracy of the estimated distance. Because Cepheids are used as standard candles to estimate distances in the universe, it is important to understand the properties of the CSEs around this kind of stars.

We report new observations with VLT/NACO in the IB_2.18 and NB_1.64 filters using a fast imaging mode in order to characterize the extended emission around RS Puppis (HD 68860) at sub-arcsecond angular resolution (≈0.01 pc at RS Pup distance). Data and reduction methods are described in Sect. 2. In Sects. 3 and 4 we present two different data processing analyses: the shift-and-add method, and a method based on a statistical analysis of the speckle noise to enable the detection of the CSE in the speckle cloud surrounding the star’s image.

This paper is focused on RS Pup’s envelope, which has been known for many years (Westerlund 1961; Havlen 1972) for it is very extended: it has a spatial extension of around 2′ in V band, which corresponds to about one parsec at RS Pup’s distance (about 2 kiloparsecs; Kervella et al. 2008). Its origin is not precisely known; it is especially unclear whether it is made of material ejected by the star or of a material that predates the formation of the star. Its actual shape has been widely discussed: discrete nebular knots have been invoked by Kervella et al. (2008) to estimate the distance using light echoes, whereas Feast (2008) argued that a disk-like CSE will bias this kind of estimate. Alternately, a bipolar nebula has been considered by Bond & Sparks (2009), making a ratio between two images of different epochs. We were able to spatially resolve the CSE close to the star. The morphology of the CSE is discussed in Sect. 4.3.

2. Observations and data reduction

RS Pup was observed with the NACO instrument installed at the Nasmyth B focus of UT4 of ESO VLT. NACO is an adaptive optics system (NAOS, Rousset et al. 2003) and a high-resolution near IR camera (CONICA, Lenzen et al. 2003), working as imager or as spectrograph in the range 1–5 μm. We used the S13 camera for our RS Pup observations (FoV of 13.5″  ×  13.5″ and 13.26 ± 0.03 mas/pixel, Masciadri et al. 2003) with the narrow band filter NB_1.64 (λ₀ = 1.644 μm with a width Δλ = 0.018 μm) and the intermediate band filter IB_2.18 (λ₀ = 2.180   μm with a width Δλ = 0.060 μm). We chose the cube mode in order to apply the shift-and-add technique. We chose windows with 512 × 514 pixels and 256 × 258 pixels respectively for the narrow and the intermediate band filters. Each cube (10 for RS Pup and 4 for the reference star) contains 460 images at 1.64 μm and 2000 at 2.18 μm. Each frame has the minimum integration time and readout type for this mode, i.e. 109 ms and 39 ms respectively for the 1.64 and 2.18 μm cubes. With these DITs and filters we are in a quasi-monochromatic regime with a very good speckle contrast (Roddier 1981). The large number of exposures allows us to select the frames with the brightest image quality.

Data were obtained on 2009 January 7, with the atmospheric conditions presented in Table 1. The PSF (point spread function) calibrator stars were observed just before or immediately after the scientific targets with the same instrumental configurations. Each raw image was processed in a standard way using the Yorick¹ language: bias subtraction, flat-field and bad pixel corrections. The negligible sky background was not subtracted.

To compare the AO system behavior and artefacts on RS Pup, we used the star Achernar and its PSF reference star observed the same night. These data were also acquired in cube mode with the same camera in the NB_1.64 and NB_2.17 (λ₀ = 2.166   μm with a width Δλ = 0.023   μm) filters. Only one cube was acquired per filter, containing 20 000 frames for Achernar and 8000 for the calibrator. The window detector size is 64 × 66 pixels in both filters with an integration time of 7.2 ms. The slightly different second filter and the shorter exposure time compared with RS Pup is not critical for our analysis because we are interested in flux ratios. The smaller window size cuts the AO PSF wings, but we will show in Sect. 3 that this does not affect our conclusion significantly.

Once the pre-processing was over, we carried out a precentering as well as a sorting according to the maximum intensity of the central peak to reject the images for which the correction was the least efficient. On the one hand, we applied the shift-and-add technique to obtain the best angular resolution possible. On the other hand, we applied a statistical study aimed at extracting a possible diffuse (i.e. completely resolved) component inside of the speckle cloud of the central core (the star). These methods are explained in Sects. 3 and 4.

3. Shift-and-add process

In this section we use the shift-and-add method and then compare the RS Pup images with another star observed on the same night for which we do not expect an envelope.

3.1. RS Pup

Both stars, RS Pup and HD 74417, are unresolved by the telescope: diffraction limits for an 8.2 m UT are 55 mas and 41 mas in filters IB_2.18 and NB_1.64 respectively, while angular diameters are approximately one milli-arcsecond (1 mas) for both stars.

The shift-and-add method was originally proposed by Bates & Cady (1980) as a technique for speckle imaging. Applied to NACO’s cube mode, it enhances the Strehl ratio by selecting the frames that are the least altered by the atmospheric turbulence. It also reduces the halo contribution when used with adaptive optics by selecting the best AO-corrected frames. The efficiency of the cube mode vs. the standard long exposure mode is discussed in Kervella et al. (2009b). Girard et al. (2010) showed for the RS Pup’s reference star at 2.18 μm that by using this method we have an absolute gain of 9% in Strehl ratio.

Our processing steps are as follows: we select the 10% best frames according to the brightest pixel (as a tracer of the Strehl ratio) in our 10 data cubes (4 for the calibrator). We then spatially resample them by a factor 4 using a cubic spline interpolation and co-align them using a Gaussian fitting on the central core, at a precision level of a few milliarcseconds (this method is described in detail in Kervella et al. 2009a). Each cube is then averaged to obtain ten final average images (4 for the calibrator). We finally compute the mean of the 10 average images of RS Pup (or the 4 for the calibrator).

The result of this method is shown in Fig. 1 for the two filters. The two upper frames are HD 74417 and RS Pup as seen through the IB_2.18 filter, while the two lower are in the NB_1.64 filter. The spatial scale in all images is 3.3 mas/pix and we chose a logarithmic intensity contrast with adjusted lower and higher cut-off to show details in the halo. On the 1.64 μm image of RS Pup compared with the HD 74417 image we can see that there is a resolved circumstellar emission up to  ~1″ from the Cepheid. This feature is fainter in the 2.18 μm filter. In both filters, because of this emission, diffraction rings are less visible around RS Pup than around HD 74417.

Fig. 1 Final frames from shift-and-add process. In the two upper panels the PSF calibrator star HD 74417 and RS Pup are shown in the IB_2.18 filter. The bottom panels show the same stars in the NB_1.64 filter. The scale intensity is logarithmic.

For each final co-added frame (i.e. the image shown in Fig. 1) we computed a ring median and normalized to unity. The resulting curves are presented in Fig. 2 for two filters. The error bars are the standard deviation of values inside rings and are plotted only for few radii for clarity. We notice a significant difference between the two curves, more pronounced in the NB_1.64 filter.

3.2. Achernar

We carried out the same processing for Achernar and its reference star δ Phe (HD 9362). These two stars are unresolved by the UTs (angular diameters are approximately 2 mas). With this window size (64 × 66 pixels) the field of view is four times smaller than in the RS Pup’s frames, and even if we can also see some Airy rings in all the images, the background is dominated by the AO halo.

We plot in Fig. 2 the radial profiles of the stacked images, normalized to 10 (for clarity in the graph).

Fig. 2 On the left are presented the radial profiles at 2.18 μm for RS Pup, Achernar and their reference stars. On the right the same at 1.64 μm. Radial profiles were computed using the median value on one-pixel thick circular rings. The gray color denotes Achernar and δ Phe profiles, shifted by a factor 10 for clarity. The scale is logarithmic for both axes.

3.3. Halo fitting

We now investigate the presence of a circumstellar envelope by studying the radial profile of the stacked images we obtained for both RS Pup and its PSF calibrator and compare them with the observations of Achernar. The main idea behind this study is to compare the variations of encircled energy (the ratio of the core’s flux to the total flux) to the expected variations caused by atmospheric condition changes. If RS Pup has indeed a circumstellar envelope, its encircled energy will be significantly lower than that of our reference stars.

We use an analytical function to fit the residual light out of the coherent core (i.e. the halo). We used a turbulence degraded profile from Roddier (1981)${\mathit{I}}_{\mathrm{halo}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathit{f}\frac{\mathrm{0.488}}{{\mathit{\rho }}^{\mathrm{2}}}{\left[\mathrm{1}\mathrm{+}\frac{\mathrm{11}}{\mathrm{6}}{\left(\frac{\mathit{r}}{\mathit{\rho }}\right)}^{\mathrm{2}}\right]}^{\mathrm{-}\frac{\mathrm{11}}{\mathrm{6}}}\mathit{,}$(1)where the fitting parameters are the FWHM ρ and f proportional to the flux of the envelope. In this case, the parameter ρ corresponds to the equivalent seeing disk for long exposure. NACO’s resulting core, even for an unresolved source, is fairly complex and suffers from residual aberrations, this is why we decided to fit the halo at a large radial distance of the core. We will later define the core as the residual of the fit of the halo function to our data.

We show the resulting fits in Fig. 4. They were applied for r > 0.26″ at 2.18 μm and 2.17 μm while the range was for r > 0.2″ at 1.64 μm for RS Pup (and its calibrator) and r > 0.26″ for Achernar (and its calibrator). The black lines are our data (ring-median), while the gray lines represent the adjusted curves. The vertical scale is logarithmic and normalized to the total flux. We also evaluated the normalized encircled energy (E_enc), i.e. the ratio of the flux in the core (i.e. within 1.22λ/D) to the total flux. We can already see from these curves that the flux in the halo is different in the images of RS Pup and Achernar. These differences can be linked to the emission of the CSE around RS Pup or a modification of the seeing between the two observations. This point will be discussed bellow with a statistical study of the encircled energy.

For RS Pup in the NB_1.64 filter we found from our fit comparing the two FWHMs, i.e. (ρ_RS   Pup − ρ_{HD   74417})/ρ_{HD   74417}, a change of halo size of 11% and encircled energies E_enc(RS   Pup) = 0.327 and E_enc(HD   74417) = 0.452 (see Table 2). This yields an encircled energy ratio ε ~ 72%. The local DIMM station registered a seeing variation of around 26%. Different values can be explained by the station location, which is not inside the telescope dome, but a few tens of meters below, rendering it more susceptible to ground layer effects that are known to affect the seeing (Sarazin et al. 2008). As a result the UT seeing is consistently better than the one measured by the DIMM. For IB_2.18 we show the fitting curves in Fig. 4 and the resulting fit parameters in Table 2. Our fit gives a seeing variation of around 9%, while the DIMM measurements give a 29% variation, leading to the same conclusion as before about the accuracy of the DIMM station seeing estimates. The encircled energies measured in IB_2.18 are E_enc(RS   Pup) = 0.417 and E_enc(HD   74417) = 0.517, giving an encircled energy ratio ε ~ 80%.

Fig. 3 Histogram of the encircled energy computed for 10% best frames. The solid lines represent a Gaussian fit.

We evaluated the uncertainties according to the variations of the encircled energy in the cubes. Because we already selected the 10% best frames, we kept these samples and plotted the histogram of the encircled energy to check the distribution function (Fig. 3). At the first order, the E_enc distribution can be approximated by a Gaussian distribution. We then computed the relative standard deviation $\mathit{\sigma }\mathrm{\left(}\mathit{E}{}_{\mathrm{enc}}\mathrm{\right)}\mathit{/}\overline{{\mathit{E}}_{\mathrm{enc}}}$ for both stars. We assessed the uncertainty of the encircled energy ratio ε with $\frac{\mathit{\sigma }\mathrm{\left(}\mathit{\epsilon }\mathrm{\right)}}{\mathit{\epsilon }}\mathrm{=}\sqrt{\left[\frac{\mathit{\sigma }\mathrm{\left(}\mathit{E}\mathrm{enc}\mathrm{\right)}}{\overline{\mathit{E}\mathrm{enc}}}\right]\begin{array}{}\mathrm{2}\\ \mathrm{RS} \mathrm{Pup}\end{array}\mathrm{+}\left[\frac{\mathit{\sigma }\mathrm{\left(}\mathit{E}\mathrm{enc}\mathrm{\right)}}{\overline{{\mathit{E}}_{\mathrm{enc}}}}\right]\begin{array}{}\mathrm{2}\\ \mathrm{HD} \mathrm{74417}\end{array}}\mathit{,}$(2)We assume that the relative standard deviation of the encircled energy is a good approximation of the seeing variations. We measured a relative standard deviation $\mathrm{\left[}\mathit{\sigma }\mathrm{\right(}\mathit{E}{}_{\mathrm{enc}}\mathrm{)}\mathit{/}\overline{\mathit{E}{}_{\mathrm{enc}}}{\mathrm{]}}_{\mathrm{RS} \mathrm{Pup}}\mathrm{~}\mathrm{8}\mathrm{\%}$ at 2.18 μm and  ~10% at 1.64 μm while $\mathrm{\left[}\mathit{\sigma }\mathrm{\right(}\mathit{E}{}_{\mathrm{enc}}\mathrm{)}\mathit{/}\overline{\mathit{E}{}_{\mathrm{enc}}}{\mathrm{]}}_{\mathrm{HD} \mathrm{74417}}\mathrm{~}\mathrm{4}$% and  ~6% respectively. We also evaluated the statistical uncertainties with a bootstrap method that gave σ ~ 0.5% for both stars. These statistical uncertainties do not take into account the seeing variations. We therefore chose Eq. (2) as a reliable estimate of the final uncertainty. This gives σ(ε)/ε ~ 9% and  ~12% respectively in IB_2.18 and NB_1.64.

Fig. 4 Halo fit for all filters with logarithmic ordinates. The fluxes are normalized to the total flux. Left panel: halo fit from the radial profile of images in the two filters for RS Pup. Right panel: halo fit from the radial profile of images in the two filters for Achernar. The fit in IB_2.18 was for r > 0.26″ and in NB_1.64 for r > 0.2″ for RS Pup and for r > 0.26″ for Achernar. The dashed lines denote the reference stars, while the thick lines denote the Cepheid and Achernar. The fitted curves are plotted in gray.

For Achernar, the fit gives a seeing variation  ~3% both in NB_1.64 and IB_2.17. The encircled energies measured in IB_2.18 are E_enc(Achernar) = 0.495 and E_enc(δ   Phe) = 0.505, giving an encircled energy ratio ε ~ 98% (see Table 2). In NB_1.64 we have 0.374 and 0.376 respectively, which gives ε ~ 99%. The encircled energies are nearly identical, giving ratios close to 100%. Using Eq. (2) as previously, we found a relative standard deviation  ~15% in IB_2.17 and  ~11% in NB_1.64. Because Achernar does not have an extended circumstellar component in NB_1.64 and IB_2.17, we attribute these variations to the Strehl variations and the AO correction quality. The encircled energy ratio seems stable between the two stars at the 2% level in IB_2.18 (see Table 2).

We can notice that the normalized curves of Achernar and its calibrator almost overlap, and this is what is expected for unresolved stars without envelope taken in same atmospheric and instrumental conditions. It is not the case for RS Pup (see Figs. 2 and 4). The additional emission coming from the envelope of the Cepheid adds a component that is clearly not negligible. In the same atmospheric conditions and AO correction between the star and its calibrator, the encircled energy ratio should be close to 100%. That is what we obtain from Achernar’s data. Scaling RS Pup’s images to the same window size (i.e. r ≲ 0.42″) slightly changes E_enc because we normalized to the total flux, but we still have a discrepancy of 100% for the encircled energy ratio (ε ~ 79% and ε ~ 87% respectively in NB_1.64 and IB_2.18). A smaller window on RS Pup makes the total flux lower because we do not integrate all the envelope flux, consequently the normalized encircled energy is larger and therefore the E_enc ratio as well.

The drop in E_enc of RS Pup compared with the statistics over Achernar’s data and the fact that this drop is larger than the one directly estimated from the seeing halo variations, motivate the following section where this drop is analyzed under the assumption that it is caused by CSE emission.

3.4. Encircled energy variation

In this part we evaluate the flux of the envelope using the fact that the theoretical encircled energy ratio should be 100% and that an additional component to the flux should decrease this ratio.

The encircled energy of RS Pup can be written as ${\mathit{E}}_{\mathrm{enc}}\mathrm{=}\frac{{\mathit{F}}_{\mathrm{coh}}\mathrm{+}\mathit{\alpha }{\mathit{F}}_{\mathrm{env}}}{\mathit{F}\mathrm{+}{\mathit{F}}_{\mathrm{env}}}\mathit{,}$where F is the total flux of the star, without the extended flux, F_coh is the flux in the coherent core, F_env is the emission of the environment and α is the fraction of the envelope flux that lays in the core. For simplicity, we will assume α ≪ 1, which means that we assume that the CSE is significantly larger than the core. Accordingly we can write ${\mathit{E}}_{\mathrm{enc}}\mathrm{\approx }\frac{{\mathit{F}}_{\mathrm{coh}}\mathit{/}\mathit{F}}{\mathrm{1}\mathrm{+}{\mathit{F}}_{\mathrm{env}}\mathit{/}\mathit{F}}\mathrm{·}$Assuming that the encircled energy is stable in NACO for our observing conditions, we can calibrate the previous formula to obtain the extended emission flux: $\mathit{\epsilon }\mathrm{=}\frac{\mathit{E}\mathrm{enc}\mathrm{\left(}\mathrm{Cepheid}\mathrm{\right)}}{\mathit{E}\mathrm{enc}\mathrm{\left(}\mathrm{reference}\mathrm{\right)}}\mathrm{\approx }\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathit{F}\mathrm{\left(}\mathrm{env}\mathrm{\right)}\mathit{/}\mathit{F}\mathrm{\left(}\mathrm{Cepheid}\mathrm{\right)}}\mathrm{·}$The relative emission F(env)/F(Cepheid) can be estimated with the encircled energy ratio from our measured values (Table 2). We can also take the uncertainty values assessed in the previous section, defined as being the relative standard deviation on 10% of best frames in the cubes. We found a relative extended emission around RS Pup contributing to 38  ±  17% of the star’s flux in the NB_1.64 band and 24  ±  11% in IB_2.18.

4. Morphological analysis

We now present a qualitative analysis of the envelope morphology. In this analysis, based on the statistical properties of the noise owing to speckle boiling, we presented the ratio $\mathit{\sigma }\mathrm{\left(}\mathit{F}\mathrm{\right)}\mathit{/}\overline{\mathit{F}}$ as a function of the factor (g). This parameter is intrinsic to speckle temporal and spatial structures, is not expected to depend on azimuth, and is proportional to the envelope flux. We mapped this factor, and any azimuthal asymmetry should be attributed to the morphology of the envelope.

4.1. Analysis method

The principle of the method is to recover the flux of the envelope in the short exposures. The difference between a simple unresolved star and a star surrounded by an envelope is that the envelope will prevent the flux between the speckles from reaching zero.

Statistical descriptions of stellar speckles have been proposed by many authors (Racine et al. 1999; Fitzgerald & Graham 2006), who showed that speckle noise is not negligible. The distribution of speckle statistics in AO were also studied by Canales & Cagigal (1999). For example, Racine et al. (1999, hereafter Ra99) showed that the speckle noise dominates other noise sources within the halo of bright stars corrected with AO systems. We know that the total variance in the cube short exposures after the basic corrections (i.e. corrected for bias, bad pixels, flat-field and recentered) is ${\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{\sigma }}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}{\mathit{\sigma }}_{\mathrm{ph}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}{\mathit{\sigma }}_{\mathrm{RON}}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{Strehl}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{,}$(3)where ${\mathit{\sigma }}_{\mathrm{s}}^{\mathrm{2}}$ is the variance of the speckle noise, ${\mathit{\sigma }}_{\mathrm{ph}}^{\mathrm{2}}$ the variance of the photon noise, ${\mathit{\sigma }}_{\mathrm{RON}}^{\mathrm{2}}$ the variance of the readout noise and ${\mathit{\sigma }}_{\mathrm{Strehl}}^{\mathrm{2}}$ the noise owing to Strehl ratio variations. We neglected the variance of the sky photon noise, which is negligible in our wavelength for bright targets. Because the Strehl ratio variations occur in the coherent core and because we are particularly interested in the halo part, we omit these variations in the rest of this section.

We present for one calibrating cube the variance observed in the top panels of Fig. 5 for both filters. We plotted the radial profile computed using median rings on all azimuths. The dashed curve denotes the variance of the readout noise. We can estimate some of the components in Eq. (3) from our data. For example, from the stacked image that we obtained by summing the cube, we obtain the core and halo profile, as we did in Sect. 3. We then obtain the photon noise variance using its Poissonian property, knowing that the variance is equal to the mean. We define for the stacked image the radial intensity profile I(r) and we estimate the term ${\mathit{\sigma }}_{\mathrm{ph}}^{\mathrm{2}}$ of Eq. (3) with ${\mathit{\sigma }}_{\mathrm{ph}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathit{\gamma }\mathit{I}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathit{\gamma }\frac{\mathit{I}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{N}}\mathit{,}$where N is the number of frames in the cube (1999 frames in IB_2.18 and 457 in NB_1.64) and γ the gain of the detector in order to convert noises in electron unit. The readout noise is estimated from the outer part of the images (see in Fig. 5 the top panels). We estimate the speckle noise from the total variance of the cubes as ${\mathit{\sigma }}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{-}{\mathit{\sigma }}_{\mathrm{ph}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{-}{\mathit{\sigma }}_{\mathrm{RON}}^{\mathrm{2}}\mathit{,}$

Fig. 5 Left panel: at the top we plot the variance (circles) and the average value (triangle) of one data cube of HD 74417 at 2.18 μm. Bottom: contribution of the different noise sources. Right panel: the same for the 1.64 μm filter. The scale for all plots is logarithmic in ordinates. The dotted curves show the 3σ limits of the readout noise, while the gray curves denote the halo fit.

Fig. 6 Mean g function for each star. The error bars correspond to the RMS. The bottom plot represents the absolute relative difference. We can see a larger difference in the core owing to the Strehl variations. From r > r1 the functions are similar at a 9% level in NB_1.64 and at 6% in IB_2.18.

In the bottom panels of Fig. 5 we present the relative importance of the different noise sources. We see that the different variances show a different radial behavior, and we see in particular that the photon noise of the core never dominates. More precisely, we can define two regimes:

• σ(r) ≃ σ_s(r) for r ≲ r₁, in this regime, the speckle noise dominates the other variances (σ_s > 3   σ_RON)

• $\mathit{\sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{\simeq }\sqrt{{\mathit{\sigma }}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}{\mathit{\sigma }}_{\mathrm{RON}}^{\mathrm{2}}}$ for r ≳ r₁

with r₁ ≃ 0.22″ in NB_1.64 band and r₁ ≃ 0.24″ in IB_2.18. The fairly weak (∝λ ^− 12/5) wavelength dependency of the speckle noise may explain the same order of magnitude between the limit radius in IB_2.18 and NB_1.64.

Indeed there is another regime in the range where the Strehl variations occur. The PSF core is mainly affected by these variations. As said previously, we are interested in the PSF halo to detect the envelope between the speckles, which is why we neglect this term.

4.2. Speckle noise invariant

We use an interesting property of the speckle noise, namely that its variance is proportional to the squared flux in the cube (Ra99), F/σ is therefore an invariant. More precisely: ${\mathit{\sigma }}_{\mathrm{s}}\mathrm{=}\mathit{g}\mathrm{\left(}\mathit{r,\lambda ,}\mathrm{\Delta }\mathit{t,S,}\mathit{...}\mathrm{\right)}\hspace{0.17em}{\mathit{F}}_{\mathit{\star }}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathit{,}$where the function g is a function depending on the wavelength λ, the integration time Δt, the Strehl ratio S and atmospheric parameters.

We checked this property by plotting in Fig. 6 the mean g function for each star, i.e. we evaluated a g function for each cube that we then averaged (the error bars correspond to the standard deviation between the cubes). We took care to remove the flux contribution from the circumstellar envelope with the values found in Sect. 3. Obviously there is no clear difference between the two stars except in the core. This can be explained by the Strehl ratio variations we omitted before in the core. We can assess these variations by evaluating the encircled energy for each frame of each calibrating cube (here we take 100% of frames, without selection) and by using it as an estimator of the Strehl ratio in function of time. By fitting a Gaussian distribution to E_enc for each star and both filters, we estimate the relative standard deviation for RS Pup to be $\mathit{\sigma }\mathrm{\left(}{\mathit{E}}_{\mathrm{enc}}\mathrm{\right)}\mathit{/}\overline{{\mathit{E}}_{\mathrm{enc}}}\mathrm{~}\mathrm{15}\mathrm{\%}$ in IB_2.18 and  ~17% in NB_1.64. For HD 74417 we found  ~8% and  ~9% respectively at 2.18 μm and 1.64 μm.

The total invariance of the g function, i.e. for the complete radius, depends on the Strehl ratio variations. Nevertheless, because these variations occur in the coherent core, we can neglect them if we go far away from the central part. For r > r₁ we found a mean absolute relative difference (see Fig. 6) of 5% in IB_2.18 and 8% in NB_1.64. So for r > r₁ the invariance of the g function is verified with a good accuracy.

We can measure g using our reference star. Assuming it is stable in time, we can define the observable Γ as the average of the flux of a given pixel in the cube, divided by the standard deviation. Defining F_cal = αF_RS   Pup = αF _⋆  and making the difference between Γ_sci and Γ_cal (where the index sci and cal refer to RS Pup and HD 74417 respectively) we have $\mathit{\xi }\mathrm{=}{\mathrm{\Gamma }}_{\mathrm{sci}}\mathrm{-}{\mathrm{\Gamma }}_{\mathrm{cal}}\mathrm{=}\frac{{\mathit{F}}_{\mathit{\star }}\mathrm{+}{\mathit{F}}_{\mathrm{env}}}{\sqrt{{\mathit{g}}_{\mathit{\star }}^{\mathrm{2}}{\mathit{F}}_{\mathit{\star }}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{RON}}^{\mathrm{2}}}}\mathrm{-}\frac{\mathit{\alpha }{\mathit{F}}_{\mathit{\star }}}{\sqrt{{\mathit{g}}_{\mathit{\star }}^{\mathrm{2}}{\mathit{\alpha }}^{\mathrm{2}}{\mathit{F}}_{\mathit{\star }}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{RON}}^{\mathrm{2}}}}\mathit{,}$(4)All variables in Eq. (4) are functions of the pixel coordinates (x,y), but we omit this notation for clarity.

The parameter α can be estimated with the results from Sect. 3. The flux ratio between the two final images of RS Pup and HD 74417 is a function of α and F(env)/F _⋆ : $\frac{\mathit{F}\mathrm{\left(}\mathrm{RS} \mathrm{Pup}\mathrm{\right)}}{\mathit{F}\mathrm{\left(}\mathrm{HD} \mathrm{74417}\mathrm{\right)}}\mathrm{=}\frac{{\mathit{F}}_{\mathit{\star }}\mathrm{+}{\mathit{F}}_{\mathrm{env}}}{\mathit{\alpha }{\mathit{F}}_{\mathit{\star }}}\mathrm{=}\frac{\mathrm{1}\mathrm{+}{\mathit{F}}_{\mathrm{env}}\mathit{/}{\mathit{F}}_{\mathit{\star }}}{\mathit{\alpha }}\mathrm{·}$(5)By measuring the total flux in our final images (Sect. 3) and using F_env/F _⋆  estimated in Sect. 3.4 we get α = 1.06    ±    0.13 at 1.64   μm and α = 0.96    ±    0.08 at 2.18   μm. We also used the total flux variations in our 10 (for RS Pup) and 4 (for the reference) averaged images to assess the uncertainties, but they turn out to be small (<2%) compared to F_env/F _⋆  uncertainties. As a first approximation we can simplify Eq. (4) with α ~ 1 to obtain $\mathit{\xi }\mathrm{=}\frac{{\mathit{F}}_{\mathrm{env}}}{\sqrt{{\mathit{g}}_{\mathit{\star }}^{\mathrm{2}}{\mathit{F}}_{\mathit{\star }}^{\mathrm{2}}\mathrm{+}{\mathit{\sigma }}_{\mathrm{RON}}^{\mathrm{2}}}}\mathrm{·}$(6)ξ is thus directly proportional to the flux of the CSE. We expect this function to be null if there is no envelope. The denominator is a radial parameter, therefore any spatial asymmetry in ξ will be linked to an asymmetry of the CSE.

From our different data cubes we computed a mean ξ by subtracting to each Γ_sci a mean Γ_cal, i.e.: $\overline{\mathit{\xi }\mathrm{\left(}\mathit{x,y}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{1}}{\mathit{n}}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathit{n}}\mathrm{(}{\mathrm{\Gamma }}_{\mathrm{sci}\mathit{,}\mathrm{j}}\mathrm{-}\overline{{\mathrm{\Gamma }}_{\mathrm{cal}}}\mathrm{)}\mathrm{and}\overline{{\mathrm{\Gamma }}_{\mathrm{cal}}}\mathrm{=}\frac{\mathrm{1}}{\mathit{m}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{m}}{\mathrm{\Gamma }}_{\mathrm{cal}\mathit{,}\mathrm{i}}\mathit{,}$where n and m are the number of cube for RS Pup and HD 74417 respectively. This process provides us with an average image proportional to the flux of the CSE in each filter.

4.3. Symmetry of the CSE

In Fig. 7 we show the ξ function for the NB_1.64 and IB_2.18 band. For radii larger than r₁, we expect the function to be null if there is no CSE (see Eq. (6)), whereas we clearly obtain a positive image, indicative of a likely detection of a CSE. The central part of the image with r < r₁ (white circle for r = r₁) denotes the part where the previously defined function g is not invariant because of Strehl ratio fluctuations. It is difficult to estimate and therefore does not supply reliable informations close to the star.

This method is useful for studying the morphology of the CSE, in particular whether or not it has a central symmetry, like one expects for a uniform shell structure or a face-on disk; or an asymmetrical shape, like in an inclined disk-like structure for instance. Despite adaptive optics artefacts, the images in Fig. 7 seem to present a uniform intensity distribution in agreement with a shell or a face-on disk structure. More precisely, we estimated a symmetry level using a residual map obtained by subtracting 90 degree rotated versions of themselves from the images presented in Fig. 7. The result of this subtraction is shown in Fig. 8. In the initial image we assessed the average value in a ring, while in the subtracted image (0° − 90°) we compute an average value in a small window (see Fig. 8). Making the ratio of these two averages (i.e. the average value in the small window over the average value in the ring), we exclude a departure from symmetry larger than 7% in NB_1.64 and 6% in IB_2.18.

Fig. 7 ξ parameter (Eq. (6)) proportional to the envelope flux at 2.18 μm (left) and at 1.64 μm (right). The white circle represents the r1 limit of Eq. (7) separating the two regimes. The scale intensity is logarithmic.

Fig. 8 Estimate of the symmetry level. Top panels: on the left is plotted the ξ parameter at 2.18   μm with a ring of radius r1 < r < r2 = 0.92″ where we computed a mean value. On the right we plot the difference of the initial (0°) and 90° rotation images where we cut a small window to assess a mean residual value. The ratio between the mean residual value and the mean value give us a level of symmetry. Bottom panels: same process at 1.64   μm with a ring of radius r1 < r < r2 = 1.1″. The scale intensity is logarithmic.

5. Conclusion

We present a study of the close environment (1000–2000 AU) of the Cepheid RS Pup using near-infrared AO-assisted imaging. We used two different techniques to 1) measure the photometric emission due to the CSE and 2) obtain an image proportional to the CSE through a statistical analysis of the large number of short exposures in which speckles are still visible.

Using a simple image modeling method, based on the unresolved PSF calibrator, we estimated (Sect. 3) a likely photometric emission of 38  ±  17% and 24    ±  11% of the star flux in the NB_1.64 and IB_2.18 bands respectively for RS Pup’s CSE. Apart from the actual interest to better understand the mass loss in Cepheids, the presence of such an extended emission may have an impact on the distance estimate with the Baade-Wesselink method.

Using an original statistical study of our short exposure data cubes, we qualitatively showed (Sect. 4) that the envelope in the two bands (NB_1.64 and IB_2.18) has a centro-symmetrical shape (excluding a departure from symmetry larger than 7%), suggesting that the CSE is either a face-on disk or a uniform shell.

Our observations were obtained in two filters isolating the Brackett 12–4 (1.641   μm, NB_1.64 filter) and Brackett γ 7–4 (2.166   μm, IB_2.18 filter) recombination lines of hydrogen. We therefore propose that our observations show the presence of hydrogen in the CSE of RS Pup. The star is naturally relatively faint in these spectral regions, because of deep hydrogen absorption lines in its spectrum. This appears as a natural explanation to the high relative contribution from the CSE compared to the stellar flux. Another contribution to the CSE flux may also come from free-free continuum emission, or Rayleigh scattering by circumstellar dust. However, this contribution is probably limited to a few percents of the stellar flux. The high relative brightness of the detected CSE therefore implies that the spectral energy distribution of the detected CSE is essentially made of emission in hydrogen spectral lines instead of in the form of a continuum. RS Pup’s infrared excess in the broad NB_1.64 and IB_2.18 bands is limited to a few percents (Kervella et al. 2009a), which excludes a continuum CSE emission of 40% in the NB_1.64 filter and 24% in the IB_2.18 filter as measured in our NACO images.

Our discovery is also remarkably consistent with the non-pulsating hydrogen absorption component detected in the visible Hα line of RS Pup by Nardetto et al. (2008), together with a phase-dependent component in emission. The extended emission we observe in the infrared is therefore probably caused by the same hydrogen envelope as the one that creates the Hα emission observed in the visible. It should also be mentioned that our result is consistent with the proposal by Kervella et al. (2009a) that the envelope of RS Pup is made of two distinct components: a compact gaseous envelope, which is probably at the origin of the observed NACO emission, and a very large ( ≈ 2′) and cold ( ≈ 45 K) dusty envelope (Kervella et al. 2008) that is mostly transparent at near-infrared wavelengths owing to the reduced scattering efficiency. These two envelopes appear to be largely unrelated because of their very different spatial scales. However, it should be noted that Marengo et al. (2009) recently detected large envelopes around several Cepheids using Spitzer. The presence of hydrogen close to the Cepheids of their sample could be investigated at high angular resolution using the observing technique presented here.

Acknowledgments

We thank Drs. Jean-Baptiste Le Bouquin and Guillaume Montagnier for helpful discussions. We also thank the ESO Paranal staff for their work in visitor mode at UT4. We received the support of PHASE, the high angular resolution partnership between ONERA, Observatoire de Paris, CNRS, and University Denis Diderot Paris 7. This work made use of the SIMBAD and VIZIER astrophysical database from CDS, Strasbourg, France and the bibliographic informations from the NASA Astrophysics Data System. Data processing for this work have been done using the Yorick language, which is freely available at http://yorick.sourceforge.net/.

References

All Tables

All Figures

	Fig. 1 Final frames from shift-and-add process. In the two upper panels the PSF calibrator star HD 74417 and RS Pup are shown in the IB_2.18 filter. The bottom panels show the same stars in the NB_1.64 filter. The scale intensity is logarithmic.
In the text

	Fig. 2 On the left are presented the radial profiles at 2.18 μm for RS Pup, Achernar and their reference stars. On the right the same at 1.64 μm. Radial profiles were computed using the median value on one-pixel thick circular rings. The gray color denotes Achernar and δ Phe profiles, shifted by a factor 10 for clarity. The scale is logarithmic for both axes.
In the text

	Fig. 3 Histogram of the encircled energy computed for 10% best frames. The solid lines represent a Gaussian fit.
In the text

Fig. 4 Halo fit for all filters with logarithmic ordinates. The fluxes are normalized to the total flux. Left panel: halo fit from the radial profile of images in the two filters for RS Pup. Right panel: halo fit from the radial profile of images in the two filters for Achernar. The fit in IB_2.18 was for r > 0.26″ and in NB_1.64 for r > 0.2″ for RS Pup and for r > 0.26″ for Achernar. The dashed lines denote the reference stars, while the thick lines denote the Cepheid and Achernar. The fitted curves are plotted in gray.

In the text

	Fig. 5 Left panel: at the top we plot the variance (circles) and the average value (triangle) of one data cube of HD 74417 at 2.18 μm. Bottom: contribution of the different noise sources. Right panel: the same for the 1.64 μm filter. The scale for all plots is logarithmic in ordinates. The dotted curves show the 3σ limits of the readout noise, while the gray curves denote the halo fit.
In the text

	Fig. 6 Mean g function for each star. The error bars correspond to the RMS. The bottom plot represents the absolute relative difference. We can see a larger difference in the core owing to the Strehl variations. From r > r1 the functions are similar at a 9% level in NB_1.64 and at 6% in IB_2.18.
In the text

	Fig. 7 ξ parameter (Eq. (6)) proportional to the envelope flux at 2.18 μm (left) and at 1.64 μm (right). The white circle represents the r1 limit of Eq. (7) separating the two regimes. The scale intensity is logarithmic.
In the text

Fig. 8 Estimate of the symmetry level. Top panels: on the left is plotted the ξ parameter at 2.18   μm with a ring of radius r1 < r < r2 = 0.92″ where we computed a mean value. On the right we plot the difference of the initial (0°) and 90° rotation images where we cut a small window to assess a mean residual value. The ratio between the mean residual value and the mean value give us a level of symmetry. Bottom panels: same process at 1.64   μm with a ring of radius r1 < r < r2 = 1.1″. The scale intensity is logarithmic.

In the text