© ESO 2020

1. Introduction

Many important open issues related to the final stages of stellar evolution can only be addressed once we comprehend the geometry of circumstellar material close to the star. For massive stars (8 − 30 M_⊙; Wang & Wheeler 1998), which eventually explode as core-collapse supernovae (SN), the problem in understanding mass-loss rate and mass-loss geometry is particularly acute (Heger & Langer 1998). Mass-loss events impact angular momentum evolution and final mass, and thus fate of the individual massive star, and these events create the circumstellar environment with which SN ejecta may interact. Following a core-collapse explosion, the mass-loss history of the progenitor continues to be relevant because of its repercussion on light-curve morphology and probably also on observed aspherical SN remnants (e.g. Patat et al. 2011; Moriya et al. 2014). The processes responsible for shaping these spatial asymmetries are largely unknown.

This paper aims to provide insight into the physical structure, geometry, mass-loss episodes, and properties of the host star, which are amongst the main unknowns in theoretical models on stellar evolution. A suitable class of stars to study mass-loss events are post-red supergiants (post-RSG), but only a few such objects are known (Oudmaijer et al. 2009). The A- to K-type yellow hypergiants (YHGs), such as IRAS 17163−3907, IRC+14020, and HD 179821, that show evidence of circumstellar dust and high mass-loss rates, are great post-RSG candidates, and therefore laboratories to study the mass-loss events that take place during the post-RSG evolution.

In this study we focus on IRAS 17163−3907 (hereafter IRAS 17163), also known as Hen 3-1379 and WRAY 15-1676, which is one of the brightest infrared sources in the sky. It is embedded in a dusty circumstellar envelope, which was dubbed the Fried Egg Nebula by Lagadec et al. (2011) because of its peculiar morphology in mid-infrared (mid-IR) images taken using the VLT Imager and Spectrometer for mid-Infrared (VISIR), showing a double detached shell and timescales of ejection that reach a few hundred years.

The distance to IRAS 17163, and therefore its class, has been under debate in recent decades. Lebertre et al. (1989) classified this object as a proto-planetary nebula (post-AGB star) based on its estimated luminosity (∼10⁴ L_⊙) and observed spectral characteristics at a distance of 1 kpc. This distance was a rough estimate based on a radial velocity of 5.3 km s⁻¹ local standard of rest (LSR), which was derived using the Hα emission; the luminosity was derived based on observed rather than dereddened photometry. Subsequent studies found this object to be an early A-type star at a distance of 4 kpc; this was estimated using the radial velocities of the interstellar K I doublet absorption, yielding a bolometric luminosity of 5 × 10⁵ L_⊙, i.e. a supergiant (Lagadec et al. 2011). This finding confirmed the status of the Fried Egg Nebula as a massive star (estimated mass 25–40 M_⊙), which is most likely evolving quickly towards rather than away from the blue part of the Hertzsprung–Russell diagram (HRD; see e.g. Meynet et al. 2015), and in all likelihood with a final fate as a SN.

Lagadec et al. (2011) suggest that IRAS 17163 is one of the very few examples of a star in transition from the RSG stage to the Wolf-Rayet or luminous blue variable (LBV) phase (IRC+14020, HD 179821, ρ Cas, HR 8752). We note that Smith et al. (2004) suggested that YHGs might be missing LBVs, also known as S Doradus variables. Such objects change their underlying temperatures and wind properties on timescales of years to decades (Humphreys & Davidson 1994; Vink & de Koter 2002). The central star of the Fried Egg Nebula shows Hα profiles similar to P Cygni and slow winds similar to those of YHGs and S Dor variables. A thorough discussion on how to tell those objects apart is presented in Humphreys et al. (2017).

Wallström et al. (2015) point out the difficulty in determining the distance to the star and reported a range of 1–7 kpc depending on the method used. In the current work we used the new parallactic measurements from Gaia Data Release 2 (DR2; Gaia Collaboration 2018) and provide a new distance estimation (Sect. 4), which brings the object as close as 1.2 kpc; we adopt this new distance in the entire paper. This is three times closer than found in the most recent studies, and therefore a revised analysis of the object is necessary. In this paper we show that despite its closer distance, when the photometry is treated properly accounting for a visual extinction of approximately ten magnitudes, a previously neglected effect, IRAS 17163 maintains its luminosity and therefore its class.

IRAS 17163 is ideally located in the sky to observe with the Very Large Telescope Interferometer (VLTI) in order to study its innermost structures and provide immediate information on this enigmatic evolutionary phase. The dust shells observed in the infrared (IR) are located at large distances from the star (> 1″) and they are seemingly circular symmetric, but inhomogeneous (e.g. Herschel Space Observatory, VLT/VISIR; Hutsemékers et al. 2013; Lagadec et al. 2011). The warm (10 μm) emission traces the most recent mass-loss episodes, while the origin and spatial distribution of the 2 μm emission is being investigated in this study. In the case of the Fried Egg Nebula, the outflowing material tends to become more spherically distributed as the distance from the source grows. Recent observations with the Atacama Compact Array (ACA) of the Fried Egg Nebula reveal the presence of a red-shifted spur (∼20″), which may trace an unidirectional ejection event (Wallström et al. 2017). Investigating the present day mass-loss provides us with information on the geometry of the stellar wind. The only other object for which this property was investigated in milli-arcsecond (mas) scales (i.e. present-day wind), is the post-RSG IRC+10420 (Oudmaijer & de Wit 2013). From near-infrared (NIR) interferometric data taken with the Astronomical Multi-BEam combineR (AMBER) the presence of an ionised, bi-polar flow was inferred, but sparse uv coverage made the result model dependent.

The key evolutionary phase, its brightness and roundness, the sparsity of observed high-mass post-RSGs, and its excellent location on the sky have led us to initiate an in-depth, high spatial resolution study of IRAS 17163 and its associated Fried Egg Nebula using spectra and images from VLT and VLTI.

In Sect. 2 we describe our optical and NIR spectroscopic (X-shooter), spectropolarimetric, and interferometric (GRAVITY, AMBER) observations and in Sect. 3, we present the observational results. In Sect. 4 we present the fundamental parameters of IRAS 17163, including the revised distance, its spectral type (i.e. effective temperature), and the dereddened photometry. In Sect. 5 we spatially resolve the inner parts of the nebula, and we trace geometries of the material down to a few astronomical units from the central object. To do so we apply simple geometric modelling, which allows us to constrain the size and geometry of the 2 μm continuum emission as well as the Brγ emission and the Na I doublet emission using our VLTI/GRAVITY and AMBER observations. In addition, in the same Sect. 5 we perform and present for the first time the model-independent image reconstruction of the object at scales so close to the star (milli-arcsecond scales; few astronomical units). In Sect. 6 we apply a 2D radiative transfer model to simultaneously fit the dereddened available photometry and the radial profiles of VISIR images (8.59 μm, 11.85 μm, and 12.81 μm; Lagadec et al. 2011) for the first time. Lastly, the discussion and the summary of the paper are presented in Sects. 7 and 8, respectively.

2. Observations and data reduction

2.1. Spectroscopy and spectropolarimetry

2.1.1. X-shooter spectroscopy

The observations of IRAS 17163 ($\mathrm{RA}={17}^{\mathrm{h}}{19}^{\mathrm{m}}49\stackrel{\text{s}}{.}3$, $\mathrm{Dec}=-39°10\prime 37\stackrel{″}{.}9$ [J2000]) took place during four observing sessions in 2017 using X-shooter at ESO’s VLT. Two short visits were made in April 2017 and two longer visits took place in May 2017. X-shooter captures the spectrum of celestial objects from 0.3 to 2.4 μm in one shot, using three arms which cover the UV-blue arm (UVB 300–560 nm), visible (550–1020 nm), and NIR (1020–2480 nm) parts of the spectrum, respectively (D’Odorico et al. 2006); the technical overview of these observations is presented in Table A.2. The observations were conducted in nodding mode with a nodding step of 6″. The slit width was 1.0 arcsec in the UVB, and 0.4 arcsec in both the visible and NIR arms. The set-up results in a spectral resolution of order 10 000 in the VIS and NIR setting and slightly lower, 9000 in the UVB. The data were reduced using the standard ESO pipeline.

Given the extreme colours of IRAS 17163, that is very faint in the blue (B ∼ 17), moderately bright at V (∼13), and very bright at NIR wavelengths, K ∼ 2.4 (2MASS; Cutri et al. 2003), the observing strategy was adjusted accordingly. Observing the blue spectrum required multiple ∼20 min exposures, while the NIR spectra were taken with short exposure times going down to the minimum possible exposure times of 0.1 s. The total exposure time in the UVB setting is of order 6000 s split over four dates. The longest exposures were taken in May 2017. The VIS setting was observed on all days as well, but the main data presented on this work are multiple 60 s exposures obtained in May 2017. Most NIR data are saturated, even those with the shortest exposure times. After visually inspecting all (hundreds of) exposures, we identified a handful of spectra that were not affected. These are mostly from the May 29 run.

Considering that this is one of the brightest infrared sources in the sky, it is interesting to note that the data discussed in this work are not only the first high signal-to-noise (S/N) data in the blue part of the spectrum, but surprisingly, they also constitute the first high-resolution, JHK-band, NIR spectra of the object.

2.1.2. Spectropolarimetry

The linear spectropolarimetric data were obtained with the Intermediate-dispersion Spectrograph and Imaging System (ISIS) on the 4.2 m William Herschel Telescope (WHT), La Palma, during the night of 2015 August 4. Most of the data from that run were published in Ababakr et al. (2017), to which we refer for more details regarding the observations. We briefly summarise these below. A1200R grating centred at 680 nm with a spectral range of 100 nm and a slit width of 1.0 arcsec providing a spectral resolution of ∼35 km s⁻¹ was used. The polarisation optics, which consist of a rotating half-wave plate and a calcite block, were employed in order to perform linear polarisation observations. The total integration time was 1800 s for IRAS 17163. Polarised standard stars and zero-polarised standard stars were observed to calibrate for the instrumental polarisation and angle offset, both of which were found to be negligible.

The data reduction was carried out using IRAF (Tody 1993). The extracted spectra were imported into the TSP package (Bailey 1997) to compute the Stokes parameters. Finally, the percentage of linear polarisation, P, and polarisation angle, θ, were obtained. The spectrum covers the wavelength range 635–720 nm. Owing to a combination of the redness of the source and response curve of the optics, the S/N of the spectrum ranges from ∼100 blueward of Hα to larger than 300 at redder wavelengths. As the count rates increase very steeply in the first 10 nm, the polarisation values at these wavelengths are less precise than those longward of around 650 nm.

2.2. Interferometry

2.2.1. GRAVITY observations

IRAS 17163 was observed for a total of nine hours during ten nights between April and August 2017 with GRAVITY (Gravity Collaboration 2017; Eisenhauer et al. 2011) using the four 1.8 m Auxiliary Telescopes of the VLTI. GRAVITY is an interferometric instrument which operates in the K band; it combines signal from the beams of the four telescopes, delivering interferometric observables at a range of spectral resolutions. To achieve an optimal image reconstruction, the data were taken at three different baseline configurations: small (A0-B2-D0-C1), medium (K0-G2-D0-J3, A0-G2-D0-J3), and large (A0-G1-J2-J3, A0-G1-J2-K0) resulting in a good filling of the uv-plane for a four-element interferometer. The uv-plane coverage is shown in Fig. 1.

Fig. 1. Coverage of uv-plane of VLTI/GRAVITY and AMBER observations of IRAS 17163.

The observed wavelength range covered the NIR K band (1.99–2.45 μm) and all the observations were performed in high-resolution (HR) combined mode, which provides a spectral resolution of R ∼ 4000. This corresponds to a velocity resolution of 75 km s⁻¹. The observed configurations combined provided projected baseline lengths B between ∼10 m and 130 m, corresponding to a maximum angular resolution of λ/2B ∼ 1.7 mas at the wavelength of Brγ, which corresponds to a spatial scale of ∼2 au at 1.2 kpc. The position angles (PAs) of the observed baselines covered the entire range between 0 and 180°. The technical overview of the observations including the integration times and atmospheric conditions (i.e. coherence time, seeing) is given in Table A.1.

