© ESO, 2016

1. Introduction

The luminous blue variable (LBV) stage represents a crucial and relatively short phase in the massive star evolution. The limited proportion of stars in this stage indicates that this evolutionary phase lasts only for a short time (~ 10⁴ yrs). Although a clear definition of LBVs is missing, these objects exhibit a high luminosity (≥ 10^5.5 L_⊙) and a photometric variability, called S Dor variability, with different timescales and different amplitudes (~ 0.1 mag). Giant eruptions in which their brightness increases significantly (>2 mag) for only a few years can also occur (see e.g., Humphreys & Davidson 1994; Weis 2011, for further details). This variability is responsible for the change of their spectrum. In a hot phase, the spectrum of an LBV is similar to the spectrum of a blue supergiant (BSG) whilst in a cool phase, it looks like that of an A or an F star. LBVs also exhibit episodes of high mass-loss rate (~ 10^-5−10^-4M_⊙ yr^-1). The material lost during extreme mass-loss episodes frequently produces a circumstellar nebula constituted of stellar ejecta. These nebulae indeed generally display a nitrogen overabundance, inferring a stellar origin. They can be created through giant eruptions and then shaped by wind-wind interactions. The stellar wind can become slower and denser, proceeds into a cavity carved by the massive star wind and partially fills it up with ejected material (van Marle et al. 2007).

These nebulae have sizes that are generally estimated to about 1−2 pc. They exhibit expansion velocities that are generally between ten and a few 100 km s^-1 but can reach several 1000 km s^-1 (such as the outer ejecta of η Car, Weis 2012, and the references therein). These nebulae can also have different morphologies. Weis (2011) yielded that about 50% of them are bipolar, 40% spherical and only 10% are of irregular shape. The bipolar nebulae can display either an hourglass shape like SWB1, MN18 or Sher 25, or bipolar attached caps like WRAY 15-751. Even though the formation of LBV nebulae is not yet fully understood, possible causes of bipolarity have been advanced, such as the rotation of the star, a density gradient in the stellar wind, or different mass-loss episodes in which the wind changes from equatorial to polar directions during the passage of the bistability jump.

The circumstellar environment thus reveals the mass-loss history of the central star in the previous evolutionary stages. In some cases, more than one ejection can occur. Jiménez-Esteban et al. (2010) reported multiple shells in the nebula surrounding the LBV G79.29+0.46, and Vamvatira-Nakou et al. (2013) showed the presence of a second nebula around WRAY 15-751.

In the present paper, we focus on the candidate LBV HD 168625 (or IRAS 18184–1623). This object is a very luminous hypergiant B star (M_V = −8.6, van Genderen et al. 1992). Its spectral type seems to vary from B 2 (Popper & Seyfert 1940) to B 8 (Morgan et al. 1955) by going through B 5.6 (Chentsov & Luud 1989) even though no dramatic light-variation has been reported in the literature in the last 40 yrs (van Genderen et al. 1992 and Sterken et al. 1999). Besides this lack of variation, preventing the classification of HD 168625 as LBV, Chentsov & Luud (1989) noticed an infrared excess at λ> 5 μm that could be related to circumstellar dust. Several years later, Hutsemékers et al. (1994) provided clues of a high mass-loss history since they observed, for the first time, the nebula surrounding this star. They pointed out from a Hα image that this shell consists of two regions: an inner 10″ × 13″ elliptical ring detached from the star with a major axis oriented ESE-WNW and an outer horn-shaped (or ansa-shaped) region oriented perpendicularly to the inner region and suggesting a bipolar morphology. From a deeper picture of the Hα nebula, Nota et al. (1996) detected faint filaments in this bipolar structure that extend its size to 16″ × 21″. Pasquali et al. (2002) estimated the size of the nebula between 15.5″ × 23.5″ from a 4 μm image and 31″ × 35.5″ from an 11 μm image, both obtained with the Infrared Space Observatory (ISO). In the 4 μm image, the nebula is relatively faint in comparison with the central star whilst in the 11 μm image, the star is fainter which allows to see two bright regions of dust emission, separated by 15″ and symmetrically located to the northwest and the southeast with respect to the central object. The torus-shaped dust emission was modeled by O’Hara et al. (2003) who found a total dust mass of 2.5 ± 0.1 × 10^-3 M_⊙. The question of the origin of a filamentary bipolar structure in the nebula around HD 168625 was risen up by Smith (2007). A first explanation resulting from the interaction between the faster wind of the star and slower ejections was given but the possibility of other scenarios involving either the evolution of one single star or interactions between two massive stars have also been advanced. In this context, Aldoretta et al. (2015) detected the presence of a companion through interferometric observations. This detection has been confirmed in a very recent study of adaptative optics images by Martayan et al. (2016) that reveals a wide-orbit companion with a projected separation of 1.15″. We note that this companion was previously detected by Pasquali et al. (2002), who suggested already that HD 168625 might be a wide binary system. The role of this companion could affect the shape of the nebula or the evolutionary process, but, to date, this remains conjecture. Furthermore, no hint exists that this companion is gravitationally bound to the central star.

The distance at which the star is located at was also heavily debated. Hutsemékers et al. (1994) and Nota et al. (1996) used for their respective analysis 2.2 kpc (the estimated distance of the neighboring Omega nebula M17) but the former reported that this distance could be underestimated and that HD 168625 is probably situated behind M17 and its visual neighbor LBV HD 168607. Robberto & Herbst (1998) determined a distance of 1.2 kpc for HD 168625. Pasquali et al. (2002) estimated from the Galactic kinematic rotation model of Brand & Blitz (1993) a distance of 2.8 kpc to HD 168625 that seems more realistic. We thus consider this last value in our analysis.

To constrain the evolutionary path of HD 168625 and possible also to determine the influence of its putative companion during its evolution, we investigate, in this paper, the difference of CNO abundances between the nebula and the star on the basis of Herschel/PACS data and of high-resolution optical spectroscopy, respectively. We also discuss the global morphology of the nebula through a large set of Hα and infrared images. The paper is organized as follows. In Sect. 2, the observations and the data reduction procedures are presented. In Sect. 3, we model the central star by using the CMFGEN atmosphere code (Hillier & Miller 1998) to determine its fundamental parameters and the surface abundances of the elements. In Sect. 4, we present the morphology of the nebula and we analyse its Herschel/PACS spectrum. We discuss the results and constrain the evolutionary stages of HD 168625 in Sect. 5 and finally we provide the conclusions in Sect. 6.

2. Observations and data reduction

2.1. Infrared observations

2.1.1. Imaging

The far-infrared observations were obtained, with the PACS photometer (Photodetector Array Camera Spectrometer, Poglitsch et al. 2010) onboard the Herschel spacecraft (Pilbratt et al. 2010), in the framework of the Mass-loss of Evolved Stars (MESS) guaranteed time key program (Groenewegen et al. 2011). The imaging data were taken on March 30, 2011 which corresponds to the observational day (OD) 686 of Herschel. The scan map mode was used. In this mode, the telescope scans at constant speed, corresponding to 20″/s in this case, along parallel lines to cover the required area of the sky. A total of eight scans were obtained (two orthogonal scan maps at 70 μm and at 100 μm, four scans obtained simultaneously at 160 μm), over four distinct observations of 157 s each, carrying the following observation numbers (obsID): 1342217763, 1342217764, 1342217765 and 1342217766. The Herschel Interactive Processing Environment (HIPE, Ott 2010) was used for the data reduction up to level 1. Subsequently, the Scanamorphos software (Roussel 2013) was used to further reduce and combine the data. The Herschel/PACS point spread function (PSF) full widths at half maximum (FWHMs) are 5.2″, 7.7″ and 12.0″ at 70 μm, at 100 μm and at 160 μm, respectively.

