© ESO, 2011

1. Introduction

The supergiant A[e] star (Lamers et al. 1998) HD 62623 (HR 2996, 3 Puppis, l Puppis), is a key object for understanding the processes at the origin of aspherical shells in massive evolved stars (Humphrey & Davidson 1994) and supernovae (Kirshner 1987). Indeed, HD 62623 is surrounded by a dense gaseous and dusty disc (Meilland et al. 2010), a structure more often found in young stellar objects and post-AGB stars (van Winckel 2006). Discs are known to govern accretion or mass-loss in these lower-mass objects, but their origin and structure remain highly debated for massive stars (Porter 2003).

Fast rotation of the star leads to an expanding disc-like wind in the case of very massive stars (Zickgraff 1985; Lamers & Pauldrach 1991) or, when viscosity becomes dominant in less massive stars, it leads to a rotating disc (Lee et al. 1991). The presence of a companion star could also lead to a rotating disc (Plets & Trams 1995).

Fast rotation, or the presence of a companion star, could be responsible for the breakup of spherical symmetry of mass-loss around the massive, hot and luminous object at the core of HD 62623. To test these two hypotheses, one needs to access unprecedented combined spatial and spectral resolutions. Here, we report the first continuum image and velocity-resolved images in a circumstellar disc using a new multi-wavelength optical interferometry imaging method, which allows us to spatially disentangle the dust and gas emissions of HD 62623.

2. VLTI/AMBER observations, and data reduction

To unambiguously resolve the close environment of HD 62623, both spatially and spectrally, we acquired data in early 2010 (Table 1) using the VLTI/AMBER instrument (Petrov et al. 2007). Each AMBER measurement consists of three visibilities, one closure phase, three wavelength-differential phases, and one flux measurement, each of them spectrally dispersed on about 500 narrow-band channels, with a spectral resolution power of 12 000 close to the Brγ line. We used the standard AMBER package (Tatulli 2007) to reduce the data, complemented with calibration scripts (Millour et al. 2008). The dataset (≈54 000 visibilities,  ≈54 000 differential phases, and  ≈18 000 closure phases, presented Fig. 2) and resulting UV coverage, presented in Fig. 1, are noticeably the largest obtained so far with the VLTI.

Fig. 1 The UV coverage of our observations. Different colours represent different nights. The ellipse has dimensions 82 × 128 m.

Our HD 62623 squared visibilities (top-middle of Fig. 2) decrease with increasing spatial frequencies with values close to 0 at high frequencies, meaning that we resolved the close environment of HD 62623. In addition, we detect non-zero closure phases, indicating asymmetries in the environment of HD 62623 (top-right panel of Fig. 2). In addition, we clearly resolve changes in the visibilities, closure phases and differential phases in the Brγ line (bottom panels of Fig. 2).

Fig. 2 Top-Left: the spectrum. The double-peaked emission line is the Brγ line. Top-Middle: the squared visibilities as a function of spatial frequency, Top-Right: the closure phases as a function of spatial frequency, Bottom-Left: the differential visibilities (visibilities normalized to 1 in the continuum), Bottom-Middle: the closure phase, and Bottom-Right: the differential phase, the three former as a function of wavelength. In all plots, different colours represent different nights, and the black lines (resp. the black dots for the top panels) represent our best-fit model. “cycles/as”. or “cyc./as”. means “cycles per arc-seconds” or “arc-seconds-1”. For the closure phase plot, the spatial frequencies are the ones of the longest baseline. In the squared visibilities panel, one can see that the disc is clearly resolved at the largest baselines.

We then used the MIRA software (Thiébaut & Foy 2003) to recover a spectrally dispersed image cube of HD 62623, using in a first attempt squared visibilities and closure phases. We did not try other image reconstruction software, based on the conclusion of Millour et al. (2009). In addition, we developed a dedicated self-calibration method (presented in Sect. 3) to add the wavelength-differential phases to the image recovery process, MIRA being used this time with squared visibilities, closure phases and phases. This additional step contribute to link the astrometry of all narrow-band images together, allowing us to measure astrometric changes in the emission line compared to the continuum.