The red giant star HD 151078 (RA = 16^h46^m48^s, Dec = −39°22′37$\stackrel{″}{.}$1 [J2000]), was observed as a standard calibration object within 30 min of the science object and under relatively stable weather conditions. The spectral type of HD 151078 is K0III, the K magnitude is 3$\stackrel{\text{m}}{.}$2, and it is characterised by a uniform disc diameter of 1 mas (JMMC SearchCal; Bonneau et al. 2011). This star was selected as calibrator for all configurations and observing nights. This object was also used as a telluric standard and during the spectra normalisation process. For the reduction and calibration of the observations the GRAVITY standard pipeline recipes (as provided by ESO, version 1.0.5) were used with their default parameters. Before correcting the telluric absorption features in the spectrum of HD 151078, we removed the absorption features including the CO band in the spectra of the calibration standard object. This could be done by dividing the standard spectra with a telluric corrected spectrum of a standard object of same spectral type (HD 105028) after having it rebinned and having applied a velocity shift to match the resolution and velocity of our observations. The K-band template spectrum of HD105028 was obtained from the NIFS sample of the on-line library of the GEMINI observatory¹ and covers a wavelength range of 2.04–2.43 μm, which is almost the entire range of our spectra.

2.2.2. AMBER observations

The interferometric AMBER observations were performed in visitor mode on the night of 2012 June 22. AMBER is a decommissioned instrument of the VLTI that used to operate in the NIR H and K bands as spectro-interferometer, combining beams from three telescopes. Our observations were performed using the two baseline configurations, D0-I1-H0 and D0-H0-G1, delivering projected baselines between ∼30 and ∼70 m and PAs between ∼5 and ∼170° (see Fig. 1 and Table A.3). The observations were made using the high spectral resolution mode (R = 12 000). The obtained spectra cover a wavelength range between 2.15 and 2.19 μm, which includes the Brγ in emission similar to our GRAVITY observations but not the Na I and Mg II emission. The technical overview of the observations including the integration times and atmospheric conditions (i.e. coherence time, seeing) is given in Table A.3.

The dataset was reduced using the amdlib software (version 3.0.9 Tatulli et al. 2007; Chelli et al. 2009). During the data reduction we applied frame selection criteria to limit the data quality degrading effects from which AMBER data generally suffer. We therefore chose to work with the best 20% of the observed frames; this approach was previously suggested to provide consistent visibilities (Malbet et al. 2007). We selected a piston cut-off of 100 μm since for HR data only large piston values (> 100 μm) affect the visibility measurements. This process is described in Tatulli et al. (2007).

Owing to overall adverse weather conditions (seeing > 1″, coherence time < 5 ms; see Table A.3), an absolute calibration could not be performed by only making use of the calibration star as a standard method. Instead, we used our GRAVITY observations of the science object to estimate the emitting size of the continuum at 2 μm and we used the predicted continuum visibilities of a uniform disc at the corresponding AMBER baselines and PAs to calibrate the AMBER visibilities. We note that the observations are 5 yr apart, and during this period a typical YHG is characterised by significant quasi-period photospheric pulsations, resulting in differences in its radius by up to 50% (e.g. Lobel et al. 1994; de Jager 1998). As we show later in the paper (Sect. 5.1), the 2 micron emission originates from the stellar object directly, corresponding to the stellar size. Therefore, some degree of uncertainty and bias is to be expected during the calibration process of the AMBER dataset.

Although the observations are ∼5 yr apart, the continuum size at 2 μm corresponds to the stellar size and is not expected to vary significantly within this period. Therefore, the time difference alone does not introduce significant uncertainties and bias in this calibration process.

3. Observational results

3.1. K-band spectrum

We present a new NIR spectrum of IRAS 17163. Our observed wavelength range includes the emission lines of the Brγ hydrogen recombination at 2.167 μm, the Na I 2.206 μm and 2.209 μm doublet, and the Mg II 2.137 μm and 2.144 μm, while there are no signs of CO lines in emission or absorption. The final averaged and corrected spectrum of IRAS 17163 centred on the lines of interest (Brγ, Na I and Mg II) is presented in Fig. 2, overplotted with the X-shooter spectrum.

Fig. 2. Normalised spectra of IRAS 17163 (GRAVITY and X-shooter) around the 2 μm continuum level. The presented wavelength range contains all the observed Brγ, Na I, and Mg II emission lines. The dotted line represents the normalised telluric absorption spectrum applied to correct the spectrum of the object. The observed difference in line-to-continuum ratio is due to the higher spectral resolution of the X-shooter.

Thanks to the GRAVITY observations, we have a time series of five months which we can directly compare. Brγ shows a P-Cygni profile and its peak intensity gradually drops with time, while Na I emission shows only emission and does not vary significantly given the uncertainties of ∼13% (Fig. 3). Given that our spectra show clearly a P-Cygni profile we can conclude that Brγ is spectrally resolved. In particular, for the first observing date (2017-04-11) the emission part of Brγ is characterised by a line-to-continuum ratio of 1.32 at a peak velocity of 57 km s⁻¹ (LSR) and an observed full width at half maximum (FWHM) of ∼72 km s⁻¹, which is close to the instrumental resolution (∼75 km s⁻¹). The dataset shows a gradual drop in the peak of the Brγ emission (up to ∼67%) over the observing period of four months. In particular, during the last observing date (2017-08-28), Brγ emission is characterised by a line-to-continuum ratio of 1.11 at a higher peak velocity of 68 km s⁻¹ and a FWHM of ∼64 km s⁻¹. The strongest line of the Na I doublet is broader by almost a factor of two (∼130 km s⁻¹) compared to Brγ, with a line-to-continuum ratio of 1.21 at a peak velocity of 27 km s⁻¹. The strongest line among the Mg II shows a FWHM of ∼100 km s⁻¹ and a line-to-continuum ratio of 1.13 at a peak velocity of 24 km s⁻¹, similar to Na I. The observed lines in X-shooter spectrum are narrower (∼55 km s⁻¹ for the atomic lines) and have a slightly larger line-to-continuum ratio than in the GRAVITY spectrum, which can be explained by the lower spectral resolution of the GRAVITY data. The line equivalent widths (EW) of the two datasets are similar to within 10%, which is the typical error in EW determinations. This suggests that either no extended flux is lost in the GRAVITY observations and that therefore the object is fully within the field of the 1.8 m AT (roughly 200 mas) or that any flux loss is the same for the continuum and the lines.

Fig. 3. Spectra of the Brγ (left) and Na I doublet (right), normalised to the continuum level, for the different observing dates and configurations in different colours. The dataset shows a gradual drop in the flux of the Brγ emission (up to ∼67%) over the observing period of four months, while Na I emission does not vary significantly.

The resulting peak velocities of Na I and Mg II are 9 km s⁻¹ and 6 km s⁻¹ higher than the systemic velocity reported by Wallström et al. (2015) (LSR; 18 km s⁻¹), who use the Fe II emission lines in the optical. Wallström et al. (2017) use the H30α recombination line, but at poorer spectral resolution and find a very similar velocity of 21 ± 3 km s⁻¹. The observed difference between our measurements and those of Fe II lines are most likely due to the low spectral resolution of GRAVITY compared to HERMES (High Efficiency and Resolution Mercator Echelle Spectrograph) spectrograph on the 1.2 m Mercator telescope (4 km s⁻¹, ∼75 000; Raskin et al. 2011). The larger observed velocity shift of Brγ is most likely due to its P-Cygni profile which results in a blue-shift absorption and a peak emission at longer wavelengths. Therefore the comparison of Brγ peak velocities with those derived from the other lines has little meaning. The spectrum of Brγ shows some indications for absorption at longer wavelengths, which needs further investigation. The observed line profiles in velocity space can be seen in Fig. 4.

Fig. 4. Line profiles of the observed Brγ, Na I, and Mg II emission towards IRAS 17163. The dotted line represents the telluric absorption spectrum applied to correct the spectrum of the object. The vertical line represents the systemic velocity as reported by Wallström et al. (2015) using the optical Fe II emission.

Previous studies have explained spectroscopic variability in P-Cygni Hα profiles as a result of variability in mass-loss rate (de Groot 1987) accompanied by possible variations in the velocity field (Markova et al. 2001). In particular, a lower mass-loss rate would lead to a more transparent envelope and therefore would give us access to the hotter regions close to the star. With significant variations in mass-loss rates though, we would expect variations in Na I emission profiles, which is not what we observe in the time range of our GRAVITY observations (four to five months). The first and thus far only NIR spectrum of the object was published by Lebertre et al. (1989) and has a spectral resolution of R ∼ 1600 in the wavelength window of 2.00–2.40 μm. As reported and seen in Fig. 9 from Lebertre et al. (1989), the only clear emission was Brγ. A closer look of that spectrum reveals some weak Mg II and Na I emission. Mg II may appear weaker than what we see in our spectrum because of the lower spectral resolution. The fact that the ratio of Brγ over Na I was clearly higher 30 years ago, however, and therefore Na I emission has become stronger, could indicate the presence of significant variations in mass-loss rates, i.e. a recent mass-loss episode. This variability appears to be more significant in a timescale of few decades but not within few months.

The Na I doublet has been previously observed in emission towards similar objects (e.g. IRC+10420; Thompson & Boroson 1977; Oudmaijer & de Wit 2013) while Mg II has not been previously observed towards objects of this class (e.g., Oudmaijer & de Wit 2013). Clark & Steele (2000) presented K-band spectroscopy of Be stars, where they reported the presence of Brγ and Mg II emission with the absence of He I and, in only few cases, also the presence of the Na I doublet emission for stars of B2-B4 spectral type. Understanding the origin of the Na I doublet emission has been a main challenge in various studies. In Sect. 7.1 we investigate several scenarios to explain the observed Na I doublet emission.

3.2. Continuum polarisation

To arrive at broadband overview of the optical spectropolarimetry, we binned the spectrum in steps of 10 nm from 640 nm onwards and show the resulting polarisation and position angle in Fig. 5. IRAS 17163−3907 shows a very strong continuum polarisation of order ∼10.5% at PA 25°. The polarisation declines by almost 1% over the spectral range, whereas the PA varies within 1°.

Fig. 5. Polarisation (bottom panel) and polarisation angle (top panel) observed towards the Fried Egg Nebula. The black dots represent the broadband values, computed in steps of 10 nm, or 100 Å, from the spectropolarimetric data. The red triangles represent the Lebertre et al. (1989) data. The dashed line is the best-fit Serkowki law for interstellar polarisation.

The polarisation can be fitted with a Serkowski law (Serkowski et al. 1975). This describes the empirical wavelength dependence of the polarisation due to interstellar dust with only three parameters: the maximum polarisation, the wavelength at which this polarisation is maximum, and, k, a measure for the shape of the polarisation behaviour. Keeping the latter at its typical value of k = 1.15, we arrive at our best fit for P_max = 12.12 ± 0.4% at λ_max = 469 ± 23 nm, which is shown in the Fig. 5 as a dotted line. The Serkowski law represents the observed polarisation very well. Given the large interstellar extinction towards the object and the fact that the observed polarisation can be described with the Serkowski law, we can infer that a large fraction of the observed polarisation is due to foreground material. As pointed out by for example Oudmaijer et al. (2001), 1% in polarisation corresponds to around one magnitude in A_V, so the polarisation implies an interstellar A_V of order 10 magnitudes. This is consistent with the value we determine in Sect. 4.3 based on the spectral energy distribution (SED). We note that the large interstellar extinction is supported by both CO observations, which reveal the existence of a foreground dark cloud (e.g. Dame et al. 2001; Wallström et al. 2015), and the strong diffuse interstellar bands (DIB) presented in Lagadec et al. (2011).