We also used other space-borne images from Spitzer Infrared Array Camera (IRAC; Fazio et al. 2004) and Multiband Imaging Photometer of Spitzer (MIPS; Rieke et al. 2004), from Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), as well as the 11−18 μm VLT/VISIR (VLT Mid-Infrared Imager Spectrometer) images studied by Umana et al. (2010), in order to better constrain the nebular environment.

2.1.2. Spectroscopy

The far-infrared spectrum of the nebula surrounding HD 168625 was taken on September 25, 2011 (OD 865), with the PACS integral-field spectrometer that covers the wavelength range from 55 μm to 220 μm in two channels that operate simultaneously in the blue 55−98 μm band (B2A and B2B) and the red 102−220 μm band (R1A–R1B). It provides simultaneous imaging of a 47″ × 47″ field-of-view, resolved in 5 × 5 square spatial pixels (called spaxels). The two obsIDs of these observations are 1342229740 and 1342229741. The resolving power of the spectrum is between 1000 and 3000 depending on the wavelength. HIPE was also used for the data reduction. The standard reduction steps were followed and in particular the subtraction of the background spectrum obtained through chopping and nodding. However, a special treatment was brought to the [C ii] 158 μm line. This line was indeed found in absorption in the output spectrum after standard reduction. It appeared that a significant contribution to this line comes from the background spectrum. We therefore estimated manually the emission only produced by the background and we removed it properly from the nebular spectrum. The errors on the flux are thus larger for this line than for the other spectral lines detected in the spectrum.

2.2. Visible observations

2.2.1. Imaging

We retrieved from the ESO archives the coronographic images taken in Hα and a continuum filter (λ_c = 6506 Å, FWHM = 35 Å) with the SUSI camera mounted on the New Technology Telescope (NTT) at La Silla on May 23, 1995. These data were presented by Nota et al. (1996). We followed the same steps for the data reduction and refer to these authors in this respect. The occulting bar has a width of 2″ and the pixel size of the image is 0.1″.

2.3. Spectroscopy

We retrieved high-quality optical spectra of HD 168625 available in the ESO archives. A first set of data was obtained with the échelle Fiber-Fed Extended Range Optical Spectrograph (FEROS) mounted successively on the 1.52 m (before 2003), then on the 2.2 m telescope at La Silla (Chile). This instrument has a resolving power of 48 000 and the detector was a 2k × 4k EEV CCD with a pixel size of 15 μm × 15 μm. Exposure times between 4 and 30 min were used to reach signal-to-noise ratios, measured on the 4800−4825 Å wavelength region, between 50 and 100, respectively. FEROS provides 39 spectral orders, covering the 3700−9200 Å wavelength domain. For the data reduction process, we used an improved version of the FEROS pipeline working under MIDAS environment (Sana et al. 2006). We combined the spectra obtained consecutively the same night to increase the signal-to-noise ratio.

A second set of data was composed of three UVES (Ultraviolet and Visible Spectrograph) spectra taken at ESO Very Large Telescope (UT2, Paranal, Chile). One spectrum was taken in the dichroic mode with the DIC 1 346+580 and the DIC 2 437+860 setups and an 0.8″ slit for all setups. The other spectra were taken with the RED 580 nm setup and an 0.8″ slit. The resolving power of this instrument is 80 000. The data reduction was performed with the standard reduction pipeline.

A third set of two spectra was acquired with the X-shooter spectrograph mounted on the ESO VLT at Paranal. The data were obtained in nodding mode. Each observation covers the UBV (3000−5500 Å), the VIS (5500−10 000 Å) and the NIR (10 000−25 000 Å) arms. The reached resolving power is about 6000, 10 000 and 8000 in the three arms, respectively. The data were also reduced with the standard pipeline.

The 16 spectra belonging to the entire dataset were normalized by fitting polynomial functions of degree 4−5 on carefully chosen continuum windows with a width of about 200 Å. The orders were then merged to form a single spectrum. The barycentric correction was applied when it was required. The full journal of observation is provided in Appendix A.

3. Modeling of the stellar spectrum

To determine the fundamental parameters of HD 168625, we use the non-LTE model atmosphere code CMFGEN. This code solves the radiative-transfer equation for a spherically symmetric wind in the co-moving frame under the constraints of radiative and statistical equilibrium. The hydrostatic density structure is computed from mass conservation and the velocity structure is constructed from a pseudo photosphere structure connected to a β-velocity law of the form v = v_∞(1−R/r)^β where v_∞ is the wind terminal velocity. Our final model includes the following chemical species: H i, He i-ii, C ii-iv, N i-v, O i-vii, Ca ii, Na i, Mg ii, Ne ii-iii, Ti ii, Al ii-iii, Cr ii-iii, Si ii-iv, S ii-v, P iv-v, Fe ii-vii, and Ni ii-iv with the solar composition of Grevesse et al. (2010) unless otherwise stated. CMFGEN uses the super-level approach to reduce the memory requirement. In the present analysis, we include 1969 super-levels for a total of 7902 transitions. Once the emerging spectrum is generated, we convolve it first by a rotational profile determined from the Fourier transform method (Simón-Díaz & Herrero 2007) and then we convolve the resulting spectrum by a Gaussian profile to include the atmospheric macroturbulent velocity (v_mac).

The observed spectra do not show any double-lined spectroscopic (SB2) signature that could be related to the putative companion. We measure, from our entire dataset (FEROS, UVES, and X-shooter spectra), a mean radial velocity (RV) of 19.7 km s^-1 by fitting Gaussian profiles on the line profiles and a RV shift of about 15.4 km s^-1, but without determining a clear period. This RV shift is computed on a timescale of about 1000 days but our dataset does not allow us to be more accurate so that it must be considered as an upper limit. As for the possibly related star Sher 25 (Taylor et al. 2014), the amplitude of the RV shift does not allow us to claim that HD 168625 is a binary system, but we cannot exclude this assumption. Another difference between the optical spectra comes from Hα profiles (see Fig. 1). Many evolved massive stars have a Hα profile that varies without being in a binary system because of variations of the wind parameters (Ṁ, v_∞, β and/or the clumping filling factor f) as a function of time. These variations do not affect the main indicators useful to determine the stellar parameters (T_eff, log g, and/or log (L/L_⊙)) nor the abundances. All in all, we therefore consider that the optical spectrum of HD 168625 is that of a single star. We thus determine the stellar and wind parameters from the spectrum taken on June 29, 2002. This spectrum has been chosen because of its high S/N ratio and the P-Cygni profile of the Hα line. The investigation of the spectral variability will be carried out in a future paper.

Fig. 1 Comparison between the FEROS spectrum taken on June 29, 2002 (black line), the UVES spectrum taken on July 10, 2001 (red line) and the X-shooter spectrum taken on June 2, 2014 (blue line). The He i 4471 and the Mg ii 4481 lines are shown as effective temperature diagnostic, the N ii lines between 4600 and 4625 Å as nitrogen content indicator, and the Hα line for the stellar wind variations.

To estimate the main parameters of HD 168625, we focus on the following diagnostics:

• Effective temperature: the estima-tion of T_eff is done from theHe i 4471/Mg ii 4481line ratio and from ionization balances. The optical spec-tra of HD 168625 allow us to only con-sider the Si iii/Si ii andFe iii/Fe ii ratios.