To calibrate the phases sign (and, hence, the on-sky orientation of the images), we also performed reconstructions using MIRA on the θ  Ori C dataset from Kraus et al. (2009)¹.

3. The self-calibration method

3.1. The wavelength-differential phase

Wavelength-differential phase in optical/infrared stellar interferometry is a partial measure of the object Fourier phase φ_object(λ). Indeed, the phase measured on a ground-based long-baseline interferometer is affected by randomly variable perturbations. These perturbations are composed of, by decreasing magnitude (in the near-infrared):

• δ, a variable achromatic optical path difference (OPD) between the two telescopes;

• δ_dry, a chromatic dry air OPD;

• δ_wet, a chromatic wet air OPD.

The measured phase ϕ_measured takes therefore the following form: (1)The differential phase computation algorithm (see Millour et al. 2006, for the application to AMBER) aims at removing the term 2π [δ(t) + δ_dry(t,λ) + δ_wet(t,λ) ] /λ from this equation. This term can be described by the following first-order Taylor-expansion: (2)A very simple model-fitting to the observed real-time phase data yield α and β time sequences. A subsequent subtraction and averaging is then performed to compute the wavelength-differential phase. In this fitting procedure, any phase term from the object that could mimic an atmospheric OPD will be removed. Therefore, at first order, differential phase is the object’s phase where both an offset and a slope have been subtracted: (3)In practice, removing this offset and slope from the object phase makes the wavelength-differential phase free of any absolute astrometric information from the object. On the other hand, the wavelength-differential phase does give relative astrometry information at one wavelength compared to the one at another wavelength, giving access to spectro-astrometric measurements (see Millour et al. 2007, for a detailed example of application).

3.2. Self-calibration applied to optical/infrared interferometry

Compared to the radio-astronomy phase, this optical/infrared wavelength-differential phase bears many similarities. The biggest difference is that its average is set to zero and that its average slope is also set to zero. A previous work (Millour 2006, pp. 63–69, especially Fig. 3.7) demonstrated, in theory, the potential of wavelength-differential phase for use in interferometric image synthesis. Inspired by a radio-astronomy review article by Pearson & Radhead (1984), we realise here the practical application of this development, by implementing an equivalent to the Hybrid Mapping algorithm (also called self-calibration algorithm), using the MIRA software instead of the CLEAN algorithm to produce the hybrid maps, and the wavelength-differential phases instead of radio-phases. The main differences in our case compared to radio-interferometry are:

• the measured phases are adjusted with a linear-fit to the trial onesinstead of just offset. This is due to the way thewavelength-differential phase is obtained.

• the first image-reconstruction step is made using only squared visibilities and closure phases, and each narrow-band image is centred to its brightest point.

A sketch explaining this method is shown in Fig. 3, and a summary of the procedure is described below:

• 1.

  we start the image reconstruction process for each wavelength from a set of squared visibilities and closure phases;

• 2.

  the images are centered on their brightest point to limit the astrometric changes between these independent images;

• 3.

  then, the self-calibration process itself is applied, by comparing the observed wavelength-differential phases (solid line in the middle-top panel of Fig. 3) to the phase computed from the image cube (hybrid map, dashed line in Fig. 3). The whole dataset is modified so that both the phase offset and slope match the image cube ones;

• 4.

  then, the image reconstruction is performed again, adding the synthetic phases data to the squared visibilities and closure phases;

• 5.

  the result of this image reconstruction step can then be input again into the self-calibration process (step 3) or used as the final image cube.

In our case, we used two of these iterative steps to reach the final result. Using more iteration does not improve image quality, and at some point introduces images artifacts, like would do any other image-recovery process.

Fig. 3 The principle of the self-calibration algorithm applied to optical/IR interferometry. Solid arrows show processes that may be repeated several times, while dashed arrows show processes carried out once.