Lebertre et al. (1989) also measured a high level of polarisation of 12.7%, 12.1%, and 11.5% at angles 27°, 26°, and 24° in the V, R_C, and I_C bands. It would appear that the polarisation at 650 nm (the central wavelength of the R_C Cousins band) has changed by ∼1.5%, but the angle stayed constant between 1988 and 2015. In both datasets the polarisation decreases towards longer wavelengths, which is expected from Serkowski’s law. However, we note that the PA clearly changes with wavelength in the 1988 data which cover a large wavelength range. Our R-band data seem to have the same slope at the longer wavelengths. This rotation is an indication of an additional polarising agent at work in addition to the interstellar polarisation, and this could be potentially be intrinsic to the source. This notion is reinforced by the fact that the polarisation is variable, although, ideally, multi-epoch measurements should be made with similar equipment to investigate this properly.

3.3. Line spectropolarimetry

To examine the polarisation across the Hα emission line, we present the polarisation spectrum in Fig. 6. This figure shows the intensity and Stokes Q, U vectors as a function of wavelength. The data are adaptively binned such that the statistical error in every wavelength step is the same. Normally, we would plot the polarisation and polarisation angle as function of wavelength², but in this particular situation, the polarisation changes across the Hα line are substantially more clearly defined when the Stokes Q, U vectors are shown. This is because the wavelength independent (across Hα) interstellar polarisation adds a constant scalar to the Q, U values, such that they retain their original appearance across Hα, whereas the polarimetric line profile is more affected (see equation in footnote; or e.g. Figs. 2 and 7 in Ababakr et al. 2017 and Oudmaijer et al. 1998 respectively).

Fig. 6. Hα spectropolarimetry of IRAS 17163. The data are presented as a triplot of the Q, U, and intensity Stokes vectors (from top to bottom). The data are rebinned to a constant error in polarisation of 0.16%. The excursion across the Q vector alone indicates an asymmetric scattering geometry at a projected angle of 90° across the sky.

The continuum values of Q and U change in such a way that the total polarisation decreases with wavelength, as outlined above. Stokes Q displays a clear and significant change across the Hα emission line, while Stokes U does not show a clear effect; the excursion that may be present in the few pixels around Hα reaches at most the 2 sigma level. A change in Q alone would mean a change in polarisation along a 0° angle on the sky, which seems to be the case here. Formally the intrinsic polarisation angle must then be computed using the change in Q and U, and can be derived using Θ = 0.5 arctan(ΔU/ΔQ); however, the relatively large uncertainty in measuring any change (if at all) in the U vector across Hα means that we cannot measure ΔU with certainty and a large formal error bar on the angle of 0° would result. A much longer exposure would be necessary to arrive at a more accurate value. As shown below the crude measure agrees with what we observe in the interferometry.

The fact that the polarisation across Hα is different from the neighbouring continuum, excludes any explanations for the polarisation that vary smoothly with wavelength, such as scattering by dust. The typical explanation for this so-called line effect takes into account that the star is surrounded by free electrons in ionised gas. The continuum photons emerging from the stellar photosphere are polarised by these free electrons. The Hα emission line photons emerge from a larger volume and as a consequence they encounter fewer free electrons resulting in a lower net polarisation. If the geometry of the ionised region is circular on the sky, both line and continuum photons have zero net polarisation. When the geometry deviates from spherical symmetry, a depolarisation across the emission line can be observed. The resulting polarisation angle is then perpendicular to the major axis of the asymmetric region on the sky. As the polarisation is largest when the electron densities are highest, the line effect is sensitive to scales often of order stellar radii or less. The technique has been successfully used to probe asymmetrical structures at small scales around many types of stars, for example classical Be stars (Clarke & McLean 1974; Poeckert & Marlborough 1976), Herbig Ae/Be stars (Oudmaijer & Drew 1999; Vink et al. 2002; Ababakr et al. 2017), and evolved stars (Trammell et al. 1994; Davies et al. 2005; Patel et al. 2008).

If electron scattering is the explanation for the continuum polarisation and change in polarisation over the Hα line, we can infer the orientation of the asymmetric material, which would be perpendicular to the intrinsic polarisation angle of 0°. In most of the references above, this asymmetric material reflects the presence of discs or outflows. However, current spectropolarimetric data do not allow us to distinguish between the two scenarios. We can say that the spectropolarimetry suggests the presence of an electron-scattering region at scales of order stellar radii and an elongation of this material in the east-west direction.

3.4. Visibilities and phases

We extracted the visibilities, differential, and closure phases from all observed configurations and PAs around the Brγ (GRAVITY and AMBER), Na I, and Mg II (only GRAVITY) towards IRAS 17163. For direct comparison, Fig. 7 presents the interferometric observables of GRAVITY as a function of baseline length of Brγ, Na I and Mg II emission, only for the large configuration. Similar plots for all configurations can be found in Figs. A.1–A.4. Figure 8 is similar to Fig. 7, but only around the Brγ, as the wavelength range of our AMBER observations did not cover the Na I and Mg II wavelengths. Therefore, our AMBER observations can inform us regarding Brγ emission and compare it with the continuum, while the GRAVITY observations also inform us about the additional emission from Na I and Mg II. We also note that AMBER combines three telescopes instead of four resulting in a single closure phase.

Fig. 7. Representative example of relative flux, visibility, differential, and closure phase as a function of wavelength around the Brγ and Na I doublet emission towards IRAS 17163 using GRAVITY on the large configuration. The 6 baselines and 4 triplets of A0-G1-J2-K0 configuration are plotted with different colours. Both Brγ emission and Na I show a smaller visibility than the continuum while Brγ shows a larger drop. Mg II does not show changes in visibility or phases. A differential phase change is also detected towards Brγ but not towards Na I.

Fig. 8. Representative examples of relative flux, visibility, differential, and closure phase as a function of wavelength for the AMBER D0-H0-G1 and D0-I1-H0 configurations around the Brγ emission towards IRAS 17163.

Brγ shows a different drop in visibility value compared to the continuum for different baselines. In particular it shows deeper drops in visibilities compared to the continuum and the other lines, while it also shows changes in differential and closure phases along the majority of the baseline PAs. The detected changes are indicative of an emitting region that is more extended than the continuum-emitting region, while the changes in closure phases suggest that Brγ emission stems from an asymmetric region. The line-emitting region (Brγ) shows a decrease in visibility of 28–35% for PAs between 80 and 100° (E–W direction) compared to the other observed directions for the longest baseline (∼130 m). This decrease is an order of magnitude larger than the corresponding errors (∼1–5%, including the uncertainty due to the calibrator and transfer function), and therefore we attribute the observed difference in geometrical effects, and in particular an elongation of the emission towards E–W direction. No asymmetries in the ionised gas were previously reported towards the object, but this interferometric finding is consistent with our spectropolarimetric results in Sect. 3.3.

The observations reveal significant differences between the Brγ emission and the line emission from Na I and Mg II. Na I shows a smaller drop in visibilities compared to Brγ, but no changes in differential phases and closure phases for the majority of the baselines and PAs. In addition the visibility drop increases as a function of baseline length for the continuum and the Na I emission and independently of PA orientation, indicative of symmetric emission along both axes. The observations suggest that Na I stems from a somewhat larger region compared to the continuum but is a smaller region compared to Brγ, and the emission is mostly symmetric at the resolution of the observations. The absence of changes both in visibilities and differential phases along Mg II indicate that the line-emitting region follows the distribution of the continuum-emitting region, while the absence of changes in closure phases point towards a symmetric structure of the emitting area. Therefore, we do not include Mg II in our more detailed analysis assuming that our geometrical results about the continuum emission also apply for Mg II.

4. Fundamental parameters of IRAS 17163

4.1. Spectral type

Before proceeding to a detailed modelling of the Fried Egg Nebula, it is important to put constraints on the effective temperature of its central star, IRAS 17163. Thus far, it had not been possible to determine the spectral type in a classical manner for the object; optical spectra up until now were only usable beyond around 500 nm, where the spectrum is dominated by emission lines, while it is the properties of the absorption lines that determine the spectral type. Moreover, spectral classifications are traditionally carried out in the blue spectral region. The blue part of the optical spectrum obtained in this work is dominated by absorption and can serve as a basis to determine the spectral type the object, while at the same time facilitate a comparison with the spectrum of the known post-RSG IRC+10420, whose emission spectrum had been shown by Wallström et al. (2015) to be very similar to that of IRAS 17163 at longer wavelengths. Lastly, the 2017 spectrum should also serve as a reference for future studies into any temperature evolution.

To carry out the spectral type classification, we downloaded high spectral resolution data of a number of supergiants from the Bagnulo et al. (2003) spectral database. The spectrum of IRAS 17163 and those of a number of A-type supergiants are presented in Fig. 9 in the typical spectral type classification wavelength range, 380–460 nm. Apart from the strong DIB at ∼440 nm betraying a large interstellar extinction, emission filling in the hydrogen absorption lines, and the lines longward of 440 nm affected by emission, the spectrum fits in naturally with the standard-type supergiants. Indeed, it can be seen that the object mostly follows the clear progression of the spectral types when placed between the A6- and A3-type supergiants with a preference for a slightly earlier than A6Ia spectral type. With temperatures for A3I and A5I supergiants listed by Straizys & Kuriliene (1981) of 8892 K and 8298 K, respectively, we would infer a temperature of about 8300–8500 K for the central star of the Fried Egg Nebula.

Fig. 9. Spectra of IRAS 17163 with a number of spectral standards taken from the UVES spectral archive (Bagnulo et al. 2003). The spectra are normalised to a continuum value of 1 and offset, for clarity, in multiples of 0.9.

A further clue to the temperature of IRAS 17163 can be gleaned from the comparison with the well-known post-RSG IRC+10420. Figure 10 shows a part of the blue absorption spectra of both objects (the IRC+10420 spectrum taken from Oudmaijer 1998). It can be seen that the stars are remarkably similar. Klochkova et al. (1997) determined a temperature of 8500 ± 250 K for IRC+10420, based on intensity of Fe I absorption lines. In turn we can infer a similar temperature for IRAS 17163 as well. In the following sections we adopt a temperature of 8500 K for the central star of the Fried Egg Nebula.

Fig. 10. Overlay of part of the blue spectrum of IRAS 17163 (black line) and IRC+10420 (red line). There is a tantalising similarity between the two objects, as previously also noticed by Wallström et al. (2015) for the optical spectrum above 5000 Å, which is dominated by emission lines.

4.2. Distance

The distance of IRAS 17163 has been a debated topic in recent decades (see also Sect. 1), but its accurate determination is crucial to constrain the fundamental parameters of both the stellar object (size, luminosity, class) and its surrounding shells (size, mass-loss rates, kinematic timescales). In this work we take advantage of the new Gaia parallactic measurements (DR2; Gaia Collaboration 2018) and derive a distance of $d=1.2_{-0.2}^{+0.4}$ kpc, which we adopt in the entire paper. To do so we applied the method proposed by Bailer-Jones et al. (2018), which uses a Galactic model to weakly inform a distance prior to be used in the determination of the distances, and not just the inversion of the parallax (0.83 ± 0.17 mas), as the latter is found to introduce more biases, failing to treat properly nonlinear effects. For the calculation we made use of the ARI Gaia services³.

Given that the specific assumptions of this method do not necessarily apply to this object, we also determined the distance using the direct parallax inversion and found it to be 1.2 kpc, which is in perfect agreement with the distance found using the method suggested by Bailer-Jones et al. (2018).

Gaia parallax presumes that the measured star is a single and point-like source. Our analysis in Sect. 5.1 suggests that the K-band continuum and hydrogen line emission are confined to the scale of 1–1.9 mas, and thus point-like with respect to the Gaia resolution (0.1″). It also appears that most, if not all, of the K-band light is contained within the field of view. Thus the resolved sizes actually support a point-like appearance for Gaia, which gives us a further confidence in using the parallactic distance estimation.