• Gravity: the wings of the Balmer lines are used to determine the log g. We use Hϵ, Hδ, Hγ and to a lesser extent Hβ, although Hβ and Hα are affected by the stellar wind, making them less reliable diagnostic lines for gravity.

• Stellar luminosity: the stellar luminosity is computed from the absolute magnitude M_V, through the following formula: $\begin{array}{ccc}{\mathit{M}}_{\mathit{V}}\mathrm{=}\mathit{V}\mathrm{-}\mathrm{(}{\mathit{R}}_{\mathit{V}}\mathrm{\ast }\mathit{E}\mathrm{(}\mathit{B}\mathrm{-}\mathit{V}\mathrm{\left)}\mathrm{\right)}\mathrm{-}\mathrm{\left(}\mathrm{\right(}\mathrm{5}\log \mathit{D}\mathrm{)}\mathrm{-}\mathrm{5}\mathrm{)}\mathit{,}& & \end{array}$where D is the distance to the star. We use V = 8.37 and (B−V) = 1.41 (Ducati 2002). The value of (B−V)₀ is taken from Schmidt-Kaler (1982) for a spectral type B 6 close to that of HD 168625. From these values, we obtain E(B−V) = 1.46 in agreement with the value given by Hutsemékers et al. (1994). Assuming R_V = 3.1, and a distance of 2.8 kpc (Pasquali et al. 2002), we estimate M_V to be equal to − 8.39. If we apply the bolometric correction provided by Schmidt-Kaler (1982) for a star with the same T_eff, we obtain a stellar luminosity of log (L/L_⊙) = 5.69. With the bolometric correction provided by Martins & Plez (2006), this becomes log (L/L_⊙) = 5.46. For the rest of this work, we consider the average of both values, i.e., log (L/L_⊙) = 5.58 ± 0.11.

• Terminal velocity: given the spectral range of our spectra, the main diagnostic line to determine v_∞ is Hα. According to Martins (2011), if the line is in pure emission, v_∞ can be measured from the line width. If the line has a P-Cygni profile, v_∞ can be estimated from the blueward extent of the absorption part.

• Mass-loss rate: as for v_∞, the Hα line is the only available diagnostic line to measure Ṁ. We focus on the emission part of the line to determine this quantity. We emphasize however that this value is directly bound to other parameters such as v_∞, or the clumping law.

Once the fundamental parameters are constrained, we build a grid of models changing, one by one, the He and the CNO abundances. We perform a χ² analysis to estimate these abundances, following a method described in detail by Martins et al. (2015). To briefly summarize this method, the χ² is calculated for each chemical element by varying its abundance value. The resulting χ² function is then renormalized so that the minimum value is equal to 1.0. We use the same lines than those quoted by these authors. The 1σ error is set to the abundances where χ² = 2.0 (see Martins et al. 2015, for further details). The fundamental parameters resulting from our CMFGEN analysis are listed in the upper panel of Table 1 and the best-fit model is displayed in Fig. 2.

Fig. 2 Best-fit CMFGEN model for HD 168625 (red line) compared to FEROS spectrum taken on June 29, 2002 (black line). The lines that are not modeled are diffuse interstellar bands (DIBs). The Hα line has been rescaled for clarity.

The star shows an enrichment of the helium and nitrogen abundances and a depletion in its carbon and oxygen contents. This confirms its evolved nature. Its projected rotational velocity is 60 km s^-1 which is rather slow. Its temperature range agrees with previous investigations found in the literature (see, e.g., García-Lario et al. 2001) whilst its log g is smaller. A comparison between the FEROS spectrum and the UVES and X-shooter spectra is done in Fig. 1. No significant effective temperature change can be seen. Furthermore, the nitrogen lines are mostly the same as a function of time. The small variations do not impact on the value of the N abundance, and are taken into account in the 1σ error provided on this value. We emphasize that the line-profile variations are not strong enough to characterize the variability associated by an S Dor cycle.

4. Analysis of the nebula

4.1. Morphology of the nebula

The nebula around HD 168625 has been observed from the Hα to 160 μm wavelength range. All these images are displayed in Fig. 3. The coronographic Hα image, taken from the data analyzed by Nota et al. (1996) shows the ionized nebula. This nebula is composed of two distinct structures:

• a bright inner elliptical torus with a major axis orientedESE-WNW and with a size of12″ × 16″, corresponding to 0.16 pc × 0.22 pc at 2.8 kpc;