In complement to the described self-calibration algorithm, we used the following parameters in MIRA for the image reconstruction: a quadratic regularization, and a Gaussian prior (to limit the effective reconstruction field of view to approximately 20–30 mas). These parameters set less constraints on the image reconstruction process than in other works, where a severe limitation of the reconstruction field of view (Monnier et al. 2007; Zhao et al. 2009) or assumptions on the circular-symmetry of the object (Le Bouquin et al. 2009; Chiavassa et al. 2010) were applied to reconstruct the images.

In the last processing steps, we are averaging all the continuum images to produce the continuum map and subtract the continuum image from each spectral channel images to produce the emission-line images. We tried the overall process using random initial images to check the unicity of our solution, and checked also on simulated datasets the consistency of the whole image reconstruction process.

4. Results

Figure 4 shows the obtained model-free image where we clearly see a ring of material of radius  ≈6 mas (4 AU at 650 pc) encircling a bright inner spot of radius  ≈2 mas (1.3 AU at 650 pc). The ring corresponds to the region of the previously-detected dusty disc by Meilland et al. (2010) where dust is heated at the sublimation temperature (inner rim), and the bright spot is a joint contribution from the star, and free-free and free-bound emission from an inner gaseous disc.

Fig. 4 Top-Left: composite aperture-synthesis image of HD 62623. The convolving beam of 1.75 × 2.76 mas is shown in the lower-left box at the same scale. The dashed line represent the polarization measurement from Yudin & Evans (1998). Top-Right: our best-fit model plot at the same scale and convolved with the same beam. Bottom-left: line images of HD 62623. The line intensity is plot as colours for different gas velocities. The intensity in the continuum is plot as gray contours (80%, 40%, 20%, 10%, 5%, 2.5% and 1.25% of the peak intensity). Bottom-Right: our best-fit model line images with the same representation and the same convolving beam.

The continuum image does not show a clear elongation of the ring because of the elongated beam in the N-S direction, but clearly shows that the ring is brighter on the West side, compared to the East side. The inclination angle of the dust disc must therefore be intermediate (30 to 60°) and the position angle of the major axis of the disc must be close to zero degrees, as we see a brighter emission zone opposite to the star and a dimmer, hidden emission zone closer to us. This is compatible with the values predicted from mid-IR observations (position angle of 10 ± 10°, inclination of 60 ± 10°) and also with the measured polarization angle by Yudin & Evans (1998): 95 ± 5°, perpendicular to the major axis of the disc.

Confirmation that the inner spot is partly composed of rotating ionized hydrogen (observed here in Brγ) comes from the line emission images, where we see the red-shifted emission appearing on the North side of the star, and the blue-shifted emission on the South side. At the line center, the emission is roughly symmetric and the elongation is directed to the East–West. This signal is a typical signature for a rotating disc. The Brγ emission in the disc vanishes rapidly with increasing radius. We argue that this is not only a consequence of the decreasing gas density, but also of the low (<10 000 K) star temperature, and, thus, of the larger neutral hydrogen fraction with increasing radius. No Brγ emission comes from the ring, which means the gaseous hydrogen mixed with the dust is not ionized in this region.

The central star is unresolved by the VLTI in the continuum image, but we resolve the central gaseous disc. The lack of line emission in the very central region (see Fig. 6 and Sect. 6) implies that the VLTI resolved an apparent “cavity” (radius not exceeding 1 mas, i.e. 0.65 AU) between the core regions including the A3Ia photospheric stellar source (producing an absorption component) and the emitting disc of plasma. This confirms that the star’s wind do not contribute significantly to the Brγ emission, and provides an upper limit for its radius in accordance to its spectral type, i.e. 65–80 R_⊙ (for d = 650 pc, 0.3 AU).

5. Modelling

To match our dataset in the continuum, we used a simple geometrical model composed of four components: the star (described by an uniform disc with an absorption spectrum), the free-free and free-bound emission from the gas (described by an elliptical Gaussian disc), the inner rim of the dusty disc (described by an elliptical skewed ring) and the outer regions of the disc (described here by an extended background).

In addition, to match the data in the emission-line region, we add an elliptical Gaussian-brightness component which is assigned the following combination of velocity fields (see Stee et al. 1995; Meilland et al. 2007, for typical applications of these velocity fields):