4.3. Dereddened photometry

The photometry of IRAS 17163 shows strong excess emission beyond 5 μm. A black-body emission assuming a temperature of 8500 K (see Sect. 4.1) can reproduce the observed dereddened photometry at wavelengths < 5 μm, for which the observed flux is solely due to the central heating source. We fit a Bosz stellar model (Mészáros et al. 2012) to the photometry of 2.2 μm and below (including bands: B, V, J, H, K, and L; Epchtein et al. 1987; Lebertre et al. 1989; Cutri et al. 2003). For the Bosz model we adopted an effective stellar temperature T_eff, ⋆ of 8500 K, a surface gravity log g of 2, which is the lowest available value in the models for the adopted temperature, and a solar composition (metallicity [M/H], carbon, and alpha-element abundances). This temperature is an upper limit given the spectral type of IRAS 17163 and the interferometric continuum size (Sect. 4.1). The observed photometry is affected by extinction even at shorter wavelengths and therefore needs to be dereddened. We used a least-squares minimisation method and found a total extinction of $A_V=12.0_{-0.33}^{+0.8}$ (Fig. 11) after applying the selective visual extinction R_V-dependent law from Cardelli et al. (1989) for R_V = 3.1, which corresponds to a Galactic interstellar medium (ISM) dust (Fitzpatrick & Massa 2009), and zero point magnitude from Bessell (1979). The observed and dereddened photometry is presented in Table A.4.

Fig. 11. Dereddened IRAS 17163 photometry with total extinction AV = 12.0 to fit a Teff = 8500 K Bosz stellar model.

Our dereddened photometry, when compared to that by Hutsemékers et al. (2013), shows that the integrated energy flux per unit wavelength of the star in that study is smaller by a factor of ∼10 at a distance of 4 kpc, and therefore either the luminosity of the object is ∼10 higher than previously thought, or its distance is ∼3 times closer. This was an oversight that was propagated in more studies, since both Lagadec et al. (2011) and Hutsemékers et al. (2013) use the photometry presented in Lebertre et al. (1989) as a reference, who had not corrected for reddening. In conclusion, also based on Sect. 6.4, integrating the best-fit SED after treating the observed photometry (i.e. dereddening) properly in combination with the closer distance of IRAS 17163 at 1.2 kpc, results in a luminosity of L = 5 × 10⁵ L_⊙, and therefore IRAS 17163 maintains its class as a YHG. Given its luminosity and effective temperature, we estimate its size to be R = 325 R_⊙. Assuming a lower mass limit of M = 25 M_⊙ (M = 25 − 40 M_⊙; Lagadec et al. 2011), the resulting surface gravity log g is ∼0.9, which is a typical value for a YHG. Although this is lower than the adopted value of log g = 2, during our fitting process we found that for a given effective temperature, different values of surface gravity did not significantly affect the colours. Therefore the adopted log g does not affect the reliability of the derived extinction within its uncertainties.

5. Geometry of the emission

5.1. Size estimations

The present GRAVITY dataset includes long baselines up to ∼130 m and it can be used to estimate the sizes of the continuum emission, and those of the neutral and ionised gas as traced from the Na I and Brγ emission, respectively. For this purpose, we used the calibrated visibilities as measured at the spectral channel of the peak of the emitting lines and calibrated visibilities of the continuum. To derive and model the actual visibility of the lines (Na I and Brγ), deprived from the continuum contributions, we used Eq. (1), which describes the total visibility (V_{line + cont}) in terms of continuum and line visibilities (V) and fluxes (F) when dealing with multi-component sources, and solving for V_line as follows:

$$\begin{array}{c}\hfill V_{\mathrm{line}+\mathrm{cont}}=\frac{V_{\mathrm{cont}}\times F_{\mathrm{cont}}+V_{\mathrm{line}}\times F_{\mathrm{line}}}{F_{\mathrm{cont}}+F_{\mathrm{line}}}·\end{array}$$(1)

For consistent comparison of the sizes of the line and continuum emission, we fit the observables with simple geometrical models of a Gaussian-shaped emitting region for a range of sizes (0.2–5 mas). We find that the best-fit results in a FWHM of 1.09 ± 0.01 mas for the continuum, 1.20 ± 0.01 mas for the Na I and 1.90 ± 0.01 mas for Brγ. Figure 12 shows the observed values overplotted with the corresponding best-fit models for the three cases.

Fig. 12. Visibilities at the peak of the Brγ, Na I emission lines and continuum emission towards IRAS 17163 overplotted with simple geometrical Gaussian models. Na I and continuum emission follow very similar distribution. Na I shows a smaller drop in visibilities compared to Brγ, indicating that Na I stems from a smaller emitting area compared to the Brγ. The AMBER dataset was rebinned to match the spectral resolution of GRAVITY. A period of 5 yr indicates a shrinkage of the Brγ emitting region to half its size.

Taking into account the extreme upper temperature limit of T_dust ∼ 2000 K (pyroxene; Kobayashi et al. 2011) at which dust sublimation occurs, and the expected radius and temperature of IRAS 17163, we find that dust cannot survive at distances < 7 mas ($d=R_{\star }{T}_{\star }^2/2T_{\mathrm{d}}^2$) which is larger than the size of Na I and continuum emission (∼1 mas) by a factor of 7. We point out that dust sublimation of crystalline and oxygen-rich dust occurs at lower temperatures (∼1200 K; Blasius et al. 2012), which is a more realistic approach to IRAS 17163. Therefore 7 mas is only a lower limit for the dust sublimation radius. We conclude that the continuum emission at these wavelengths (2.0–2.4 μm) originates directly from the stellar object. This finding is also supported by the photometry findings at shorter wavelengths in Sect. 4.3, where a black body could fit the photometry at wavelengths < 5 μm.

Based on our findings in Sect. 4.3 the expected size of the central star is 2.5 mas (diameter; ∼650 R_⊙) at the distance of 1.2 kpc. The estimated size of the continuum is 1.09 ± 0.01 mas and therefore at first there is a notable discrepancy between the calculated and measured size. Given the geometrical complexity of stellar atmospheres of evolved stars, this difference is most likely due to our simple model of Gaussian distribution, which was selected for a consistent comparison with the line emissions. When we adopt a uniform disc to model the interferometric observables, the continuum size is 1.73 ± 0.02 mas, which is closer to the expected. The intensity profile that characterises such objects is not a uniform disc nor a Gaussian distribution, and it is strongly characterised by limb-darkening signatures, and even more complex structures are predicted by theory (e.g. convective patterns Chiavassa et al. 2011).

Having inferred the effective temperature in Sect. 4.1, a robust determination of the physical size of the atmosphere would lead to a more accurate estimation of the luminosity of the star. This is something we could not achieve with our current geometrical models, but could be achieved with the use of more advanced atmospheric models, similar to what has been done for RSGs in the past (e.g. Paumard et al. 2014). These kinds of models are currently lacking for YHGs and this is beyond the scope of this work. For a more robust comparison, we should also consider the uncertainties in the luminosity of the source and the Gaia distance. In addition, we should note that the reported uncertainties in the size estimations represent the goodness of fit and are rather underestimated. In conclusion, considering all the aforementioned uncertainties, the measured angular diameter is in reasonable agreement with the photometric diameter.

5.2. Model-independent image reconstruction

The VLTI is a very powerful interferometer providing high angular resolution in the visible and NIR. The measurements coming from the observations are not direct images however and therefore using image reconstruction algorithms is crucial to make a full use of its imaging capabilities.

We perform and present for the first time the image reconstruction of IRAS 17163 and its close vicinity at milli-arcsecond resolution. To do so we used the Multi-aperture Image Reconstruction Algorithm software (MIRA; Thiébaut 2008). During this process MIRA fits the observed squared visibilities (power spectrum) and closure phases (bispectrum) with models, using some a priori image properties (positivity, normalisation, and regularisation) and a cost-estimate optimisation. For the cost-estimate optimisation we used the L2-L1 regularisation, which assumes the pixel positivity and compactness of the source. With only four telescopes the above constraints are necessary to interpolate properly the gaps in the coverage of the uv-plane. The maximum number of evaluations of the cost function was set to 200 with a global weight of regularisation of 10⁻¹⁰. The regularisation type was set to entropy. In our approach an exact Fourier transforms algorithm was set. A detailed description of the chosen parameters is provided in Thiébaut (2008).

The computed field of view was ∼75 mas for B_min = 11 m (2.44 λ/B_min; 4 × Θ_max to ensure we do not cut out emission), where B_min is the shortest projected baseline, and λ = 2.2 μm. Although previous studies have shown that the image reconstruction result did not significantly depend on the required starting image (Millour et al. 2009), we chose to work with multiple simple randomised start images of the same dimensions that host a point or more extended source at the centre. We ran the algorithm for a range of pixel sizes (0.1–1 mas) and we found that for a pixel size > 0.5 mas the resulting image does not change significantly, while for smaller pixels (0.1–0.5 mas) finer structures could be seen. To avoid bias introduced by the particular image model (e.g. pixel shape), the spatial resolution of the model should be well beyond the resolution limit of the observations (< 1.7 mas). As a general approach, a super resolution is usually used and the pixel size is determined by Δθ <  λ/4B_max (∼0.8 mas). To avoid the introduction of artefacts a pixel size of 0.5 mas was chosen. A bandwidth of 4 × 10⁻⁴ μm was set for the Brγ and Na I image reconstruction corresponding to the width of one spectral channel, while the bandwidth of the continuum was set to two orders of magnitude higher (0.03 μm).

Separate image reconstruction was performed for the continuum and around the central channel of the Brγ and Na I emission (Fig. 13). The image reconstruction reveals a compact symmetric emission of the continuum and Na I, while exhibiting a more extended Brγ emission along with asymmetric features, thereby confirming what we already discussed in Sect. 3.4. In particular the image reconstruction for Brγ also reveals a northern and southern component.

Fig. 13. Image reconstruction of the continuum emission (greyscale) overplotted with Brγ in red contours at [0.1, 0.3, 0.6, 0.75, 0.95] × peak and Na I in blue contours at [0.2, 0.5, 0.8, 0.95] × peak towards IRAS 17163. Brγ covers a more extended region compared to the continuum and Na I, which appear to follow very similar distributions.

In conclusion we present a model-independent image reconstruction of this object at 2 μm for the first time. Although we treated this delicate procedure with extra care, the observed flux changes of the Brγ at different epochs introduce more uncertainty in the reconstructed image by possibly affecting the observed visibilities. The channel maps of the Brγ emission, together with the combined (AMBER and GRAVITY) image reconstruction, using three different software packages is less robust and is presented in Figs. B.1 and B.2. Furthermore, Fig. 12 shows that the object is formally only partially resolved (i.e. does not reach the first zero of visibility), and therefore at the smallest observed scales the image suffers from higher uncertainty. Despite the uncertainties introduced by the special features of this unique dataset, we find that the main spatial structures confirm our findings when applying geometrical models to the interferometric observables, showing a more extended and asymmetric Brγ emission than that of the continuum and sodium.

6. Modelling the shells of the Fried Egg Nebula

6.1. Modelling the warm dust emission: Introduction

We model the warm dust emission in various wavelengths to obtain quantitative information on the recent mass-loss history as traced by the dust and to derive the dust mass and morphology of the distinct dust shells surrounding IRAS 17163. In previous works, both Lagadec et al. (2011) and Hutsemékers et al. (2013) assume a distance of 4 kpc and fit the available photometry with radiative transfer models to estimate the mass within their observed detached shells. Lagadec et al. (2011) focusses on fitting a small volume encompassing two warm shells observed in mid-IR wavelengths using VISIR (the Fried Egg Nebula) with a radius of ∼0$\stackrel{″}{.}$45–2$\stackrel{″}{.}$5. Hutsemékers et al. (2013), on the other hand, focusses on fitting a larger cooler volume encompassing of both the Fried Egg Nebula as an unresolved single shell, and a much larger shell revealed at the far-infrared, 70 μm, and 160 μm images using PACS/Herschel with a radius of ∼18″–40″. We note that although both studies achieve a good fit to the SED (apart from Hutsemékers et al. (2013) to the optical part), the resulting total dust mass for the Fried Egg Nebula differs by an order of magnitude.