• a perpendicular bipolar/biconal structure (with a major axis oriented NNE-SSW with a northern lobe clearly defined whilst the southern one is weaker. The size of this entire structure is estimated to 16″ × 21″, i.e., 0.22 pc × 0.29 pc at 2.8 kpc.

The inner torus is tilted of about 40° with respect to the plane of the sky. The annular form composing the northern lobe of the bipolar/biconal structure has a size of about 8.8″ × 16.5″, i.e., 0.12 pc × 0.22 pc at 2.8 kpc, which therefore infers to the bipolar/biconal structure an inclination of about 60° with respect to the plane of the sky. From the Hα image, we also see a small cavity between the star and the inner torus. The thickness of the shell is estimated to 1.1″ (i.e., 0.02 pc at 2.8 kpc). The inner shell is also clearly observed in the continuum image whilst the outer structure is not observed. As outlined by Hutsemékers et al. (1994), it suggests that the nebula reflects the stellar luminosity by scattering due to large dust grains or shows non-equilibrium dust emission. The ionized nebula was also observed in the radio domain, at 8.64 GHz, by Leitherer et al. (1995). The size of 20″ measured from the north to the south and of 18″ from the east to the west agrees with that estimated from the Hα image. In the radio image, the two strongest emissions are detected at 4.5″ to the west and at 4.5″ to the south with respect to the central star. These locations are marked by crosses in the Hα image of Fig. 3.

Another larger bipolar/biconal structure oriented NNE-SSW is detected in the 8 μm Spitzer image and discussed by Smith (2007). The northern lobe has a size of 34″ × 70″, i.e., 0.46 pc × 0.95 pc at 2.8 kpc. The southern one is barely visible but seems to have the same size as the northern one. The inclination of these lobes infers to the bipolar/biconal structure an angle of 60° from the plane of the sky.

Fig. 3 Top: from left to right: Hα image taken with SUSI on NTT (left), the 8 μm Spitzer image, and the 12 μm VLT/VISIR image. Middle: the MIPSGAL 24 μm image, the PACS 70 μm image and the PACS 100 μm image from left to right. Bottom: the PACS 160 μm images. For all the images, north is up and east is left. The white arrow indicates the location of the LBV HD 168607. In the Hα image of Nota et al. (1996), the crosses indicate the locations of the maxima of emission detected in the radio image of Leitherer et al. (1995). We emphasize that the 12 μm image is from Umana et al. (2010) and that the inner part of the 24 μm image is saturated.

Mid-infrared images have been taken by Robberto & Herbst (1998) and by Umana et al. (2010) to map the dust distribution in the nebula. The nebula at these wavelength is slightly larger than the ionized nebula. It is composed of an inner torus with a size of 14.5″ × 20″, i.e., 0.19 pc × 0.27 pc at 2.8 kpc, giving to it an inclination of about 40°, in agreement with the inner torus detected in the Hα image. We also see that the dust is highly condensed in the west-northwest and in the east-southeast part of the torus. In addition, we detect a more diffuse emission oriented perpendicularly to the dust torus, and probably related to the visible bipolar/biconal structure. We roughly estimate the size of the northern lobe to 9″ × 20″, i.e., 0.12 pc × 0.27 pc at 2.8 kpc. From these values, the inclination of about 60° agrees with that measured for the bipolar/biconal structure seen in the Hα and the 8 μm Spitzer images. According to Umana et al. (2010), there exists a trend in the dust emission to be more prominent in the southern part of the equatorial torus for wavelength smaller or equal to 13 μm whilst the dust emission comes preferably from the northern part of the dust torus at larger wavelengths. O’Hara et al. (2003) interpreted this trend as a result of variations of the grain size in the different regions of the nebula or of a gradient in the dust/gas ratio. Mid-IR images also allowed Skinner (1997) and Umana et al. (2010) to detect both polycyclic aromatic hydrocarbon (PAH) and silicate features in the dust shell. The latter authors also pointed out the presence of a photodissociation region (PDR) around the ionized nebula. Furthermore, an analysis by Blommaert et al. (2014) of the 69 μm line observed with Herschel shows the presence of pure forsterite grains containing no iron around HD 168625.

The 22 μm WISE (see Martayan et al. 2016) and the MIPSGAL 24 μm images (Fig. 3) should provide a relatively similar morphology of the nebula around HD 168625. The comparison of both images shows that the inner torus-like nebula is more extended in the WISE image than in the MIPSGAL one, because of a larger width of the PSF in the WISE data. Two transversal cuts have been performed in the ESE-WNW and in the NNW-SSE directions to determine the size of the nebula. From the better quality MIPSGAL image, we estimate it to 67″ × 85″, i.e., 0.91 pc × 1.15 pc at 2.8 kpc, much larger than what we determine until now. At these wavelengths, the perpendicular bipolar/biconal structure is not detected. The rings detected by Smith (2007) from the 8 μm image coincide with the border of the nebula observed at 24 μm.

The far-infrared PACS images outline the cold dust distribution in the nebula around HD 168625. These images confirm the morphology already observed from the mid-IR images. We see a torus-like nebula oriented ESE-WNW where the southeast and the northwest parts are the brightest ones. This could suggest a temperature and/or dust grain size gradient in the nebula around HD 168625. Its size at these wavelengths is estimated to 29″ × 36″, i.e., 0.37 pc × 0.49 pc at 2.8 kpc. The inclination is thus about 40° as all the structures detected in the ESE-WNW axis. An "arm"-like feature is detected in the northern part (and to a lesser extent in the southern part of the nebula) on the three far-infrared images. This structure seems related to the bipolar nebula oriented along the NNE-SSW axis. However, the filamentary loops detected at 8 μm are not visible in the PACS images. Once again, it does appear that the northern lobe is the brightest one in the bipolar/biconal structure. We also note a more diffuse emission around the inner dust torus smaller than the size of the nebula estimated from the 24 μm image but with a size of 45″ × 59″, i.e., 0.61 pc × 0.80 pc at 2.8 kpc. The physical relation of this emission as a part of the nebula related to the star is nevertheless questionable.

The three-color image (Fig. 4) shows that HD 168625 is embedded in a cavity in the interstellar medium probably formed by the star through previous mass-loss episodes. We also see that HD 168607 that is located at 1′ of HD 168625 does not show any associated nebula. No cavity seems to exist around this confirmed LBV from the PACS images whilst a clear cavity was observed by Hutsemékers et al. (1994) in their Hα image.

The global morphology of the nebula around HD 168625 is thus extremely complex. To summarize, it is composed of torus-like shells visible at different wavelengths and oriented in the ESE-WNW direction. These structures seems to be tilted of 40° with respect to the plane of the sky. In addition, a perpendicular bipolar/biconal structure, oriented in the NNW-SSE direction, is also detected in the Hα, 8 μm, and 13 μm images. This structure seems to have an inclination of 60° with respect to the plane of the sky. Furthermore, the northern lobe is brighter than the southern one.

Fig. 4 Three-color image (70 μm in blue, 100 μm in green and 160 μm in red) of the HD 168625 nebula. North is up and east is to the left.

4.2. The Infrared nebular spectrum

Figure 5 illustrates the footprint of the PACS spectral field-of-view on the image of the nebula at 70 μm. This field-of-view is composed of 25 (5 × 5) spaxels, each corresponding to a different area of the nebula. The nebula seems to be covered by the entire spectral field-of-view but to be sure that no flux is lost by the detector, we compare the flux of the spectral continuum to the photometric flux densities measured from the 70 μm, the 100 μm and the 160 μm images (for further details on the method, see Vamvatira-Nakou et al. 2013). For HD 168625, a good agreement exists and no correction factor must be applied. The analysis of the dust Spectral Energy Distribution including PACS data will be discussed in another paper although the results are relatively close to those provided by O’Hara et al. (2003).

The integrated spectrum of the nebula over the 25 spaxels exhibits the [O i] 63 μm, [N ii] 122 μm, [O i] 146 μm, [C ii] 158 μm, and [N ii] 205 μm lines (Fig. 6). Furthermore, this spectrum also contains the band of forsterite at 69 μm that was discussed by Blommaert et al. (2014). Given the spectral type of the central star, it is obvious that the spectrum of the surrounding nebula cannot display any higher ionization lines such as [N iii] 57 μm or [O iii] 88 μm. We emphasize that the [N ii] 205 μm line has a problematic calibration in PACS. This issue has already been mentioned by Vamvatira-Nakou et al. (2013, 2015). To summarize, the flux measured for this line is incorrect because of a light leak that superimposes a significant amount of light from 102.5 μm at that wavelength. The relative spectral response function (RSRF) used to reduce the data suffers from the same light leak. Consequently, when the RSRF is applied during the data reduction process, the signal at wavelengths larger or equal to 190 μm is divided by a too high number. Therefore, the flux of the [N ii] 205 μm line must be corrected by a factor of 4.2. An error of 25% is assumed for the corrected final flux of the [N ii] 205 μm line.

To measure the intensies of the lines from the whole nebula, we perform a Gaussian fit on the different line profiles of the integrated spectrum. The measurements are listed in Table 2. We also do the same measurements on all lines present in the different spaxels to estimate the flux distribution across the nebula. As we already mentioned, the errors on the flux measurements of the [C ii] 158 μm line are larger than for the other lines. The flux measurements on the different spaxels are given in Appendix B.

Fig. 5 Footprint of the PACS spectral field-of-view on the image at 70 μm of the nebula. The spaxels are represented in gray. North is up and east is left.

Fig. 6 PACS spectrum of the HD 168625 nebula, integrated over the 25 spaxels. The different bands are indicated with different colors. The continuum and the line at 205 μm are not calibrated, the inset is expressed in arbitrary units. The spectrum between 95 and 102 μmis affected by an instrumental artefact that makes the spectrum scientifically unusable in this wavelength domain.

4.3. Electron density and abundances

Nota et al. (1996) determined an upper limit to the electronic temperature T_e of 7000 K as well as the electron density distribution n_e across the nebula. The density measured at 3″ to the north shows a distribution between 200 and 1600 cm^-3 whilst n_e reaches between 400 and 3000 cm^-3 at 2″ to the south. These authors inferred a mean value of 1000 cm^-3 from all their measurements.

In the optical domain, the electron density is computed from the ratio [S ii] 6717 Å/[S ii] 6731 Å. Hutsemékers et al. (1994) estimated this ratio to about 0.9., giving an electron density of n_e = 630 ± 300 cm^-3 for the whole nebula around HD 168625. In the infrared waveband, the ratio between [N ii] 122 μm and [N ii] 205 μm is a good indicator to compute the electron density of the nebula. Considering the values provided in Table 2, we compute a ratio [N ii] 122/[N ii] 205 = 5.7 ± 2.9 for the whole nebula. By using the nebular/IRAF package, and by assuming a T_e of 7000 K, we determine an electron density of n_e = 350 ± 270 cm^-3. We consider this estimate of the electron density in the following calculations because the electron density is best determined when it is similar to the critical density of the lines whose ratio is used as diagnostic (Rubin et al. 1994).

In the literature, we find many different estimations of the Hα flux from the nebula around HD 168625. Hutsemékers et al. (1994) estimated a reddened flux of 4.6 × 10^-13 erg s^-1 cm^-2 and a dereddened flux of 1.2 × 10^-11 erg s^-1 cm^-2 given a E(B−V) = 1.46. Nota et al. (1996) provided values of 1.8 × 10^-12 erg s^-1 cm^-2 and of 3.7 × 10^-10 erg s^-1 cm^-2 for the reddened and the dereddened fluxes, respectively, by assuming E(B−V) = 1.86. Finally, Pasquali et al. (2002) estimated a reddened and a dereddened flux of 4.0 × 10^-12 erg s^-1 cm^-2 and 9.6 × 10^-10 erg s^-1 cm^-2, respectively. It is thus clear that the final value of the Hα flux depends on the extinction (different from one study to another) and on the contamination by dust scattering of the stellar continuum and Hα emission that is difficult to correct for. Therefore, to be more accurate on the Hα flux, we use the radio flux of 14 ± 2 mJy measured at 8.64 GHz by Leitherer et al. (1995) and we convert it into Hα flux using Eq. (1) of Vamvatira-Nakou et al. (2016). We obtain F₀(Hα) = 1.84 ± 0.26 × 10^-11 erg s^-1 cm^-2. We then apply the nebular/IRAF package to determine the N/H abundance ratio of the nebula. We assume a case-B recombination (Draine 2011), with T_e = 7000 K and with n_e = 350 cm^-3, and we derive a nitrogen abundance of N/H = 4.1 ± 2.1 × 10^-4 by number, i.e., 12 + log (N/H) = 8.61 ± 0.29. This value confirms the enrichment of the nebula in nitrogen and thus that the ejecta has a stellar origin. Our value computed using n_e = 350 cm^-3 is larger than the average value of Nota et al. (1996) of 12 + log (N/H) = 8.04 computed using n_e = 700 up to 2300 cm^-3 for different regions of the nebula. The Hα flux, different in both analyses, also explains the different values of the nitrogen abundance in the nebula.

Fig. 7 Temperature-density PDR diagnostic diagram. The grid of flux ratios F[O i] 63/F[O i] 146 versus ${\mathit{F}}_{\mathrm{\left[}\mathrm{O}\hspace{0.17em}{i}\mathrm{\right]}\mathrm{63}}\mathit{/}{\mathit{F}}_{\mathrm{\left[}\mathrm{C}\hspace{0.17em}{ii}\mathrm{\right]}\mathrm{158}}^{\mathrm{PDR}}$ was calculated by solving the level population equations for a range of temperatures and densities. The dashed lines corresponds to the pressure equilibrium constraint between the H ii region and the PDR, with the corresponding errors (dash-dot lines). The horizontal dashed lines correspond to the observed ratio with the errors (see Vamvatira-Nakou et al. 2013, for more details).

4.4. Ionizing flux

From the Hα flux, we calculate the rate of emission of hydrogen-ionizing photons, Q₀ using the equations provided by Vamvatira-Nakou et al. (2013). For the nebula around HD 168625, we estimate the rate of emission of hydrogen-ionizing photons to Q₀ = 1.09 ± 0.16 × 10⁴⁶ photons s^-1 with an electron temperature of T_e = 7000 K and a distance of 2.8 kpc. From the CMFGEN model, we determine a rate of emission of hydrogen-ionizing photons from the star of Q₀ = 6.9 ± 4.0 × 10⁴⁵ photons s^-1, which is in agreement, within the uncertainties, with the value obtained from the Hα flux.

The radius of the Strömgren sphere is calculated using Q₀ determined from the Hα emission. We obtain R_S = 0.13 ± 0.08 pc (~ 9.6″ ± 6.5″ at 2.8 kpc), assuming that the star is not hot enough to ionize the helium, that the ionized gas fills in the whole volume of the nebula and T_e = 7000 K. To have a radius of the Strömgren sphere equals to the radius of the ionized nebula detected from the Hα image, we must assume that the ionized gas fills in only about 20% of the entire volume of the nebula. This also leads to the conclusion that the Hα nebula is likely ionization bounded.

The mass of the ionized gas can be estimated from the Hα flux. We consider the formula given by Vamvatira-Nakou et al. (2013) to determine this mass. Since the star has T_eff = 14 000 K, and its spectrum does not show ionized helium, the mass of the ionized nebula is computed to be M_i(Hα) = 0.17 ± 0.04 M_⊙.

4.5. Photodissociation region

The [O i] 63 μm, [O i] 146 μm and [C ii] 158 μm lines probably indicate the presence of a PDR around HD 168625. This region usually surrounds the ionized region of the nebula and is composed of neutral gas. These lines can also reveal the existence of shocks, that would be the result of interactions between the fast stellar wind and the slow expanding remnant of a previous evolutionary phase. The calculated ratios between, on the one hand, the [O i] 63 μm and the [O i] 146 μm and, on the other hand, between the [O i] 63 μm and the [C ii] 158 μm are in agreement with the PDR models of Kaufman et al. (1999) and not with the shock models of Hollenbach & McKee (1989). The [O i] 63 μm/[C ii] 158 μm ratio is, furthermore, a reliable indicator of the true origin of the lines. This value is indeed found smaller than 10 in PDRs (Tielens & Hollenbach 1985) whilst it is larger for shocks. We can thus conclude that a PDR and not a shock is responsible for the presence of those lines in the spectrum of the nebula around HD 168625. We note that this conclusion is unchanged after correcting the [C ii] for the contribution of the H ii region as detailed below. Umana et al. (2010) already detected the PDR in the nebula of HD 168625 through spectral features in the [11–18] μm range indicating the presence of polycyclic aromatic hydrocarbons (PAHs). To determine the physical conditions in the PDR, we use the three infrared lines quoted above. We must however subtract any possible contribution to the observed lines of the H ii region.

Neutral oxygen can be found only in neutral regions. Therefore, no doubt exists that the [O i] lines arise exclusively from the PDR (Malhotra et al. 2001). Neutral atomic carbon has however an ionization potential (11.26 eV) lower than that of hydrogen, and can thus be found in the PDRs as well as in the H ii regions.

The C/O abundance ratio can be estimated based on the PDR line fluxes and following a method described by Vamvatira-Nakou et al. (2013) that was used to separate the contributions of the PDR and of the H ii region to the flux of [C ii] 158 μm. In the ionized gas region, the ratio of fractional ionization is given by $\frac{\mathrm{⟨}{\mathrm{C}}^{\mathrm{+}}\mathrm{⟩}}{\mathrm{⟨}{\mathrm{N}}^{\mathrm{+}}\mathrm{⟩}}\mathrm{=}\frac{{\mathit{F}}_{\mathrm{\left[}\mathrm{C}\hspace{0.17em}{ii}\mathrm{\right]}\mathrm{158}}^{\mathrm{H}\hspace{0.17em}{ii}}\mathit{/}{\mathit{ϵ}}_{\mathrm{\left[}\mathrm{C}\hspace{0.17em}{ii}\mathrm{\right]}\mathrm{158}}}{{\mathit{F}}_{\mathrm{\left[}\mathrm{N}\hspace{0.17em}{ii}\mathrm{\right]}\mathrm{122}}\mathit{/}{\mathit{ϵ}}_{\mathrm{\left[}\mathrm{N}\hspace{0.17em}{ii}\mathrm{\right]}\mathrm{122}}}$(1)where we define ${\mathit{F}}_{\mathrm{\left[}\mathrm{C}\hspace{0.17em}{ii}\mathrm{\right]}\mathrm{158}}^{\mathrm{H}\hspace{0.17em}{ii}}\mathrm{=}\mathit{\alpha }\hspace{0.17em}{\mathit{F}}_{\mathrm{\left[}\mathrm{C}\hspace{0.17em}{ii}\mathrm{\right]}\mathrm{158}}$ with F_{[C ii] 158} being the total flux of the [C ii] 158 μm line and α as a factor to be determined. Assuming that ⟨ C⁺ ⟩ / ⟨ N⁺ ⟩ = C/N and calculating the emissivities using the nebular/IRAF package with n_e = 350 cm^-3 and T_e = 7000 K, we find $\frac{{\mathit{F}}_{\mathrm{\left[}\mathrm{C}\hspace{0.17em}{ii}\mathrm{\right]}\mathrm{158}}^{\mathrm{H}\hspace{0.17em}{ii}}}{{\mathit{F}}_{\mathrm{\left[}\mathrm{N}\hspace{0.17em}{ii}\mathrm{\right]}\mathrm{122}}}\mathrm{=}\mathrm{[}\mathrm{0.37}\mathrm{\pm }\mathrm{0.24}\mathrm{]}\frac{\mathrm{C}}{\mathrm{N}}\mathrm{·}$(2)The N/O abundance ratio is not known in the literature, even though Nota et al. (1996) assumed that it is larger than 3 × (N/O)_⊙. Since the nebula characterizes the star at an earlier epoch, it seems unlikely that the N/O ratio of the nebula is larger than the surface N/O ratio measured from the star. Therefore, we define an upper limit of N/O equal to 1.4. We first assume that N/O is equal to 1. We obtain $\log \mathit{\alpha }\mathrm{=}\log \mathrm{C}\mathit{/}\mathrm{O}\mathrm{-}\mathrm{0.767}\mathit{,}$(3)using the measured F_{[C ii]}/F_{[N ii]} ratio and C/O = C/N × N/O. Assuming that there is a pressure equilibrium between the ionized gas region and the PDR, we have ${\mathit{n}}_{{\mathrm{H}}_{\mathrm{0}}}\mathit{k}{\mathit{T}}_{\mathrm{PDR}}\cong \mathrm{2}{\mathit{n}}_{\mathrm{e}}\mathit{k}{\mathit{T}}_{\mathrm{e}}\mathrm{=}\mathrm{4.9}\mathrm{\pm }\mathrm{3.7}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{\mathrm{6}}\hspace{0.17em}{\mathrm{cm}}^{-3}\hspace{0.17em}\mathrm{K}$(4)that is used to define a locus of possible values in the diagram of Fig. 7.

The structure of the PDR is described by the density of the atomic hydrogen, n_HO, and the incident FUV radiation field, G₀, that can be computed using the following equation (Tielens 2005): ${\mathit{G}}_{\mathrm{0}}\mathrm{=}\mathrm{625}\frac{{\mathit{L}}_{\mathit{\star }}\hspace{0.17em}\mathit{\chi }}{\mathrm{4}\mathit{\pi }{\mathit{R}}^{\mathrm{2}}}\mathit{,}$(5)where L_⋆ is the stellar luminosity, χ is the fraction of this luminosity above 6 eV, equal in the case of HD 168625 to 0.4 ± 0.05 (Young Owl et al. 2002) and R is the radius of the ionized gas region surrounded by the PDR (i.e., 0.22 pc). G₀ is expressed in terms of the average interstellar radiation field that corresponds to an unidirectional radiation field of 1.6 × 10^-3 erg s^-1 cm^-2. G₀ is found to be about 6.3 × 10⁴ for the PDR around HD 168625. This value coupled with the Fig. 1 of Kaufman et al. (1999) provides a temperature of the PDR of about 1600 K. This gives log n_H0 = 3.48 (Eq. (4)). From these values and Fig. 7, we derive $\log \mathrm{(}{\mathit{F}}_{\mathrm{\left[}\mathrm{O}\hspace{0.17em}{i}\mathrm{\right]}\mathrm{63}}\mathit{/}{\mathit{F}}_{\mathrm{\left[}\mathrm{C}\hspace{0.17em}{ii}\mathrm{\right]}\mathrm{158}}^{\mathrm{PDR}}\mathrm{)}\mathrm{+}\mathrm{\left[}\mathrm{C}\mathit{/}\mathrm{O}\mathrm{\right]}\mathrm{=}\mathrm{0.5}$, where by definition [C/O] = log (C/O)−log (C/O)_⊙. Using the values of the flux densities from Table 2 and (C/O)_⊙ = 0.5 (Grevesse et al. 2010), we solve the equations to obtain α = 0.09 and C/O = 0.53. By considering N/H = 4.1 × 10^-4 (Sect. 3), we have C/H = 2.2 × 10^-4. Assuming different values of N/O going from N/O = 0.4 (see Smith 2007, and the references therein) to N/O = 10 through N/O = 1.4 which is the N/O ratio determined for the central star, the resulting values of α, C/O, and C/H are listed in Table 3. We can discard N/O values that provide us with C/H ratios larger than the solar value. Moreover, since (N/H)_neb = (N/H)_star, it is reasonable to have (C/H)_neb ~ (C/H)_star. We therefore conclude that N/O is between 0.8 and 1.4 and we adopt $\mathrm{C}\mathit{/}\mathrm{H}\mathrm{=}\mathrm{1.}{\mathrm{6}}_{-0.35}^{\mathrm{+}\mathrm{1.16}}\mathrm{\times }{\mathrm{10}}^{-4}$ in the nebula.

The total mass of the hydrogen in the PDR can be estimated using the [C ii] 158 μm line flux, derived for the PDR and the equation given by Vamvatira-Nakou et al. (2013). For C/H = 1.6 × 10^-4, we estimate a mean value for the neutral hydrogen mass to M_H = 1.0 ± 0.3 M_⊙. Most of the nebular mass is thus contained in the neutral gas and not in the ionized gas.

5. Discussion

5.1. Evolutionary status

The LBV phase is still poorly understood. How and when this evolutionary phase occurs are questions that must still be debated. Currently, more and more studies aim at showing that LBVs could be binary systems. Although the term “binary” must be clarified, η Car (Damineli et al. 1997) and MWC 314 (Lobel et al. 2013) were clearly reported as binaries whilst one has detected for HR Car (Rivinius et al. 2015), and HD 168625 (Pasquali et al. 2002; Martayan et al. 2016) close objects on wide orbits that could be related to companions. If a physical orbit is confirmed for those possible companions, we can wonder whether the LBV phase could be related to binary interactions or whether the current single star evolutionary tracks can explain this evolutionary stage, as it seems to be the case for WRAY 15-751 (Vamvatira-Nakou et al. 2013) and for AG Car (Vamvatira-Nakou et al. 2015).

In the present analysis, no clear SB2 signature was found in the optical spectrum of HD 168625 but only variations of its radial velocities. We thus considered that the spectrum only comes from one star and that these variations come from pulsations (see, Sher 25, Taylor et al. 2014). We determined its log (L/L_⊙), T_eff and log g to be 5.58 ± 0.11, 14 000 ± 2000 K and 1.74 ± 0.05, respectively. From these parameters, we can estimate the location of the star on the Hertzsprung-Russell (HR) diagram and associate this location with an evolutionary track from Ekström et al. (2012). We also determine from the abundances that we derive the location where the nebula should have been ejected. We display in Fig. 8 the current location of the star and the “predicted” location of the star at the time of the ejection of the nebula. Both locations have been compared to several evolutionary tracks computed for different masses and different rotation rates. It appears that the observations agree, within the uncertainties, with a star of initial masses between 28 and 33 M_⊙ and a critical velocity v/v_crit between 0.3 and 0.4. The tracks computed without rotation or with a small rotational rate (v/v_crit < 0.3) fail at reproducing the abundances found in the nebula and in the star (no green line is observed in Fig. 8 for those tracks).

Fig. 8 Location of HD 168625 (red) in the HR diagram and location at the time of the nebula ejection estimated from the abundances of the nebula (green). The bold green lines specify the location where, in addition, the mass-loss rates is between log Ṁ = −5.0 and log Ṁ = −4.6. Evolutionary tracks are from Ekström et al. (2012).

Fig. 9 Evolution of the N abundance (in number) as a function of mass-loss rate for a 30 M⊙ star of solar metallicity and an initial rotational rate of 0.35. The evolutionary track comes from Ekström et al. (2012). The green line corresponds to the nitrogen abundance determined for the nebula from the Herschel/PACS data and the bold green line specifies where log Ṁ is between − 5.0 and − 4.6. The red line represents the location of the star according to the stellar N content and the stellar parameters from the CMFGEN model.

The N/H and the C/H abundances versus the theoretical mass-loss rate for a star with an initial mass M_init = 30 M_⊙ and a v/v_crit = 0.35 are exhibited in Fig. 9 and in Fig. 10, respectively. To create those figures, we take into account the values of the N/H and the C/H abundances that we derive from our analysis but also the stellar parameters log (L/L_⊙) and T_eff obtained from the best fit CMFGEN model. The size of the nebular shell (from the largest shell, 1.15 pc, to the ionized nebula, 0.22 pc), its total mass (M_i + M_H + M_dust ~ 0.17 M_⊙ + 0.99 M_⊙ ~ 1.17 ± 0.34 M_⊙, the mass of the dust shell being negligible) and its expansion velocity (19 km s^-1, Pasquali et al. 2002) give a mass-loss rate between log Ṁ = −5.0 and log Ṁ = −4.6. This allows us to refine the epoch of the ejection of the nebula (bold green line in Figs. 8−10). This epoch corresponds to the end of the main-sequence phase, and, from the current location of the central star, it also appears that the star has not yet reached the red supergiant (RSG) phase. We note however that possible rapid excursions toward higher mass-loss rates that are not expected by current evolutionary tracks could also explain such values.

Fig. 10 Same as for Fig. 9 but for carbon abundance.

5.2. Evolutionary scenarios

5.2.1. Single star scenario

Single star evolutionary tracks can explain the evolution of the N content between the nebula and the central star but this still does not explain the shape and the different shells observed of the nebula around HD 168625. A few objects exhibit an analogue bipolar structure to SN 1987A (Plait et al. 1995; Mattila et al. 2010). Similar polar empty rings have been observed for Sher 25 (Hendry et al. 2008), SWB1 (Smith et al. 2013), or MN18 (Gvaramadze et al. 2015) but the nature of their formation still remains unknown.

It is however worth noticing that none of these objects has reached the RSG phase. Smartt et al. (2002) indeed showed from the determination of the surface abundances of Sher 25 that this star has not gone through a RSG phase. Smith et al. (2013) estimated the N content of the ring of SWB1 to the solar value and came to the conclusion that this star has also not gone through a RSG phase. Gvaramadze et al. (2015), from a CMFGEN analysis, determined that MN18 is a redward evolving star that became a BSG only recently. Therefore, the scenario proposed by Chita et al. (2008) aiming at explaining the creation of the bipolar structure from interactions between a BSG wind and a RSG thin shell seems thus compromised. Moreover, this theory fails at representing the empty polar rings.

From our analysis, nothing prevents us to think that the larger polar rings (detected on the 8 μm image) have the same age as the equatorial ring/torus. In this case, as already mentioned by Smith (2007), they could not have been ejected in a RSG phase because they would be far too fast. Therefore, since the creation of these two structures (the equatorial ring and the polar caps) suggests a common ejection occuring during the BSG phase following the main sequence. According to Smith et al. (2013), the polar rings would be formed by a dusty ionized flow of gas photoevaporated from the equatorial dense ring/torus. Its pressure prevents the stellar wind from moving towards the equatorial shell so that it would be deviated toward the poles, creating the polar empty rings.

Rotation has also been claimed as a possible cause to explain the formation of the bipolar structure (see e.g., Gvaramadze et al. 2015). This scenario assumes that a fast stellar wind (v_∞ ~ 350 km s^-1 for HD 168625) concentrated in the equatorial plane because of the rotation collides with a slower and denser wind that moves away from the star at smaller velocities (i.e., 19 km s^-1 for HD 168625). The outflow is then collimated in the polar direction by the presence of this dense equatorial material. The disk is however expected by assuming a large rotation velocity, which is not the case for HD 168625 (vsini = 60 km s^-1) (see Bjorkman & Cassinelli 1993, for further details). Furthermore, hydrodynamical simulations were performed by Smith et al. (2013, and the references therein) and failed at representing empty polar rings, but rather complete polar caps.

Finally, it is also possible that the illumination of the different parts of the nebula is anisotropic, indicating that some parts of the nebula could block some lines-of-sight, and would expect to result in a bipolar appearance.

5.2.2. Binary system scenario

Although no clear evidence has been found that HD 168625 is a binary system, the detection of the presumably wide-orbit companion by Pasquali et al. (2002) and Martayan et al. (2016) (even though no clear motion of this “companion” is detected between the two epochs) and the possible role of a binary/multiple system in the evolution of this object remain questionable. Moreover, HD 168625 is quoted as a progenitor of SN 1987A that has been thought to be a binary merger by Morris & Podsiadlowski (2009). Therefore, we must envisage the binarity scenario to explain the ejection of this bipolar structure. The binary evolutionary tracks being complex to calculate, we were not able to compare the CNO abundances that we have determined for the central star and for the nebula to those tracks to see whether a match could be found.

As the abundances determined for the central star are in agreement with those estimated from a single star evolutionary track, it does not allow us to consider a mass transfer by Roche lobe overflow. Even though such a transfer happened, we expect that the rotational velocity of the mass gainer increases whilst that of the mass donor decreases. Given the projected rotational velocity, it seems unlikely that HD 168265 is the gainer unless it loses a significant quantity of angular momentum on a very short timescale. We indeed see that, for HD 168625, the projected rotational velocity that we have determined from the optical spectrum is not large enough to classify it as rapid rotator. If we assume an inclination of about 40° or 60° (both values estimated from the imaging analysis), the rotational velocity of HD 168625 is about 70−95 km s^-1. We also see that, for the other possible SN progenitors that show similar nebulae in the polar direction, their projected rotational velocity is also quite small. As outlined by Gvaramadze et al. (2015), the most efficient mechanisms to spin down a massive star is the magnetic braking via a magnetized wind but that mechanism requires magnetic fields of about 10³−10⁴ G and a high mass-loss rate. HD 168625 has been probed by Fossati et al. (2015) to look for magnetic field but these authors did not detect any significant magnetic field. If a transfer of material by Roche lobe overflow occurs in HD 168625, we should also see emissions in the X-ray domain. However, Nazé et al. (2012) did not detect any significant X-ray emission for this star. Therefore, we must consider another form of exchange: the binary merger. In this case, this scenario could explain the absence of X-ray but would imply rapid rotation, which is not found in HD 168625.

6. Conclusion

We have presented the analysis of the Herschel/PACS imaging and spectroscopic data of the nebula surrounding the candidate LBV HD 168625 together with Hα and mid-infrared images as well as high-resolution optical spectra. The images, covering the Hα–160 μm domain, confirm the complex morphology of the nebula, composed of torus-like shells visible at different wavelengths and of a perpendicular bipolar structure that is observed in the Hα, 8 μm and 13 μm images.

We have also seen that the ionized nebula is smaller than the entire nebular shell. We determined in the present paper that about 85% of the total mass of the nebula comes from the neutral gas (M_H ~ 1.0 M_⊙ out of 1.17 M_⊙ for the total mass) against 14% from the ionized shell (M_i ~ 0.17 M_⊙), the remaining percentage being attributed to the dust.

The far-infrared spectrum of the nebula provided us with the nitrogen and carbon abundances. We have determined an overabundance in nitrogen (N/H ≅ 4.1 × 10^-4) and a depletion in carbon (C/H ≅ 1.6 × 10^-4) in the nebula around HD 168625. We compared these values to the CNO abundances determined with the CMFGEN atmosphere models computed from optical spectra of the central star. We were able to constrain the evolutionary path of the star and to determine the moment of the nebula ejection. This star seems to have an initial mass between 28 and 33 M_⊙ and should have ejected its material during or just after the Blue Supergiant phase and before the red supergiant stage. The evolution of HD 168625 seems thus related to the evolution of a single star even though the complex structure of its bipolar/biconal rings remains unexplained under this assumption and could therefore favor the binary scenario.

Acknowledgments

We thank the anonymous referee for his/her comments and remarks. We also thank Prof. Gregor Rauw for their valuable comments on the manuscript, Pr. D. J. Hillier for making CMFGEN available and Dr. Nick Cox for providing the scanamorphos images. L.M., D.H., P.R. acknowledge support from the Belgian Federal Science Policy Office via the PRODEX Program of ESA. The Liège team also acknowledges support from the Fonds National de la Recherche Scientifique (F.R.S.-F.N.R.S.). This research was also funded through the ARC grant for Concerted Research Actions, financed by the French Community of Belgium (Wallonia-Brussels Federation). PACS has been developed by a consortium of institutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF-IFSI/OAA/OAP/OAT, LENS, SISSA(Italy); IAC (Spain). This development has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX (Belgium), CEA/CNES (France), DLR(Germany), ASI/INAF (Italy), and CICYT/MCYT (Spain). Data presented in this paper were analyzed using “HIPE”, a joint development by the Herschel Science Ground Segment Consortium, consisting of ESA, the NASA Herschel Science Center, and the HIFI, PACS and SPIRE consortia. This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, as well as NASA/ADS and SIMBAD (CDS/Strasbourg) databases.

References

Appendix A: Journal of observation of the optical spectra

Appendix B: Emission line fluxes for each spaxel

All Tables

All Figures

Fig. 1 Comparison between the FEROS spectrum taken on June 29, 2002 (black line), the UVES spectrum taken on July 10, 2001 (red line) and the X-shooter spectrum taken on June 2, 2014 (blue line). The He i 4471 and the Mg ii 4481 lines are shown as effective temperature diagnostic, the N ii lines between 4600 and 4625 Å as nitrogen content indicator, and the Hα line for the stellar wind variations.

In the text

	Fig. 2 Best-fit CMFGEN model for HD 168625 (red line) compared to FEROS spectrum taken on June 29, 2002 (black line). The lines that are not modeled are diffuse interstellar bands (DIBs). The Hα line has been rescaled for clarity.
In the text

Fig. 3 Top: from left to right: Hα image taken with SUSI on NTT (left), the 8 μm Spitzer image, and the 12 μm VLT/VISIR image. Middle: the MIPSGAL 24 μm image, the PACS 70 μm image and the PACS 100 μm image from left to right. Bottom: the PACS 160 μm images. For all the images, north is up and east is left. The white arrow indicates the location of the LBV HD 168607. In the Hα image of Nota et al. (1996), the crosses indicate the locations of the maxima of emission detected in the radio image of Leitherer et al. (1995). We emphasize that the 12 μm image is from Umana et al. (2010) and that the inner part of the 24 μm image is saturated.

In the text

	Fig. 4 Three-color image (70 μm in blue, 100 μm in green and 160 μm in red) of the HD 168625 nebula. North is up and east is to the left.
In the text

	Fig. 5 Footprint of the PACS spectral field-of-view on the image at 70 μm of the nebula. The spaxels are represented in gray. North is up and east is left.
In the text

	Fig. 6 PACS spectrum of the HD 168625 nebula, integrated over the 25 spaxels. The different bands are indicated with different colors. The continuum and the line at 205 μm are not calibrated, the inset is expressed in arbitrary units. The spectrum between 95 and 102 μmis affected by an instrumental artefact that makes the spectrum scientifically unusable in this wavelength domain.
In the text

Fig. 7 Temperature-density PDR diagnostic diagram. The grid of flux ratios F[O i] 63/F[O i] 146 versus ${\mathit{F}}_{\mathrm{\left[}\mathrm{O}\hspace{0.17em}{i}\mathrm{\right]}\mathrm{63}}\mathit{/}{\mathit{F}}_{\mathrm{\left[}\mathrm{C}\hspace{0.17em}{ii}\mathrm{\right]}\mathrm{158}}^{\mathrm{PDR}}$ was calculated by solving the level population equations for a range of temperatures and densities. The dashed lines corresponds to the pressure equilibrium constraint between the H ii region and the PDR, with the corresponding errors (dash-dot lines). The horizontal dashed lines correspond to the observed ratio with the errors (see Vamvatira-Nakou et al. 2013, for more details).

In the text

	Fig. 8 Location of HD 168625 (red) in the HR diagram and location at the time of the nebula ejection estimated from the abundances of the nebula (green). The bold green lines specify the location where, in addition, the mass-loss rates is between log Ṁ = −5.0 and log Ṁ = −4.6. Evolutionary tracks are from Ekström et al. (2012).
In the text

Fig. 9 Evolution of the N abundance (in number) as a function of mass-loss rate for a 30 M⊙ star of solar metallicity and an initial rotational rate of 0.35. The evolutionary track comes from Ekström et al. (2012). The green line corresponds to the nitrogen abundance determined for the nebula from the Herschel/PACS data and the bold green line specifies where log Ṁ is between − 5.0 and − 4.6. The red line represents the location of the star according to the stellar N content and the stellar parameters from the CMFGEN model.

In the text

	Fig. 10 Same as for Fig. 9 but for carbon abundance.
In the text