• a radial expansion velocity field of the form: with v_r(0) the expansion speed at the smallest radius, v_r(∞) the expansion speed at infinite radius, and γ the exponent of the power law;

• a rotation velocity field of the form: with v_θ(0) the maximum rotation speed of the disc. When β equals –0.5, the disc is in Keplerian rotation, when β equals 1 the disc is in solid rotation, and when β equals –1, the disc has a constant angular momentum.

Parameters of the whole model range from angular sizes of the components (star, gas disc and dust disc), the system inclination angle, flux ratios between the star, the gas disc, and the dust disc, and also the rotation and expansion velocity, the exponents of the rotation and expansion law, in addition to the on-sky orientation angle.

This model produces intensity maps that are Fourier-transformed, and used to produce synthetic observables (squared visibilities, closure phases and differential phases). These synthetic observables are adjusted to the observed data using a combination of a gradient descent algorithm and a simulated annealing algorithm. Given the large amount of observed data, we fitted this model having all parameters potentially free.

In addition to our global fit, we tested several typical models (i.e. fixing few parameters and leaving the other ones free in the fitting procedure described above), including: a model with free rotation and/or expansion laws, a pure Keplerian rotation law, and a solid rotation.

Our model-fitting shows that there is no or very little expansion in the disc around HD 62623, and that rotation is indeed at the core of the disc velocity law. To check that our best obtained model is located in an absolute minimum, we produced a χ² curve, starting from our best solution, fixing the β parameter to several values, and fitting the model to the data by using the method described above.

Our final reduced χ² is equal to 1.92 for a Keplerian-rotating disc model (See Fig. 5). We find a slightly lower χ² for β = −0.7, i.e. a rotation law between Keplerian rotation and Angular momentum conservation (see Table 2). However, all χ² values, except for β ≈ 0, lie between 1.9 and 2.2. This basically means that the significance of this result is quite low, probably due to the little number of angular-resolution elements we have in the gas disc.

We conclude that, though the lowest χ² values lie close to Keplerian rotation (i.e., rotation law exponent –0.5), our observation only weakly contrain the rotation law of the gas disc. On the other hand, we can clearly rule-out expansion-dominated velocity fields (χ² > 4 for all models), confirming that rotation is indeed the main velocity field in the HD 62623 disc.

Fig. 5 χ2 curve obtained by fixing the β parameter to different values, and fitting the model to the data. Dashed vertical lines indicate several canonical rotation laws, and arrows represent χ2 values above 2.2.

6. Discussion

The spectrally and spatially resolved AMBER/VLTI observations in the continuum and Brγ line of HD62623 bring unprecedented information on the disc surrounding an A supergiant.

The central star appears unresolved by the VLTI in the continuum image, providing an upper limit for its radius in accordance for its spectral type, i.e. 65–80 R_⊙ (resp. 130–160) assuming d = 650 pc (resp. 1300 pc) solar radii. The decrease of the line emission in the very central region, giving a doughnut-like shape for the integrated line flux (see Fig. 6), implies that the VLTI resolved an apparent cavity in the innermost region of the emitting disc of plasma (radius not exceeding 1 mas, i.e. 0.65 AU), which includes the A3Ia photospheric stellar source. This implies that there is not enough emission in the vicinity of the central star to fill-in this lack of flux. This is a new and clear evidence that the star’s wind do not contribute significantly to the Brγ emission, and that the spectral type of the central source is not at odd with a normal member of its class. The disc of plasma is therefore at the origin of most, if not all of the emission lines seen the spectrum of HD62623.

Fig. 6 Line-integrated flux from HD 62623, with the same scale and convolution as Fig. 4. Solid contours are at 80, 60, 40, 20, 10 and 0% of the peak intensity. Dashed contours represent negative contributions. Note the doughnut-like shape of the line-emission.

The best conditions in terms of temperature and density to generate a Brγ flux are met close to the star, and the Brγ emission disappears rapidly with increasing radius. We can argue that this is not a only a consequence of the decreasing density, but also of the low (10 000 k) star temperature and thus the relatively larger neutral hydrogen fraction with increasing radius.