The degeneracy can be broken by not only fitting the SED (i.e. constraining the optical depth of each shell) but the intensity profiles as well, which constrains the position and thickness of the separate shells and the dust temperature. In this study we adopt a model to simultaneously fit (a) the dereddened photometry with a newly determined value of visual extinction based on our revised distance and (b) the VISIR/VLT infrared radial profiles of the star and the surrounding Fried Egg Nebula, breaking the SED degeneracy.

6.2. Dust modelling: Initialising 2-Dust

In this work we used the 2D radiative transfer code, 2-Dust (Ueta & Meixner 2003), towards the Fried Egg Nebula to fit simultaneously for the first time both the SED and the flux maps (radial profiles) of chosen wavelengths (8.59 μm, 11.85 μm, and 12.81 μm), which allows us to spatially constrain the observed shells (R_in, R_out) and their properties (i.e. opacity, temperature, and mass). The selected wavelengths are based on the VISIR images of the object, published in Lagadec et al. (2011).

The code 2-Dust solves the equation of radiative transfer tracing the heating mechanisms of a grain located at any place in the shell radiated from both the central star and from surrounding dust radiation using a 2D polar grid. The code deals with a large range of input parameters including the fundamental parameters of the central heating object (i.e. effective temperature, size, and distance) and the properties of the dust (Sect. 6.3) that surrounds it (i.e. shell size, dust composition, and opacities).

In particular, the stellar parameters and the distance of IRAS 17163 were determined in Sect. 4, and therefore were fixed to T_eff = 8500 K, R = 325 R_⊙ (L = 5 × 10⁵ L_⊙), and d = 1.2 kpc during the entire fitting process (IRAS 17163; Table 1). In our approach we refined a polar grid of 48 radial and 8 latitudinal grid points. The log-spaced radial grids for the shells are characterised by an inner and an outer shell radius (R_in, R_out). The locations and sizes of the dusty shells were adjusted accordingly to fit the observed radial profiles. In particular, the images and radial profiles of the Fried Egg Nebula at the chosen wavelengths have previously revealed two shells with further indications of a third hot shell which had not been accounted for during the modelling process (Lagadec et al. 2011). In this work we provide evidence for the existence of a third hot inner shell that we modelled, resulting in a better fit to both the SED and the radial profiles. In particular, we modelled three shells: a hot inner shell extending between 0$\stackrel{″}{.}$3 and 0$\stackrel{″}{.}$45, an intermediate warm shell extending between 0.6″ and 1$\stackrel{″}{.}$1, and an outer cooler shell between 1$\stackrel{″}{.}$9 and 2$\stackrel{″}{.}$5, which we implemented in our model during the SED and radial profile fitting process (Fig. 14). The outermost cooler shell reported in Hutsemékers et al. (2013) of a radius of ∼18″–40″was not included in our models. The r_d, min of the intermediate and the outer shells were strongly constrained at 0$\stackrel{″}{.}$6 and 1$\stackrel{″}{.}$9, respectively, by the observed peak of each of the shells within the intensity profile using the existing radial profiles of the images at mid-IR (Fig. 14). Constraining the r_d, max of each of the shells was less robust. The exact position and size of the shells are both important in determining the temperature and the 10–18 μm flux ratio of each shell.

Fig. 14. Left: VISIR/VLT mid-IR monochromatic image at 8.59 μm (Lagadec et al. 2011). Right: radial profiles at 8.59 μm, 11.85 μm, and 12.81 μm are extracted and overplotted with the modelled profiles as described in Sect. 6.2. The size and position of each shell is indicated in light blue. The point spread function (PSF) at 11.85 μm is also overplotted for reference, revealing the hot inner shell. The PSF becomes broader at longer wavelengths.

To be able to model the three distinct shells with 2-Dust, we have to modify the density distribution by applying different values of optical depths τ_λ, at each of the individual shells. To do so, we treated each shell separately, using the output flux distribution (flux density over wavelength) of the central star and the inner hot shell as the input for the intermediate shell and so forth. The optical depth τ applied to each shell was the only free parameter in our models and was fine tuned to provide a good fit to the relative peak intensity profiles of the observed 8.59, 11.85, and 12.81 μm emission in combination with a good SED fit. Details on the fitting and the final results are presented in Sect. 6.4.

6.3. Dust composition and size

The dust in IRAS 17163 is dominated by silicates, which were first observed in the IRAS spectrum of the source as strong emission around 10 μm and 18 μm (Lebertre et al. 1989), and includes pyroxenes (Mg_xFe_1 − xSiO₃). This is indicative of an oxygen-rich composition and it is similar to what is also found towards IRC+10420, Wray 15-751, and other YHGs (Oudmaijer et al. 1996; Voors et al. 2000). In this work we adopted a dust composition of a 50/50 mixture of Mg(40)/Fe(60) and olivines (Dorschner et al. 1995) with astronomical silicates (Weingartner & Draine 1999), similar to those presented in Hutsemékers et al. (2013). We adopted the so-called MRN dust grain size distribution determined by Mathis et al. (1977), which follows the power law n(α) ∝ α^−3.5 with α_min <  α <  α_max and α_min = 10⁻³ μm. The dust grains that surround IRAS 17163 can have a maximum size α_max between 1 μm and 3 μm (Hutsemékers et al. 2013). The bulk density is fixed to 3.2 g cm⁻³. Adopting a maximum grain size of 1 or 3 μm affects the absorption and scattering cross sections. In particular, for wavelengths > few microns, the scattering cross sections are found to be > 1 order of magnitude larger for an upper limit on the grain size of 3 μm compared to 1 μm. The silicate features at 10 μm and 18 μm can be seen in both cases. Smaller grains but with constant bulk density, result in lower mass grains and therefore a larger amount of grains within one shell, which increases the total flux within the shell.

Therefore, the absorption and scattering cross sections generated by the code are dependent on the dust composition and grain size distribution, altogether influencing the shapes of the resulting SED and radial profiles for each wavelength. In our case a maximum grain size of α_max = 1 μm appeared to provide a better fit of the silicate features in the SED and the intensity profiles compared to the α_max = 3 μm. Adopting different dust types and mixtures is expected to affect the modelled emission, but such an approach is beyond the scope of the current work.

6.4. Best fit: SED and synthetic radial profiles

During the fitting process, the optical depth τ of each shell, which influences the strength of the emission, was tuned with the goal to reproduce a good fit on the relative peak intensity profiles of the observed 8.59, 11.85, and 12.81 μm emission. We started the initial fitting process by trying to reproduce the resulting mass-loss rates reported in Hutsemékers et al. (2013). Then we scaled everything according to the new distance, which meant that the optical depths had to be reduced by a factor of ∼10 in order to reproduce the updated flux of the shells. We expanded our approach to three shells and varied the values of the optical depths around a central value (up to a factor of 2), which helped us to constrain the solution that appeared to provide a good fit to the observations.

Our best fit was achieved for a τ_70 μm of 3.5 × 10⁻⁵ for the inner, 3.5 × 10⁻⁴ for the intermediate, and 3 × 10⁻⁴ for the outer shell. This corresponds to A_70 μm of 7.4 × 10⁻⁴ for all three shells, which for the adopted dust properties and composition (i.e. absorption, and scattering cross sections), translates to a total circumstellar extinction of A_V, dust ∼ 0.19 (A_V/A_70 μm ∼ 257), assuming an averaged grain mass of 2.0 × 10⁻¹⁸ g. The interstellar extinction can be then determined simply by subtracting the circumstellar extinction from the total extinction derived in Sect. 4.3, resulting in an A_V, IS = 11.8.

The SED and intensity radial profile fits are presented in Figs. 14 and 15, respectively. In Fig. 15 we present the resulting SED of the three shells individually (hot inner shell, cool inner shell, and outer shell), but also the combined SED (total spec), which we assessed for the goodness of the fit. In the same plot we present the modelled SED presented in Hutsemékers et al. (2013), which provided a good fit of the IR data points (PACS and AKARI). We note that in that work the good fit starts to work only for the wavelengths beyond 10 μm, showing a very poor fit between 1 μm and 10 μm, while the optical photometry below 1 μm was totally neglected. It is beyond the scope of the current paper to fit the outermost shell of the Fried Egg Nebula, and therefore our modelled SED does not reproduce the IR photometry at wavelengths > 70 μm. Yet given the more complex nature of our models (SED, radial profiles), the dereddening, and the completeness of the photometry, we consider it an improved fit overall compared to all previous works. The resulting parameters for each shell are presented in Table 1 together with the physical parameters of the central heating source IRAS 17163.

Fig. 15. Spectral energy distribution fit to the dereddened photometry gathered from literature (symbols). The dashed line shows the 5 × 105 L⊙, 8500 K black-body (BB) input used for the model. The dotted line is the stellar and hot inner shell emission. The dot-dashed line shows the stellar and cool inner shell emission. The large spacing dashed line shows the outer shell emission contribution. We note that our models did not include the outermost cooler dust shell (18–40″) reported in Hutsemékers et al. (2013), and therefore do not reproduce the far-infrared photometry.

For the SED dereddening we used the standard interstellar extinction curve from McClure (2009), for a total-to-selective extinction of R_V = 3.1, who found that extinction effects can be significant at wavelengths as long as 10 μm; we dereddened the data for the extinction due to the ISM. The extinction curve takes into account the silicate absorption feature from the ISM, for the A_V, IS extinction of 11.8, and therefore it is corrected in the SED.

Our synthetic images were produced with the same pixel scale as the 2008 VISIR observations of 0.075″ per pixel, corresponding to the VISIR pre-upgrade DRS detector. The observed radial profiles of the images presented in Lagadec et al. (2011) can be directly compared to the synthetic radial profiles for a given pixel scale. The number of grid points (radial and latitudinal) was also carefully selected to limit the presence of artefacts in the synthetic maps and extracted radial profiles. To extract the final synthetic radial profile for each wavelength, we added the resulted images for all three shells, and the central source was only included in the inner shell map. To match the observations the synthetic images were convolved to the angular resolution of the VISIR/VLT (θ = 1.22 λ/D), which corresponds to 0$\stackrel{″}{.}$26, 0$\stackrel{″}{.}$36, and 0$\stackrel{″}{.}$39 at 8.59 μm, 11.85 μm, and 12.81 μm, respectively for the 8.2 m VLT. The observed radial profiles are asymmetric, with the left (eastern) side of the distribution showing some brighter spots compared to the right (western) side of the distribution. In this study, for the direct comparison with the perfectly symmetric modelled radial profiles, we mimic a perfect symmetry by focussing on the right side of the radial profile, which is more consistent to the azimuthal averaged intensity profile (Fig. 14).

The modelled radial profile at 8.59 μm fits the observed profile well at 8.59 μm, but the fit is less robust with the 11.85 μm and 12.81 μm emission. In particular, a larger amount of flux is needed to fit the observed profiles, which becomes more apparent near the location of the inner hot shell. Previous low resolution spectroscopic data (LSR; Epchtein et al. 1987) were used as a guide for the 10 μm and 18 μm silicate features that helped constrain the SED, which is based on observations of 30 years ago. A larger amount of flux at 11.85 μm and 12.81 μm, and therefore an improved fit at those wavelengths could possibly be achieved by taking into account a broader silicate feature at 10 μm, but that is not supported by the LSR data. This led to a compromised fit of those longer wavelengths. Furthermore, the observed images and mass loss of the Fried Egg Nebula are not perfectly symmetric but rather show some clumpiness, which is also expected to contribute to the observed differences. In conclusion, given the above constraints and the sensitivity of the models to the initial dust composition, we consider our three shell model to be able to provide a good fit to the overall shape of both the observed SED and radial profiles.

6.5. Physical properties of the shells

The sum of the dust mass within the three shells and a canonical gas-to-dust ratio of 100 results in a total mass of ∼6.5 × 10⁻³ M_⊙ within the Fried Egg Nebula, which is about three and one orders of magnitude lower than what Lagadec et al. (2011) (∼4 M_⊙) and Hutsemékers et al. (2013) (∼0.21 M_⊙) reported, respectively. This difference is mainly a consequence of the closer distance we adopt in our study, while the adopted gas-to-dust ratio of 40 in Hutsemékers et al. (2013) could only explain a higher mass by about a factor of 2. The mass of each shell, M_d, is given by M_d = (4πτ_λr_d, minr_d, max)/κ_λ. Given the closer distance to the object, each shell is also closer to the central heating source by a factor of 3.3 compared to previous works and τ_λ is smaller by a factor of 10. Therefore, the total mass for each shell is expected to be smaller by about two orders of magnitude, which confirms what we find.

We note that both the shells of the Fried Egg Nebula and the larger scale shell are closer to the central heating object than previously thought as a result of the closer distance (Sect. 4.2). The dust temperature is given by $T_{\mathrm{d}}={(L_{\star }/16\pi \sigma r_{\mathrm{d}}^2)}^{1/4}$, where L_⋆ is the stellar luminosity, r_d is the distance of the dust grain to the star, and σ is the Stefan Boltzman constant. Therefore, a physical distance of the shells closer to the star means an increase in T_d in the shells compared to what was previously reported.

7. Discussion

In this section we discuss some of our main findings, topics of debate in other studies, both theoretical and observational. In particular, we attempt to explain the observed emission spectrum in the K band, which contains the Na I 2.2 μm doublet in emission, the origin of which has remained unexplained in previous studies. In addition we discuss the varied mass-loss rates of our three modelled distinct shells, which is a very unique finding that can help to shed more light on both theory and observation of objects in enigmatic stages of stellar evolution.

7.1. Origin of the Na I emission

7.1.1. Traditional interpretations