How was such a dense disc formed? Even in the case the central star would be considered a fast rotator, its rotation rate is lower than the nearby rotating plasma (Vsini ranging from 35 to 80 km s^-1, depending on the work: see for instance Abt & Morrell 1995; Levato 1972; Huang 1953; Westgate 1933).

Previous work by Plets & Trams (1995) provided evidence for a companion star to HD 62623, from radial velocity variations. In our images, no companion is detected to about 1% of the peak intensity (i.e. the typical flux level of image artifacts), an upper limit compatible with the low mass ratio inferred from the Doppler signal (mass of the companion of about 1–2 M_⊙). The angular momentum transfer from an orbiting companion may explain a Keplerian velocity field. Such an energetic argument is widely used in the context of the lower mass evolved stars (van Winckel 2003). Furthermore, HD 62623 is a relatively cool source among massive stars, exhibiting a radiative wind not strong enough to easily justify the accumulation of such dense material in the orbital plane in its vicinity.

The spectroscopic observations of Plets & Trams (1995) provide an orbital period of 131d or 161d. In the light of the VLTI maps, two options can be considered:

• one can place the companion close to the supergiant, in the detected cavity, and the plasma and dusty emission originates from a circumbinary disc;

• one can place the companion further away, perhaps close to the inner rim of the dusty disc. In this case, the disc of plasma is circumprimary, whereas the dust disc can be considered as circumbinary.

On the basis of the third Kepler’s law, one can exclude a separation larger than 3 AU, since it would imply a total mass of the system larger than 100 Msun (for a period of 161 days). For the same reasons one can also exclude a separation much smaller than 1.5 AU, since it would imply a total mass smaller than few solar masses, i.e. smaller than the current mass of the primary assuming a supergiant luminosity class. As demonstrated in Plets & Trams (1995), the mass-ratio is extreme and the primary harbors most of the mass. Hence, the probable separation has to be found within the smallest separations, in the range 1.2–1.7 AU.

Furthermore, assuming a circular orbit would imply a minimum velocity of about 250–350 km s^-1, not compatible with any observational features detected to date. In particular, the low excitation [OI] forbidden line formed close to the dust forming region shows a double peak with a separation of 55 km s^-1, suggesting a continuity of the Keplerian kinematics from the dense regions probed by Hα or Brγ toward the dusty extended ones. Furthermore, significant disturbances of the circumstellar medium would be expected that are not seen in the reconstructed images (which are, however, affected by a limited dynamic range).

We favor the scenario proposed in Plets & Trams (1995), that puts the companion very close to the star, within the Brγ emission-free cavity observed in our reconstructed images. An orbital radius of 1.2–1.7 AU corresponding to a total mass of the system of about 15–20 solar masses corresponds at a distance of 650 pc (resp. 1300 pc) to 1.8–2.6 mas (resp. 0.9–1.3), and an efficient transfer of angular momentum route may occur via the L2 co-orbiting Lagrange point, moving at a velocity of 130–170 km s^-1. Accounting for the new constraints on the system inclination, this value is in accordance to the Hα and Brγ velocities of 120 km s^-1 inferred from the double peaks separation in the emission line, both transitions emitting efficiently from this location. This scenario also provides a natural explanation for the anomalously large mass-loss of the A3 supergiant, whose external layers should experience strong tides under the influence of the close companion.

7. Conclusion

These observations represent a new step in our understanding of the formation of flat, dense, Keplerian discs around evolved stars. These excretion discs share many physical characteristics with the accretion discs encountered around young stellar sources, and are now also unambiguously detected around binary post-AGBs, i.e. systems in which a low to intermediate mass star expels its remaining envelope, in the process of becoming a hot white dwarf (van Winckel 2003; van Winckel 2006). In these low-mass binary systems, an efficient mass-transfer has occurred when one of the star gets giant, leading to the formation of potentially long-lived dusty discs. As such HD 62623 represents a missing link to the massive counterparts of these mass-losing stars, showing that a structured disc can also exist around such a highly luminous evolved star, and that this disc originates from the presence of a companion star.