The observed Na I 2.2 μm doublet in emission has been reported in more studies of well-known YHGs (IRC+10420, HD 179821, HR 8752, and ρ Cas (Lambert et al. 1981; Hrivnak et al. 1994; Hanson et al. 1996; Oudmaijer & de Wit 2013), LBVs and B[e] stars (e.g. Hamann & Simon 1986; McGregor et al. 1988; Morris et al. 1996). The origin of the Na I 2.2 μm doublet towards evolved massive stellar objects has been traditionally seen as a challenge and has remained unexplained.

We consider several scenarios to explain the origin of this neutral atomic line emission towards IRAS 17163. The neutral Na I emission arises from a region that is 10% larger than the continuum source and is in turn smaller than the region responsible for the emission component of the Brγ profile. Oudmaijer & de Wit (2013) find a similar situation for IRC+10420. We would expect a neutral metal (Na I) to stem from a region further out than that of the ionised hydrogen. Since this is clearly not what we observe, the challenge is to determine which situation can give rise to this apparent contradiction. In this section, we consider whether Brγ emission arises through recombination in a more extended, ionised zone, and Na I from the photosphere instead, or whether a mostly neutral atmosphere can account for all the emission lines.

One of the scenarios is that of a disc. A dusty disc distribution, where both Na I and continuum are distributed, allows for shielded Na I. The observed more extended Brγ emission could then be explained by the penetration of UV/optical radiation at longer distances away from the disc in a low density halo. In Sect. 5.1 we argued that dust cannot survive at distances < 7 mas, unless it is self-shielded. In that case, we would expect some additional emission from the dust present in the SED at ∼2.0 μm, which is not what we see. We conclude that the continuum emission at 2.0 μm originates from the stellar object directly and thus there are no self-shielded regions close to the star that can provide a place for Na I to exist. Therefore, the scenario of a dusty disc distribution at ∼1 mas (emitting sizes of continuum and Na I) can be ruled out.

Another scenario previously proposed to explain that kind of line emission spectrum is the presence of a pseudo-photosphere. Oudmaijer & de Wit (2013) examine the scenario of an optically thick wind in the case of IRC+10420, which forms a pseudo-photosphere and acts as a shielding mechanism for the sodium, preventing its ionisation. This scenario was also supported by Humphreys et al. (2002) based on the high mass-loss rate of ∼2 × 10⁻⁴ M_⊙ yr⁻¹ of IRC+10420. As we show below, the region within 1.09–1.90 mas around the central star cannot be dilute enough for a nebular model to apply and at the same time explain the measured line fluxes. Therefore, this scenario can probably be excluded.

7.1.2. Two-zone local thermodynamic equilibrium model

For a proper interpretation of the 2.1–2.4 μm emission line spectrum, it turns out not to be necessary to assume that the H Brγ line emission arises by recombination of H⁺ in an ionised gas. In the following we show that a simple model based upon a uniform disc or spherical shell in strict local thermodynamic equilibrium (LTE) reproduces the observations fairly well. Strict LTE means that the ionisation balance is described by the Saha equation for all the abundant elements and that the energy levels of all these atoms and ions are populated according to the Boltzmann equation at the same temperature. Such a simple model is constructed as follows. Consider a face-on disc (or spherical shell) of radius r and thickness L = ρr with an angular diameter θ = 2r/D, where ρ is a dimensionless constant and D is the distance to the source. The equilibrium is specified by a set of elemental abundances relative to hydrogen, temperature T, and number density of hydrogen in all forms (H, H⁺, H₂, H⁻), n_H in cm⁻³. For reference, solar photospheric abundances of the 17 most abundant elements are assumed (cf. Asplund et al. 2009). The abundances of all neutrals, first ions, and a few molecules (H₂, CO) are computed in LTE at the specified conditions (T and n_H). The peak optical depths and flux densities of the infrared lines of interest are then calculated, given the observed line widths. The continuum opacities and fluxes are computed at the wavelengths of the lines, taking into account electron-ion and electron-neutral free-free opacity (bremsstrahlung), H⁻ bound-free, and Thomson scattering. A series of such models spanning some range in density and temperature was computed to find the conditions that closely match the observed properties of the Na I, Mg II, and H Brγ lines in the K band. The adopted observed parameters are summarised in Table 2.

The peak fluxes $f_{\nu }^{\mathrm{obs}}$ and line widths (FWHM) were determined from Gaussian fits to the X-shooter spectra. A standard extinction curve⁴ (Weingartner & Draine 2001) was used to determine the wavelength-dependent extinction A_λ appropriate for a visual extinction A_V = 12 mag. The extinction-corrected flux, $f_{\nu }^{\mathrm{corr}}=f_{\nu }^{\mathrm{obs}}{10}^{0.4A_{\lambda }}$, is also shown in Table 2. The corrected flux is compared directly with the model. To reproduce our observed K-band spectrum we construct a two-zone model as an oversimplification of a dynamical atmosphere. In particular, models with T ≈ 6750 K, log n_H = 13.2 cm⁻³, and a ratio of thickness to radius ρ ≈ 0.1 yield peak flux densities in the NIR that are consistent with the observed Na I and Mg II emission, for an angular diameter θ ≈ 1.2 mas. The peak emission in the Brγ profile can be matched with a somewhat more extended region, θ ≈ 2.0 mas, at somewhat lower temperature, T = 5000–5500 K, and density, log n_H ≈ 11.6–12.8 cm⁻³. The optical depth in Brγ in this more extended, cooler region is high enough to account for the strong absorption covering the smaller continuum-emitting region, as seen in the P Cygni-type line profile (Figs. 2 and 16). Applying expanding wind models with velocity gradients to explain fully the exact shape of the Brγ emission line is beyond the scope of the paper.

Fig. 16. Representative model spectrum in K band (red) overlaid on the X-shooter data (black). The model is based on a static spherical shell in strict LTE. The residual spectrum is shown in green.

Both zones of line emission have low fractional ionisation in LTE, n(e)/n_H ≈ 10⁻⁸. Again, we point out that LTE conditions are sufficient to explain the observed emission lines of both neutrals (H I and Na I) and ions (Mg II) without any need for some kind of photo-ionised atmosphere (see Fig. 16). We note that we did not attempt to fit the observed data, and our models point to the range of conditions that can account for the fluxes and angular sizes of the observed IR emission lines and associated continuum. It also follows that the relative angular sizes of K continuum, Na I, and Brγ are likely the result of the fact that the optical depth for each varies with radius in the extended atmosphere. The different apparent sizes do not require distinct layers or shells. Brγ attains unit optical depth rather farther out in the atmosphere than Na I.

In the future, it would be worthwhile to develop a physically consistent model of the dynamical atmosphere and to explore additional constraints from the observations (for example fit the observed P Cygni-type line profile seen in Brγ emission).

The very similar reported mass-loss rates, luminosities, and temperatures of IRAS 17163 and IRC+10420, in combination with the presence of a more extended Brγ, with the coexistence of a more compact Na I emission in both cases, make this approach promising in explaining the K-band line emission spectrum towards more objects of this class.

7.2. Mass loss

7.2.1. Mass-loss rates

The consequences of the new distance to the Fried Egg Nebula (smaller optical depth, smaller sizes of the dusty shells) also impacts the derived mass-loss rates. In particular, the model results in a mass-loss rate of 6 × 10⁻⁷ M_⊙ yr⁻¹ for the hot inner shell, 9 × 10⁻⁴ M_⊙ yr⁻¹ for the intermediate shell, and 5 × 10⁻⁵ M_⊙ yr⁻¹ for the outer shell. The mass-loss rate in the outer shell is about two orders of magnitude lower than what was reported in Lagadec et al. (2011).

Large mass-loss rates, higher than 10⁻⁵ M_⊙ yr⁻¹, can be expected for this class of YHGs and have been previously measured towards other objects to be ∼10⁻⁴ M_⊙ yr⁻¹ (IRC+10420; Oudmaijer et al. 1996), which is comparable (higher by only a factor of 2) to what we find towards the Fried Egg Nebula. In some cases, the mass-loss rates are even as high as ∼5 × 10⁻² M_⊙ yr⁻¹ (ρ Cas Lobel et al. 2005). A revision of more YHGs, where the new Gaia parallaxes are taken into account, is necessary and it may change our current view on mass-loss rates. The mass-loss rate we report for the outer shell is of the same order of magnitude as that reported in Wallström et al. (2017) for the observed parsec scale shell ∼0.8 × 10⁻⁴ M_⊙ yr⁻¹. However, these values apply when adopting an expansion velocity of 30 km s⁻¹ (Hutsemékers et al. 2013) as opposed to 100 km s⁻¹ in Wallström et al. (2017), which differs by a factor of ∼3, thereby influencing the resulting mass-loss rate by the same factor. Wallström et al. (2017) estimate the expansion velocity using the CO 2–1 spectrum obtained with ALMA, and in particular by modelling the broad “plateau” emission attributing the emission to the stellar wind. Given that we have no information of the situation in the inner shells, being unresolved in Wallström et al. (2017), we also considered the more moderate expansion velocity.

As we already mentioned, the mass-loss rates vary among the three distinct shells. In particular, we find that the intermediate shell has the largest mass-loss rate, followed by the outer shell, while the inner hot shell has the lowest mass-loss rate, which is in fact lower than what it is measured towards other YHGs by about two orders of magnitude (Oudmaijer et al. 1996). The mass-loss rates for the individual shells differ by orders of magnitudes is probably a real effect and cannot be explained by the uncertainties introduced by the adopted parameters in our models alone (e.g. gas-to-dust ratio, and dust composition). Shenoy et al. (2016) report distinct mass-loss episodes towards a sample of hypergiants, including IRC+10420, at timescales of thousands of years. In particular, IRC+10420 also shows a decline in mass-loss rate up to an order of magnitude, while simulations predict an increase in mass-loss-rate as the object comes closer to exploding as a SN (Moriya et al. 2014). The physical mechanism responsible for such “eruptive” mass-loss events has been traditionally a topic under debate.

7.2.2. Mass-loss mechanisms

To explain the presence of shells, we discuss a pulsational driven mass-loss mechanism and a line driven mass-loss mechanism. In line with previous literature we begin with the classical interpretation of the pulsation hypothesis as the mechanism behind eruptive mass-loss events. We then introduce the bi-stability mechanism, which occurs in line-driven winds. Although is well known for hotter objects, this mechanism has thus far not been considered in YHG studies. Yet as we discuss, the bi-stability mechanism has the potential to explain the observed mass loss towards those objects as well.

Pulsational-driven mass loss. Yellow hypergiants (YHGs) are known to drive significant and unique quasi-period photospheric pulsations (e.g. Lobel et al. 1994; de Jager 1998), which under favourable conditions (T_gas ∼ 8000 K; Lobel 2001) and on timescales of decades can develop into an eruptive episode. The eruptions are characterised by twice the duration of an ordinary pulsation and by an increased mass-loss rate (e.g. ρ Cas, HR 8752; van Genderen et al. 2019; Lobel et al. 2003). Once hydrogen recombination occurs synchronously with the pulsation decompression (i.e. strong convective motions), this results in an outburst event and a subsequent drop of the temperature (van Genderen et al. 2019). Lebertre et al. (1989) showed a clear variability in the JHK data of IRAS 17163 in a timescale of 3 yr (up to 0$\stackrel{\text{m}}{.}$5 in J) but it was not possible to trace a periodic trend in those data. Because IRAS 17163 has not been followed up photometrically, it is not possible to identify the observed variability with characteristic patterns of a pulsation mechanism.

Adopting the stellar parameters (Sect. 4), the expected pulsation timescale is of order of months (see Eq. (4) in Oudmaijer & Bakker 1994), while in Sect. 7.3 we find timescales of decades for the shells towards the Fried Egg Nebula. We also note that the ejecta caused by the pulsation mechanism in the known objects are optically thick, resulting in a pseudo-photosphere which mimics a star much cooler than an F type. In the case of the Fried Egg Nebula, the shells are rather optically thin, and as we discuss in Sect. 7.1 a pseudo-photosphere is not required to explain the observed emission spectrum. Combining all the information above, we conclude that, if acting, the instabilities leading to the eruptions are less severe in the case of IRAS 17163. We therefore explore another mechanism other than classical pulsations as follows.

Line-driven mass loss. An alternative explanation would involve interacting winds that are driven by radiation pressure in spectral lines. Any change in radiative acceleration has an effect on both the wind velocity and mass-loss rates. IRAS 17163 is located in the white wall and yellow void area in the HRD. In particular Nieuwenhuijzen et al. (2012) present two regions: one at log T_eff ∼ 3.9–3.95 and one at log T_eff ∼ 4.05 − 4.1. In the first regime where IRAS 17163 is located, the ionisation of H might be responsible for increased dynamical atmospheric instabilities (Humphreys et al. 2002).

Although we would expect that the most abundant species (hydrogen, helium) are dominating the line force, Vink et al. (1999) showed that more complex ions have a higher contribution because of their very high number of transitions (thousands) in the relevant part of the spectrum where the stellar flux contribution is at its maximum. Iron (Fe) in particular, is a key wind driving element. Little drops in temperatures result in the recombination of Fe to lower ionisation stages, increasing the density of the wind, reducing its terminal velocity by about a factor of two, and increasing the mass-loss rate by several factors (Vink et al. 1999). This process is known as the bi-stability mechanism. Wind bi-stability is proposed to play a crucial role during shell ejections (Pauldrach & Puls 1990). In particular, these changes are expected at the first bi-stability jump occurring at T ∼ 21 000 K as a consequence of the recombination of Fe IV to III, and at the second bi-stability jump occurring at T ∼ 8800 K because of the recombination of Fe III to Fe II (Vink et al. 2000, 2001; Petrov et al. 2016). The evolutionary status and effective temperature of IRAS 17163, of order 8500 K, along with the presence of the surrounding multiple, rather optically thin, shells characterised by various mass-loss rates, justify the exploration of the second bi-stability jump.

The second bi-stability jump predicts escape and terminal wind velocities that are comparable to the outflow velocity of 30–100 km s⁻¹ observed towards IRAS 17163 (Lamers et al. 1995). When objects are unstable because of a proximity to the Eddington limit and approach the opacity bumps, they can enter fast and slow wind transition cycles, resulting in shells (e.g. Gräfener et al. 2012). If bi-stability is the mass-loss driving mechanism, then the number of the observed shells indicates that IRAS 17163 has undergone four such cycles, including that by Hutsemékers et al. (2013).

7.2.3. Mass-loss geometry

We show that IRAS 17163 is surrounded by spherical symmetric shells at ∼0$\stackrel{″}{.}$3–3″ scales (Sect. 6; see also, Lagadec et al. 2011). A similar morphology has been seen before in another object of its class, HD 179821 (Ferguson & Ueta 2010). This morphology has been traditionally explained as mass loss triggered by a large-scale pulsation mechanism in combination with radiation pressure on the dust.

As we presented in Fig. B.1, the channel maps (GRAVITY) of the inner regions (milli-arcsecond scales) of IRAS 17163 reveal morphological asymmetries at different velocities. This could be the result of outflowing material in different velocities and orientations, indicating non-uniform mass-loss events at small scales, which may all together lead to a more spherical and symmetric emission.

At larger scales (∼25″), Wallström et al. (2017) present kinematic evidence in CO emission of asymmetric mass loss and clumpiness connected with the outer shell reported in Hutsemékers et al. (2013). Asymmetric mass loss has been previously observed towards other YHGs, with the observed asymmetries attributed to the surface activity of a convective nature (e.g. IRC+10420; Tiffany et al. 2010).

IRAS 17163 was also imaged with the Hubble Space Telescope (HST), but with no sign of a dusty circumstellar shell (Siódmiak et al. 2008). The situation was similar for other YHGs (ρ Cas, HR 8752; Schuster et al. 2006). The lack of nebulosity in HST images has been previously connected to low mass-loss rates (e.g. μ Cep; de Wit et al. 2008), the possible faint nature of the shells at the observed wavelengths or the evolutionary status of the objects (Schuster et al. 2006).

7.3. Kinematic timescales

The shells are about three times closer to the central star compared to previous studies, and therefore the kinematic timescales, t_kin, of each shell have become shorter accordingly. In particular, the timescales of the intermediate and outer shell are found to be 100 yr and 120 yr, respectively. As a consequence, IRAS 17163 has likely passed to its YHG phase more recently than previously thought. The inner hot shell, which we report on and model in this work for the first time, is estimated to be a very recent mass-loss event with a maximum timescale of only ∼30 yr. It is striking that ∼30 yr ago Lebertre et al. (1989) presented a clear variability in the JHK data (up to 0$\stackrel{\text{m}}{.}$5 in J) towards the object. Given our new findings, this variability can very well be the result of an eruption that led to the formation of the inner third shell. It is worth mentioning that at the updated distance of 1.2 kpc the timescale of the cold, more extended, outer shell presented in Hutsemékers et al. (2013) is roughly reduced from 16 000 yr to 4000 yr. Our reported values are most likely an upper limit as they were derived for the lower expansion velocity of 30 km s⁻¹. For example, if we used 100 km s⁻¹ as done by Wallström et al. (2017), the values would be shorter by a factor of 3. With kinematic timescales of 100 yr or shorter, the star is expected to have evolved notably over the last century, which could be traced to changes also in the visual magnitude or the infrared fluxes. For future reference, it would be worth investigating if there are any archival⁵, historical photometric measurements to trace the long-term variability of IRAS 17163.

We conclude that the Fried Egg Nebula is a key object which is characterised by three distinct mass-loss events with varied mass-loss rates and maximum timescales from 30 yr up to 120 yr, with the most recent event to be the least powerful in terms of mass-loss compared to the previous two events. Given the observational restrictions, IRAS 17163 is among the very few cases where multiple mass ejections have been reported at such sort timescales (see also ρ Cas; Beardsley 1961; Shenoy et al. 2016). Such objects can shed light on fundamental questions on massive star evolution, such as how many repeating eruptions can be expected to occur for a particular object, what is the exact cycle of the eruptions, and what are the dependencies between the mass-loss rates and episodes, and the local conditions (e.g. mass, effective temperature, clumpiness, and metallicity).

8. Summary

This paper presents new optical and NIR spectroscopy (X-shooter) and high spatial resolution (down to mas) GRAVITY and AMBER observations of the YHG IRAS 17163. In particular we reported on the spectral type of IRAS 17163, the updated distance using the new Gaia measurements, and the geometry of the continuum emission and that of Brγ and Na I line emission close to the object, image reconstruction; we also presented modelling of published photometry by applying the radiative transfer code 2-Dust.

Evolved stars, with their high luminosity and extended atmospheres are great targets for optical and NIR interferometry, but yet their imaging is challenging; only few objects known to date meet the criteria for a successful image reconstruction at these wavelengths (e.g. η Car; Weigelt et al. 2016; Sanchez-Bermudez et al. 2018, Gravity Collaboration). We present the first image reconstruction of IRAS 17163 at short wavelengths (∼2 μm) and provide new insights into its close surroundings but also into the characteristics of its surrounding shells and mass-loss episodes.

The key results of our study are presented as follows:

• Based on classic spectral typing (between A3 and A6) and the analogy with IRC+10420, we infer the photospheric temperature of IRAS 17163 to be ∼8500 K.

• The continuum emission at 2 μm stems from the star and is not associated with hot dust.

• Our models reveal the existence of a third inner hot shell for the first time. This is a very recent mass-loss event with a timescale of a maximum of 30 yr.

• The three distinct shells trace three mass-loss episodes and are all characterised by variability in the mass-loss rate; the lowest mass-loss rate is that of the inner shell.

• The model-independent image reconstruction at the selected wavelengths (GRAVITY; Brγ, Na I, and continuum) shows the more extended and rather asymmetric nature of the Brγ emission (E–W direction; for example asymmetric wind origin) compared to Na I and continuum. In particular, Brγ shows also a northern and southern component, while Na I appears to follow the distribution of the photosphere. A similar situation has previously been observed by applying geometrical models towards another object of the same class (IRC+10420). Our spectropolarimetric analysis betrays asymmetries of order stellar radii along an EW position angle, which is in agreement with the proposed structure.

• The Brγ emission shows a P-Cygni profile and a drop in its peak intensity up to ∼40% in a four month period, which is indicative of a variable wind from the central YHG with variable gas mass-loss rates.

• Our two-zone LTE model can reproduce the K-band spectrum, and in particular the Na I doublet in emission, without the need of a pseudo-photosphere.

• We introduce the second bi-stability jump to explain the shells constituting the Fried Egg Nebula; this is the first approach to explain mass loss via this mechanism towards a YHG.

To summarise, our models show a newly discovered hot inner shell surrounding IRAS 17163. Follow-up mid-IR, high-angular resolution observations, for example, with the Multi AperTure mid-Infrared SpectroScopic Experiment (MATISSE/VLTI), will allow us to spatially resolve and image this shell, which will shed light on its morphology; therefore more information on the most recent mass-loss episode will be retrieved. We also show that a pseudo-photosphere is not required to explain the observed Na I doublet in emission. Instead, a spherical shell model in strict LTE can reproduce the K-band line emission spectrum of IRAS 17163. Given the similarity of the K-band spectrum among YHGs, such a model should be tested towards all known objects of this class.

Lastly, to interpret the observed spherical shells, we explored both pulsational-driven mass loss and line-driven mass-loss mechanisms. We conclude that with only one object and the observations in hand, we cannot discriminate between the two mechanisms. The bi-stability mechanism with its second bi-stability jump provides a promising and unified mechanism to explain mass-loss processes towards both hot (blue S Dor variables; first bi-stability jump) and cooler (YHGs; second bi-stability jump) stages of stellar evolution. A future study in which the pros and cons of both the pulsation and bi-stability mechanisms are thoroughly discussed towards a sample of all known YHGs will hopefully allow us to properly distinguish between the two theories on the variable mass loss at those stages.

Acknowledgments

EK is funded by the STFC (ST/P00041X/1). We thank Sundar Srinivasan, Pierre Kervella and Werner Salomons for their help when preparing the paper and useful discussions. We also thank both referees for a careful reading of the manuscript and for providing useful comments and suggestions that improved the paper. This research has made use of the AMBER data reduction package of the Jean-Marie Mariotti Center (http://www.jmmc.fr/amberdrs). Based on observations collected at the European Southern Observatory under ESO programme(s) 099.D-141 (X-shooter), 089.D-0576 (AMBER), 099.D-0664 (GRAVITY). This research has made use of the SIMBAD data base, operated at CDS, Strasbourg, France. This research has made use of the Jean-Marie Mariotti Center OImaging service (http://www.jmmc.fr/oimaging) part of the European Commission’s FP7 Capacities programme (Grant Agreement Number 312430). The GRAVITY data reduction presented in this work was undertaken using ARC3, part of the High Performance Computing facilities at the University of Leeds, UK.

References

Appendix A: Observations

Fig. A.1. As Fig. 7, but for the small configuration A0-B2-C1-D0. The drop in visibilities at ∼2.139 μm is based on a single data point and not connected to an emission line, therefore we consider it an artefact. The closure phase does not show changes in this configuration.

Fig. A.2. As Fig. 7, but for the medium configuration K0-G2-D0-J3. The changes in differential phases from positive to negative values correspond to similar baselines but different PA by > 30°, and are indicative of offsets of the photo centre projected in a specific baseline.

Fig. A.3. As Fig. 7, but for the medium configuration A0-G2-D0-J3.

Fig. A.4. As Fig. 7, but for the large configuration A0-G1-J2-J3. This is the only case where Na I also shows a change in the closure phase.

Appendix B: Image reconstruction

B.1. Channel maps

For the channel reconstruction maps of Brγ we followed the same process as in Sect. 5.2 but for the velocity channels at −42 (absorption component), 25, 66, 108, 122, and 150 km s⁻¹ (Fig. B.1). The channel maps reveal differences in morphology. In particular, at velocities < 100 km s⁻¹ the emission shows two distinct emitting clumps in W–E, while for velocities > 100 km s⁻¹ there is a change of 90° in orientation. This could indicate an outflowing material in different velocities and orientations, which may all together lead to a more spherical and symmetric emission. Higher spectral and spatial resolution observations are required to confirm the observed differences among the different velocities and will lead to a proper investigation of the substructures.

Fig. B.1. Reconstructed channel velocity maps of the Brγ emission towards IRAS 17163. The corresponding velocities are [−42, 25, 66, 108, 122, 150] km s−1 (LSR).

B.2. Combined dataset

In addition to Sect. 5.2 we combined the AMBER and GRAVITY datasets and followed the same process for the Brγ and continuum. Brγ shows time variations and its morphology is expected to be affected. Therefore although the reconstructed image of the combined dataset is a result of a more complete uv-coverage compared to the GRAVITY dataset alone, the image reconstruction is less robust due to the observed flux and size variations of the Brγ emission at the different epochs (see also Fig. 12). Given the above, the reconstructed image of the combined dataset (Fig. B.2) naturally differs from that of the GRAVITY dataset (Fig. 13).

Fig. B.2. Image reconstruction of the Brγ emission for the combined dataset (AMBER+GRAVITY) and for three different software packages, MIRA, WISARD and BSMEM, from left to right. The image is algorithm dependent, but they all show a similar global orientation of the emission. The resulting image using BSMEM shows an elongated structure around RA offset, Dec offset = (−5, 5), (−5, 5), which is perpendicular to that of the other two software packages. The observed difference is due to the different regularisers and optimisation process available for each algorithm.

A thorough comparison with other software packages was crucial to give us confidence that the image shows real structures and not artefacts. For this purpose in addition to MIRA we used BSMEM (Baron & Young 2008) and WISARD (Meimon et al. 2004, 2008). We present the resulted images of the Brγ spatial distribution for the combined dataset in Fig. B.2. The reconstructed images show similarities regarding the asymmetric and extended nature of the Brγ. In particular all algorithms reveal features in the NW–SE orientation. The resulted image using BSMEM shows an elongated central emission perpendicular to that of the other two software packages. The observed differences that can be seen in the detailed, pixel sized, structures in the combined dataset, are due to the different ways the algorithms treat the interferometric observables, such as the different available regularisers and optimisation process (i.e. gradient vs. global) per algorithm. Therefore, the observed differences among the algorithms are software-related artefacts. We conclude that the observed elongation in NW–SE orientation is real and indicative of an asymmetric recent mass-loss, but how thin or thick this emission is, is something we cannot constrain with the current dataset.

All Tables

All Figures

	Fig. 1. Coverage of uv-plane of VLTI/GRAVITY and AMBER observations of IRAS 17163.
In the text

Fig. 2. Normalised spectra of IRAS 17163 (GRAVITY and X-shooter) around the 2 μm continuum level. The presented wavelength range contains all the observed Brγ, Na I, and Mg II emission lines. The dotted line represents the normalised telluric absorption spectrum applied to correct the spectrum of the object. The observed difference in line-to-continuum ratio is due to the higher spectral resolution of the X-shooter.

In the text

	Fig. 3. Spectra of the Brγ (left) and Na I doublet (right), normalised to the continuum level, for the different observing dates and configurations in different colours. The dataset shows a gradual drop in the flux of the Brγ emission (up to ∼67%) over the observing period of four months, while Na I emission does not vary significantly.
In the text

	Fig. 4. Line profiles of the observed Brγ, Na I, and Mg II emission towards IRAS 17163. The dotted line represents the telluric absorption spectrum applied to correct the spectrum of the object. The vertical line represents the systemic velocity as reported by Wallström et al. (2015) using the optical Fe II emission.
In the text

	Fig. 5. Polarisation (bottom panel) and polarisation angle (top panel) observed towards the Fried Egg Nebula. The black dots represent the broadband values, computed in steps of 10 nm, or 100 Å, from the spectropolarimetric data. The red triangles represent the Lebertre et al. (1989) data. The dashed line is the best-fit Serkowki law for interstellar polarisation.
In the text

	Fig. 6. Hα spectropolarimetry of IRAS 17163. The data are presented as a triplot of the Q, U, and intensity Stokes vectors (from top to bottom). The data are rebinned to a constant error in polarisation of 0.16%. The excursion across the Q vector alone indicates an asymmetric scattering geometry at a projected angle of 90° across the sky.
In the text

Fig. 7. Representative example of relative flux, visibility, differential, and closure phase as a function of wavelength around the Brγ and Na I doublet emission towards IRAS 17163 using GRAVITY on the large configuration. The 6 baselines and 4 triplets of A0-G1-J2-K0 configuration are plotted with different colours. Both Brγ emission and Na I show a smaller visibility than the continuum while Brγ shows a larger drop. Mg II does not show changes in visibility or phases. A differential phase change is also detected towards Brγ but not towards Na I.

In the text

	Fig. 8. Representative examples of relative flux, visibility, differential, and closure phase as a function of wavelength for the AMBER D0-H0-G1 and D0-I1-H0 configurations around the Brγ emission towards IRAS 17163.
In the text

	Fig. 9. Spectra of IRAS 17163 with a number of spectral standards taken from the UVES spectral archive (Bagnulo et al. 2003). The spectra are normalised to a continuum value of 1 and offset, for clarity, in multiples of 0.9.
In the text

	Fig. 10. Overlay of part of the blue spectrum of IRAS 17163 (black line) and IRC+10420 (red line). There is a tantalising similarity between the two objects, as previously also noticed by Wallström et al. (2015) for the optical spectrum above 5000 Å, which is dominated by emission lines.
In the text

	Fig. 11. Dereddened IRAS 17163 photometry with total extinction AV = 12.0 to fit a Teff = 8500 K Bosz stellar model.
In the text

Fig. 12. Visibilities at the peak of the Brγ, Na I emission lines and continuum emission towards IRAS 17163 overplotted with simple geometrical Gaussian models. Na I and continuum emission follow very similar distribution. Na I shows a smaller drop in visibilities compared to Brγ, indicating that Na I stems from a smaller emitting area compared to the Brγ. The AMBER dataset was rebinned to match the spectral resolution of GRAVITY. A period of 5 yr indicates a shrinkage of the Brγ emitting region to half its size.

In the text

	Fig. 13. Image reconstruction of the continuum emission (greyscale) overplotted with Brγ in red contours at [0.1, 0.3, 0.6, 0.75, 0.95] × peak and Na I in blue contours at [0.2, 0.5, 0.8, 0.95] × peak towards IRAS 17163. Brγ covers a more extended region compared to the continuum and Na I, which appear to follow very similar distributions.
In the text

Fig. 14. Left: VISIR/VLT mid-IR monochromatic image at 8.59 μm (Lagadec et al. 2011). Right: radial profiles at 8.59 μm, 11.85 μm, and 12.81 μm are extracted and overplotted with the modelled profiles as described in Sect. 6.2. The size and position of each shell is indicated in light blue. The point spread function (PSF) at 11.85 μm is also overplotted for reference, revealing the hot inner shell. The PSF becomes broader at longer wavelengths.

In the text

Fig. 15. Spectral energy distribution fit to the dereddened photometry gathered from literature (symbols). The dashed line shows the 5 × 105 L⊙, 8500 K black-body (BB) input used for the model. The dotted line is the stellar and hot inner shell emission. The dot-dashed line shows the stellar and cool inner shell emission. The large spacing dashed line shows the outer shell emission contribution. We note that our models did not include the outermost cooler dust shell (18–40″) reported in Hutsemékers et al. (2013), and therefore do not reproduce the far-infrared photometry.

In the text

	Fig. 16. Representative model spectrum in K band (red) overlaid on the X-shooter data (black). The model is based on a static spherical shell in strict LTE. The residual spectrum is shown in green.
In the text

	Fig. A.1. As Fig. 7, but for the small configuration A0-B2-C1-D0. The drop in visibilities at ∼2.139 μm is based on a single data point and not connected to an emission line, therefore we consider it an artefact. The closure phase does not show changes in this configuration.
In the text

	Fig. A.2. As Fig. 7, but for the medium configuration K0-G2-D0-J3. The changes in differential phases from positive to negative values correspond to similar baselines but different PA by > 30°, and are indicative of offsets of the photo centre projected in a specific baseline.
In the text

	Fig. A.3. As Fig. 7, but for the medium configuration A0-G2-D0-J3.
In the text

	Fig. A.4. As Fig. 7, but for the large configuration A0-G1-J2-J3. This is the only case where Na I also shows a change in the closure phase.
In the text

	Fig. B.1. Reconstructed channel velocity maps of the Brγ emission towards IRAS 17163. The corresponding velocities are [−42, 25, 66, 108, 122, 150] km s−1 (LSR).
In the text

Fig. B.2. Image reconstruction of the Brγ emission for the combined dataset (AMBER+GRAVITY) and for three different software packages, MIRA, WISARD and BSMEM, from left to right. The image is algorithm dependent, but they all show a similar global orientation of the emission. The resulting image using BSMEM shows an elongated structure around RA offset, Dec offset = (−5, 5), (−5, 5), which is perpendicular to that of the other two software packages. The observed difference is due to the different regularisers and optimisation process available for each algorithm.

In the text