A different mechanism might lead to the formation of flattened dusty structures around other massive stars, but this finding definitely opens a promising route of interpretation of highly asymmetrical environments around massive stars in the frequent cases in which a low mass companion is suspected, albeit hardly detectable.

Detailed physical modeling of the AMBER and MIDI data and a more complete study of the self-calibration algorithm applied to optical/IR interferometry will be presented in forthcoming papers.

The presented self-calibration imaging algorithm practically demonstrates the potential of using the wavelength-differential phase in an image-reconstruction process. This self-calibration method is easily applicable to optical/infrared interferometric data and improves greatly the quality of the reconstructed images (which are now astrometrically linked from wavelength-to-wavelength).

In addition, this imaging method opens a wealth of possibilities that were considered as unfeasible before this paper: interferometric “integral-field spectroscopy”, potentially reliable images from 2-Telescopes datasets, etc. This clearly opens new fields of research both in the signal processing domain (image synthesis techniques) and in the astrophysical domain (image-synthesis on “faint” targets, where the number of combined telescopes will always be limited: YSOs and AGNs).

Acknowledgments

We thank E. Thiébaut for providing the image-reconstruction software MIRA, M. Kishimoto, L. Rolland, K. Ohnaka for fruitful discussions, A. Merand, and S. Stefl for supporting the observations. The Programme National de Physique Stellaire (PNPS), the Institut National en Sciences de l’Univers (INSU), and the Max Planck Institute for Radioastronomy (MPIfR) are acknowledged for their financial and observing time support. This research has made use of services from the Centre de Données de Strasbourg (CDS), from the Jean-Marie Mariotti Centre (JMMC), and from the NASA Exoplanet Science Institute (NExScI) to prepare and interpret the observations.

References

All Tables

All Figures

	Fig. 1 The UV coverage of our observations. Different colours represent different nights. The ellipse has dimensions 82 × 128 m.
In the text

Fig. 2 Top-Left: the spectrum. The double-peaked emission line is the Brγ line. Top-Middle: the squared visibilities as a function of spatial frequency, Top-Right: the closure phases as a function of spatial frequency, Bottom-Left: the differential visibilities (visibilities normalized to 1 in the continuum), Bottom-Middle: the closure phase, and Bottom-Right: the differential phase, the three former as a function of wavelength. In all plots, different colours represent different nights, and the black lines (resp. the black dots for the top panels) represent our best-fit model. “cycles/as”. or “cyc./as”. means “cycles per arc-seconds” or “arc-seconds-1”. For the closure phase plot, the spatial frequencies are the ones of the longest baseline. In the squared visibilities panel, one can see that the disc is clearly resolved at the largest baselines.

In the text

	Fig. 3 The principle of the self-calibration algorithm applied to optical/IR interferometry. Solid arrows show processes that may be repeated several times, while dashed arrows show processes carried out once.
In the text

Fig. 4 Top-Left: composite aperture-synthesis image of HD 62623. The convolving beam of 1.75 × 2.76 mas is shown in the lower-left box at the same scale. The dashed line represent the polarization measurement from Yudin & Evans (1998). Top-Right: our best-fit model plot at the same scale and convolved with the same beam. Bottom-left: line images of HD 62623. The line intensity is plot as colours for different gas velocities. The intensity in the continuum is plot as gray contours (80%, 40%, 20%, 10%, 5%, 2.5% and 1.25% of the peak intensity). Bottom-Right: our best-fit model line images with the same representation and the same convolving beam.

In the text

	Fig. 5 χ2 curve obtained by fixing the β parameter to different values, and fitting the model to the data. Dashed vertical lines indicate several canonical rotation laws, and arrows represent χ2 values above 2.2.
In the text

	Fig. 6 Line-integrated flux from HD 62623, with the same scale and convolution as Fig. 4. Solid contours are at 80, 60, 40, 20, 10 and 0% of the peak intensity. Dashed contours represent negative contributions. Note the doughnut-like shape of the line-emission.
In the text