Free Access
Issue
A&A
Volume 599, March 2017
Article Number A61
Number of page(s) 17
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201629585
Published online 01 March 2017

© ESO, 2017

1. Introduction

With its Of?p spectral type, HD 148937 belongs to a small class of peculiar stars that has triggered quite some interest in recent years. This category of stars was defined by Walborn (1972) to designate stars displaying some similarities to Of supergiants but also clearly different characteristics – in particular, the presence of strong C iii 4650Å emission lines. HD 148937 was amongst the first objects classified as Of?p, but was only recently the target of in-depth studies.

HD 148937 is the brightest and hottest massive star to possess a strong magnetic field (Martins et al. 2012). The effective temperature of this object reaches ~40 000 K, its gravity log g ~ 4.0, and its mass approximately 50–60 M (Nazé et al. 2008a; Martins et al. 2015). It displays chemical enrichment, but with a level comparable to O supergiants (Martins et al. 2015). It is a slow rotator, as evidenced by its narrow spectral lines (vsini ≤ 45 km s-1, Nazé et al. 2010; Wade et al. 2012; Martins et al. 2015). Its magnetic field was discovered by Hubrig et al. (2008) and was subsequently monitored by Wade et al. (2012) who derived a bipolar field strength of 1kG, confirming the pole-on geometry derived by Nazé et al. (2010).

The strong magnetic field is able to channel the stellar wind of HD 148937 towards its (magnetic) equator. As the magnetic axis is tilted with respect to the rotation axis, the dense equatorial regions are seen under different angles during the stellar rotational cycle. This explains why the star presents line profile variations of its Balmer and He i lines, which are partially formed in this confined wind region. The spectral changes detected for HD 148937 are, however, of a smaller amplitude than for other Of?p stars (Nazé et al. 2008b). This is explained by a nearly constant pole-on configuration (Nazé et al. 2010) that was confirmed later with spectropolarimetric observations (Wade et al. 2012). In view of their interpretation, these changes should be periodic and they indeed are, with a period of only 7.03 days (Nazé et al. 2008b), whilst the other Galactic Of?p stars display much longer periods (73 days for CPD28° 2561, 158 days for NGC 16242, 538 days for HD 191612, and 55 yr for HD108). Finally, HD 148937 emits an intense X-ray emission which is one order of magnitude brighter than “normal” O-stars but similar to other magnetic O-stars (Nazé et al. 2010, 2012, 2014). The derived location of the emission region (~1 R above the photosphere) as well as the X-ray strength are compatible with hot plasma arising from the head-on collision between channeled wind flows from the two hemispheres, even though the softness of the emission, compared to e.g. θ1 Ori C, remains unexplained (Nazé et al. 2012, 2014).

HD 148937 is located at the center of NGC 6164/5, a bipolar nebula that was studied by Leitherer & Chavarria (1987) and Dufour et al. (1988). The latter showed that the composition of the bipolar structure NGC 6164/5 is chemically peculiar with an overabundance in nitrogen and in helium and a depletion in oxygen and in neon.

The kinematics of the nebula was analyzed in detail in several studies and seems to be also peculiar. Catchpole & Feast (1970) found that the brightest lobes (constituting the brightest parts of NGC 6164 and NGC 6165) of the nebula were moving nonisotropically with radial velocities (RVs) of about + 21 and −43 for the NW and the SE lobes, respectively. These results were refined and confirmed by Pismis (1974) and later by Carranza & Agüero (1986) and Leitherer & Chavarria (1987). However, the kinematics of the nebula as a whole is not as simple as this would suggest because the inner parts of the nebula present an independent kinematic structure with different parts approaching and/or receding.

When one looks carefully at the interstellar environment that the star is embedded in, the star and the nebula appear further surrounded by an ellipsoidal stellar-wind-blown shell with a radius of about 12′ and, further away, by a Strömgren sphere with a radius of about 1.1° (Leitherer & Chavarria 1987). The origin of such complex motions, as well as the nature of the global morphology, are still open questions.

In the present paper, we investigate the Spitzer IRS and Herschel/PACS data of the nebula to further constrain the properties of NGC 6164/5. We also put these results in perspective through a new modeling of the central star with the CMFGEN atmosphere code (Hillier & Miller 1998). The abundances determined in both environments are compared to trace back the evolution of HD 148937 and to understand how such a nebula was created. The paper is organized as follows. In Sect. 2, the observations and data reduction procedures are presented. In Sect. 3, we describe the global morphology of the nebula and a general summary of its kinematics is provided in Sect. 4. We analyze the Spitzer IRS and the Herschel/PACS spectra obtained at nine different positions and in two different regions across the nebula, respectively, in Sect. 5. We then model the central star with an atmosphere code in Sect. 6. We discuss the results and constrain the evolution as well as the origin of the nebula in Sect. 7. Finally, we provide the conclusions in Sect. 8.

2. Observations and data reduction

2.1. Nebula observations

2.1.1. Infrared data

The far-infrared observations were obtained, with the PACS photometer (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 PACS imaging observations were carried out on September 2nd, 2011, which corresponds to Herschel’s observational day (OD) 842. The observing mode was the scan map in which the telescope slews at constant speed (in our case the “medium” speed of 20″/s) along parallel lines in order to cover the required area of the sky. Two orthogonal scan maps were obtained for each filter. The observation identification numbers (obsID) of the four scans are 1342227809, 1342227810, 1342227811, and 1342227812, each having a duration of 889 s.

We use the Herschel interactive processing environment (HIPE, Ott 2010) to perform the data reduction up to level 1. Subsequently, we apply the Scanamorphos software (Roussel 2013) to reduce and combine the data with the same wavelength. The pixel size in the final maps is 2″ in the blue channel (70 and 100  μm) and 3″ in the red channel (160  μm). The 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.

To complete this imaging data set, we also retrieved the WISE infrared images (Wright et al. 2010) at 12  μm and at 22  μm to better constrain the morphology and the properties of the nebula at shorter wavelengths.

The Herschel/PACS spectroscopic data of the nebula have been obtained with the PACS integral-field spectrometer that covers the [52–220]  μm wavelength region. This instrument operates simultaneously in two channels: the blue one ranging from 52 to 98  μm and the red one going from 102 to 220  μm. The nebula is, however, too large with respect to the spectrometer field of view. Therefore, two on-source spectra have been taken, one at the coordinates RA = 16h 33m 44.5s and δ = −48° 06′ 27.4″ and one at coordinates RA = 16h 33m 42.7s and δ = −48° 04′ 50.0″ (see the upper panel of Fig. 1). The former region will be referred to in the present paper as the Herschel-1 region (H1) whilst the latter is reported as the Herschel-2 region (H2). In addition, an off-source spectrum has been taken at the coordinate RA = 16h 30m 33.0s and δ = −48° 07′ 01.1″. These spectra allow us to derive the properties of the nebula, in particular the abundances. For each spectrum, simultaneous imaging of a 47″ × 47″ field of view is obtained, resolved in 5 × 5 square spatial pixels (i.e., spaxels). The two-dimensional field of view is then rearranged along a 1 × 25 pixel entrance slit for the grating via an image slicer employing reflective optics. The resolving power is between 940 and 5500 depending on the wavelength. The two obsIDs related to the spectrum taken in the H2 region are 1342252090, and 1342252091 whilst the two other obsIDs related to the data taken in the H1 region are 1342252092, and 1342252093. The data were reduced with HIPE version 14, using the interactive pipeline script for unchopped line spectroscopy with calibration set 77. In unchopped mode, the level of the continuum is not guaranteed and has no physical meaning. Therefore, we subtract the continuum to only represent the emission lines.

thumbnail Fig. 1

Top: general picture of the nebula NGC 6164/5 and definition of the nomenclature used in the present paper. North is up and east is left. The small white circles represent the positions of the Spitzer IRS spectra whilst the green squares give the footprints of the Herschel/PACS field of view. The NW and SE lobes are indicated by gray arrows and the inner zones of ejecta are represented in red. Bottom: Hα image of NGC 6164/5 taken from the SuperCosmos Hα Survey (Parker et al. 2005). The white boxes represent the areas of the nebula analyzed by Leitherer & Chavarria (1987) whilst the red arrows represent the positions analyzed by Dufour et al. (1988).

Nine infrared spectra of HD 148937 (AORkey 18380032, PI: Morris) and of its nebula NGC 6164/5 (AORkey 18380800, PI: Morris) were obtained with the Infrared Spectrograph (IRS; Houck et al. 2005) onboard the Spitzer Space Telescope (Werner 2005). The observation of the nebula (AORkey 18380800, #1 to #8 in the upper panel of Fig. 1) was performed in “cluster” mode in which spectra are taken at eight pointed positions whilst the observation of HD 148937 (AORkey 18380032, #9 in the upper panel of Fig. 1) was performed in “single pointing” mode. The reduced spectra were downloaded from the Cornell Atlas of Spitzer/IRS Sources (CASSIS; Lebouteiller et al. 2011, 2015). We use the high-resolution (R = λ/δλ = 600) spectra observed with the Short-High (SH) and Long-High (LH) modules, covering spectral ranges from 9.9 to 19.6  μm (SH) and from 18.7 to 37.2  μm (LH), respectively. The aperture sizes are 4.7″ × 11.3″ and 11.1″ × 22.3″ for the SH and LH spectra, respectively. The spectral extraction uses either a full aperture extraction for which the flux is simply integrated within the aperture or an optimal extraction for which the flux is found by scaling the point-spread function to the source’s spatial profile. While the former method is well suited for extended sources, the latter method is only reliable for point-like sources. For this reason, for the eight spectra of the nebula, we used the full aperture extraction, as the emission is always extended in our observations and the optimal extraction for the spectrum taken at the position of HD 148937 (# 9 in the upper panel of Fig. 1). Unfortunately, we were not able to analyze the LH spectrum in the optimal extraction mode because the background emission was too strong to allow its extraction. We therefore only use the SH spectrum at that position. For the spectra across the nebula, a scaling factor has been applied to the SH spectra in order to align it with the continuum of the LH spectrum. This scaling factor is requested to quantify the correction for aperture losses and is reported for each spectrum in Table 1.

Table 1

Line fluxes measured on the Spitzer IRS spectra.

thumbnail Fig. 2

Images of NGC 6164/5 around HD 148937. Top: the Gemini Hα + [N ii] image (left), the WISE 12  μm and 22  μm images (middle and right, respectively). Bottom: PACS images of the nebula at 70  μm, 100  μm and 160  μm from left to right, respectively. The white asterisk indicates the position of HD 148937. The dash-lined box provides the size of the Hα image because the sizes of these images are different. North is up and east is left.

The locations at which the different infrared spectra have been taken are displayed in the upper panel of Fig. 1. As a comparison, and because we use it later in the present paper, we also show the areas of the nebula analyzed by Leitherer & Chavarria (1987) and Dufour et al. (1988) in the lower panel of Fig. 1.

2.1.2. Optical data

Optical images of NGC 6164/5 have been retrieved from the Gemini Science Archives. A series of 120 s exposures were obtained with GMOS-S imager in three different filters: [S ii] (λc = 6720 Å, FWHM = 44.4 Å), [O iii] (λc = 4990 Å, FWHM = 44.3 Å) and Hα + [N ii] (λc = 6560 Å, FWHM = 71.6 Å). Unfortunately, no continuum image was taken along with those nebular images. The standard reduction has been performed.

To calibrate the flux of the nebula on the Hα image, we used the HST image in the Hα filter obtained in 1991 with the WFPC equipped with the F656N filter. This filter is centered on Hα but the neighbouring [NII] lines also contribute to the recorded emission (to the level of 30%). Reduced HST observations were downloaded from the archives. They focused only on the southeast lobe of the nebula (NGC 6165) and were reported by Scowen et al. (1995). These data were combined to make a mosaic using the task wmosaic under IRAF and the flux calibration is obtained by multiplying the recorded ADUs by the PHOTFLAM keyword and the filter width (19 Å), by dividing by the exposure time (2 × 800 s) and by correcting for 30% of the [N ii] contribution.

2.2. Stellar observations

To determine the stellar parameters and the abundances of HD 148937, we also retrieved the spectra of HD 148937 taken with FEROS and ESPaDOnS presented by Nazé et al. (2008b) and Wade et al. (2012), respectively. The analysis of the variations observed in the spectra has been done in both papers and will not be repeated in the present study. The reduction processes of both sets of data were the same as explained in the two papers.

3. Morphology of the nebula

The Hα and the infrared images of HD 148937 and of its surrounding nebula NGC 6164/5 at 12  μm, 22  μm, 70  μm, 100  μm, and 160  μm are shown in Fig. 2.

The Hα image reveals an ionized nebula, oriented to the NW-SE axis, with a bipolar or “8” shape centered on the Of?p HD 148937. This ionized nebula is composed of different visible structures: two brighter lobes located at SE and NW edges of the nebula (in blue and indicated by gray arrows in the upper panel of Fig. 1) and three fainter zones of ejecta located closer to the central star (in red in the upper panel of Fig. 1). The whole nebula is formed by NGC 6164 at NW and NGC 6165 at SE (see the large white circles in Fig. 1). In the ionized nebula, the maxima of emission in these two regions are located at about 150″ and 185″ from HD 148937 in the NW and the SE directions, respectively. Assuming that the distance between the star and Earth is about 1.3 kpc (the distance of Ara OB1, Herbst & Havlen 1977, which HD 148937 is assumed to be embedded in), the projected separation would then be equal to 0.95 pc and 1.16 pc for NGC 6164 and NGC 6165, respectively. Leitherer & Chavarria (1987) estimated the inclination between the plane in which the nebula has been ejected and the line of sight to 85° even though, according to these authors, the inclination can be assumed to be between 60° and 90°.

thumbnail Fig. 3

Three-color image of PACS photometer field of view with 70  μm in blue, 100  μm in green and 160  μm in red. North is up and east is left. The box indicates the location of the Hα image (see Fig. 2). The dashed line represents the location of the wind-blown structure at 12′ from HD 148937 (seen on the Hα image).

In the infrared images, the lobes become fainter whilst the three inner zones of ejecta become brighter until the nebular emission is not distinguishable from the external environment at wavelengths longer than 70  μm. This observation suggests that the dust emission is weaker in the lobes although the presence of cooler dust in the lobe cannot be excluded. Indeed, at the inner edges of the lobes, one might detect, from the Hα image, a zone of collision. According to van Marle et al. (2011), only large grains can pass through this zone. By combining a larger distance which the lobes are located at and a possible larger size of dust grains in the lobes, we expect lower dust emission in the lobes at wavelengths shorter than 70  μm, as observed. In the 22  μm image, we see that the zone of ejecta (in red in Fig. 1) located at west is particularly bright in comparison to the others. This zone is in fact contaminated by the presence of a dark cloud seen in projection (Peretto & Fuller 2009), that makes the analysis of the NGC 6164/5 more complex. Because of this dark cloud and of the high background emission, we cannot measure accurate flux densities from NGC 6164/5 only and thus determine an accurate value of the mass of the dust in the nebula NGC 6164/5. Finally, it is worth noting that the infrared images also show a cavity centered on HD 148937. This lack of emitting material has probably been formed through the action of the stellar wind. The shape and the size of this cavity are different from one infrared image to another: it has a horseshoe shape in the 12  μm, 22  μm, and 70  μm images but appears larger and with an unconstrained shape in the images at longer wavelengths. We can roughly estimate lower limits on the angular distances of 35″, 51″, and 64″, giving distances of 0.22, 0.32 and 0.40 pc between the star and the eastern, western and northern inner zone of ejecta, respectively, assuming a distance between the star and the Earth of 1.3 kpc.

In Fig. 3, we show a larger-scale view using a three-color image based on the PACS data (70  μm in blue, 100  μm in green and 160  μm in red). The nebula around HD 148937 is barely detectable at the bottom of this image (white box in Fig. 3) and does not constitute the brightest part of the image. Given the position of the nebula and the size of the observed field, it is not possible from this image to detect the Stromgren sphere at 1.1° observed by Leitherer & Chavarria (1987) whilst the wind-blown structure at 12′ from HD 148937 is partially in the field of view of PACS (dashed line in Fig. 3), but cannot be distinctly observed because of the background emission.

thumbnail Fig. 4

Spitzer IRS spectrum taken at position #2 (see Fig. 1). The step at 19  μm is due to the different apertures used for the SH and LH spectra.

4. Kinematics of the nebula

The kinematics of the whole ionized Hα nebula was analyzed by Catchpole & Feast (1970), Pismis (1974), Carranza & Agüero (1986) and Leitherer & Chavarria (1987). These different authors measured the RVs of the ejecta in four different regions: the NW and SE lobes as well as in the western and the eastern bright parts. Their measurements were, in general, in good agreement. Before comparing the RVs of these different parts, it is worth noting that the velocity of the center-of-mass of HD 148937 was measured to be −26 km s-1 (Nazé et al. 2008a). A part of the NW lobe (constituting NGC 6164) is moving away from the observer with RVs between + 21 and + 55 km s-1 but another part of this region is also approaching the observer with a RV of −3 km s-1. The SE lobe (forming NGC 6165) is moving toward the observer with RVs between −30 and −145 km s-1. The kinematics of the ejecta located between both lobes is more complex. Leitherer & Chavarria (1987) reported Hα profiles reproduced with four Gaussian profiles, giving four different RVs for the ejecta located at the western part of the nebula (RVs = −70, −19, + 33, and + 63 km s-1 were reported) whilst they used only two Gaussian profiles to fit the Hα profiles observed in the eastern part of the nebula, giving RVs of −79 and −13 km s-1. This complex kinematics suggests that the bipolar or “8”-shaped nebula displayed in the Hα image could thus be the result of the projection of a helical structure centered on HD 148937, as explained in further detail by Carranza & Agüero (1986).

5. Emission line spectrum

5.1. Spitzer IRS spectra

Eight mid-infrared spectra have been taken at different positions across the nebula and one on HD 148937. As previously mentioned, the exact locations of these nine spectra are displayed in Fig. 1. These data show spectral lines that can provide important constraints on the physical conditions in the nebula. As an example, we show in Fig. 4 the merged Spitzer IRS spectrum (SH+LH) taken at position #2 (see the upper panel of Fig. 1) and where the SH spectrum is not corrected for the scaling factor. In addition, because HD 148937 is a point-like source, its LH spectrum was extracted under the optimal mode. The latter is, however, not usable because the source is highly contaminated by the background level. Therefore, the measurements for HD 148937 are taken only on the SH spectrum.

We measure the emission line flux densities by fitting a Gaussian profile to the main line profiles. The following emission lines are detected in the different spectra: [S iv] 10.5  μm, Hu α 12.4  μm, [Ne ii] 12.8  μm, [Ne iii] 15.5  μm, [S iii] 18.7  μm, [Ar iii] 21.8  μm, [Fe iii] 22.9  μm, H2 28.2  μm, [S iii] 33.5  μm and [Si ii] 34.8  μm as well as PAH features at 11.3  μm. The correction factors as well as the corrected fluxes measured for the different lines are listed in Table 1. We also emphasize that the Spitzer position #8 is located outside the nebula, on an infrared source. The measured fluxes are thus contaminated by this external source and are thus not representative of the fluxes measured across the nebula.

5.1.1. Electron density

The ratio between the [S iii] 33.5  μm and [S iii] 18.7  μm lines can be used to derive the electron density in the ionized nebula surrounding HD148937.

Table 2

Electron temperatures, electron densities, and abundances across the nebula.

To determine the emission line intensities of the [S iii] 18.7  μm line, we must apply a correction factor to take the aperture losses into account. This correction factor is estimated by aligning the continuum flux of the SH spectra to that of the LH spectra in the small overlap between the two spectra, except for the region near the star (#9) where this correction is impossible (see Sect. 2.1.1). In addition to the error on the ratio between the intensities of the [S iii] lines (~5%), the global error on the determination of the electron density is dominated by the errors on the correction factor and on the electron temperature.

To measure the electron densities in the different regions, we use the nebular/IRAF package. To this purpose, we need an estimate of the electron temperature across the nebula. Leitherer & Chavarria (1987) as well as Dufour et al. (1988) investigated the electron temperature at different locations (not necessarily the same positions as the Spitzer IRS spectra, see Fig. 1) in the nebula. Even though the electron temperatures close to the lobes are in agreement between these two studies, the other values obtained in the other parts of the nebula differ from one analysis to another. Furthermore, Leitherer & Chavarria (1987) provided two upper limits for the electron temperature. We try to use these upper limits to determine the electron densities in our spectra but in some cases, these upper limits are too high to derive a physical value for the electron density from the observed line ratios. Therefore, we consider in this section as well as for our entire analysis a global electron temperature of 7000 K (compatible with values reported by these authors). We also assume a conservative uncertainty of 2500 K on the electron temperature.

From this value, we obtain very small electron densities across the nebula except for the lobes where we find ne = 300 ± 70 cm-3 and ne = 930 ± 70 cm-3 in the central parts of the NW and the SE lobes, respectively. Elsewhere, at the different locations of the Spitzer IRS spectra, the electron densities are closer to 100 cm-3. At position #4, the electron density is relatively low, at the limit of measurability. We note that, if we decrease Te to 5000 K, the electron density in this area becomes 70 ± 20 cm-3. We emphasize that we have determined the rms densities on the Hα image at the same positions where the Spitzer data were taken to test our estimates on the electron density. These rms densities can be taken as lower limits on the electron density but are dependent on the volume that we consider to measure them. As first approximation, we have assumed spherical and cylindrical volumes with radii equal to the aperture sizes of Spitzer. We note a global agreement between the rms densities and the electron densities, thus validating our estimations. The derived values of the electron densities are given in Table 2.

5.1.2. Abundances

The Spitzer IRS spectra do not provide information on CNO abundances across the nebula. They do, however, provide abundances for elements such as sulfur, argon, or neon. To take homogeneous measurements of the hydrogen fluxes, we focus on the flux estimated from the Hu α hydrogen line at 12.37  μm. We then assume a case B recombination with Te = 7000 K. Table 2 contains the abundances that we derive for the nine Spitzer IRS spectra.

To derive the abundances, we use the nebular/IRAF package. The total elemental abundances are difficult to derive because certain elements do not have measurable lines in every ionization state in the Spitzer wavelength domain. Therefore, we consider that the ionic abundances such as Ar++/H+, or S++/H+ are lower limits to the real abundance values. We nevertheless obtain abundances for argon, neon and sulfur compatible with those provided by Dufour et al. (1988). Furthermore, the abundances of these three elements are relatively constant across the nebula, with Ar/H compatible with the solar abundance whilst Ne/H and S/H could be considered as depleted in both the inner and outer regions of the nebula.

5.2. PACS spectra

Given the size of the nebula, only two areas of the nebula (reported as H1 and H2 in Fig. 1) have been observed by PACS. The footprints of the PACS spectral field of view on the nebula as well as the off-source position are shown in Fig. 1. Each PACS spectral field of view is composed of 25 (5 × 5) spaxels, each corresponding to a different part of the nebula. We measure the emission line intensities, in each one of the 25 spaxels present in the two on-source, as well as in the off-source spectra (see the tables in Appendices A, B and C), by fitting a Gaussian profile to every line profile. We also measure the intensities of these lines on the integrated spectra over the 25 spaxels. The final values of the line intensities, measured on these integrated spectra correspond to the subtraction between the on-source and the off-source emission line fluxes. They are given in Table 3. The resulting spectra are displayed in Figs. 5 and 7 for the H1 and for the H2 regions, respectively, whilst the comparisons between the on-source and the off-source spectra are displayed in Figs. 6 and 8 for the H1 and the H2 regions, respectively.

The following forbidden emission spectral lines are detected in both spectra: [N iii] 57  μm, [O i] 63  μm, [O iii] 88  μm, [N ii] 122  μm, [O i] 146  μm, [C ii] 158  μm and [N ii] 205  μm. The presence of high ionization lines indicates that the nebula around HD 148937 is highly ionized whilst the lowest ionization lines seem to reveal the existence of shocks or of a photodissociation region (PDR) around the central star.

thumbnail Fig. 5

PACS spectral lines of the H1 region in the nebula around HD 148937, integrated over the 25 spaxels and background corrected, using the off-source region. The continuum has been removed for clarity.

thumbnail Fig. 6

PACS spectral lines of the H1 region (in solid line) and of the off-source region (in dashed line) integrated over the 25 spaxels. The continuum has been removed for clarity.

thumbnail Fig. 7

As for Fig. 5 but for the H2 region.

thumbnail Fig. 8

As for Fig. 6 but for the H2 region.

5.2.1. The H1 region

Electron density.

As we already mentioned in Sect. 5.1.1, the electron temperature has been assumed to be 7000 K through the entire nebula. To determine the electron density in the H1 region, we use the nebular/IRAF package. We compute a ratio of 2.33 ± 0.53 between the [N ii] 122  μm and the [N ii] 205  μm lines on the integrated spectrum (Table 3). In our spectral range, these lines are good indicators of the electron density. Assuming an electron temperature of 7000 K, we determine an electron density of ne = 65 ± 25 cm-3. We note that ne increases to 70 cm-3 if we consider the upper limit on the electron temperature (16 600 K) reported by Leitherer & Chavarria (1987) for the H1 region. This remains, however, within the error bars. The electron density found in the H1 region is in agreement with the value obtained from the Spitzer data at position #6.

Table 3

Line fluxes measured on the integrated spectra of the H1 and H2 regions of NGC 6164/5, after subtracting the off-source values.

Hα flux.

To determine the Hα fluxes across the nebula, we use the Hα image taken with Gemini, recalibrated using the HST data (see Sect. 2). We use the photometric aperture of the same size as the PACS spectrometer field of view. The reddened flux measured in the H1 region is estimated to F0(Hα) = 1.06 × 10-11 erg cm-2 s-1. By assuming (BV) = 0.343 (Maíz-Apellániz et al. 2004) and (BV)0 = −0.27 (given the spectral type O 5.5f?p of HD 148937, Martins & Plez 2006), we compute E(BV) = 0.61. From this value, we determine de-reddened fluxes of F0(Hα) = 4.38 ± 0.5 × 10-11 erg cm-2 s-1 using the extinction law of Cardelli et al. (1989).

Nebular abundances.

To compute the N++/O++ ratio, it requires to determine the volume emissivities of the [N iii] 57  μm and the [O iii] 88  μm lines. In this context, we use the nebular/IRAF package by assuming ne = 65 cm-3 and Te = 7000 K (see Eq. (4) in Vamvatira-Nakou et al. 2016). From the measured line fluxes and their uncertainties (Table 3), the N++/O++ abundance ratio is 1.06 ± 0.04. In order to see whether this ratio is representative of the overall N/O ratio, we determine the N++/N+ ratio from the [N iii] 57  μm and the [N ii] 122  μm lines. This ratio is equal to 2.9. As mentioned by Stock & Barlow (2014), when this value is larger than unity, this indicates that N++ is the dominant ionization state of nitrogen. Therefore, given the similarity of their ionization potential, we can assume that O++ is also the dominant ionization state for oxygen in the observed field. In conclusion, we can assume that the N++/O++ is representative of the global N/O ratio. This value is much higher than the solar value, estimated to be (N/O) ~ 0.14 (Grevesse et al. 2010).

An estimate of the nitrogen abundance in the H1 region can also be made. In this context, we use the Hα flux that we have determined, as well as the fluxes measured on the [N iii] 57  μm, the [N ii] 122  μm, and the [N ii] 205  μm lines. We assume a case-B recombination with Te = 7000 K (Draine 2011). The ionic abundances N+/H+ and N++/H+ have been determined with the nebular/IRAF package, and then summed to give the final value of N/H = 5.4 ± 1.4 × 10-4, by number. If we consider Te = 16 600 K (the upper limit given by Leitherer & Chavarria 1987) and ne = 70 cm-3, the N/H ratio is equal to 3.1 ± 1.2 × 10-4 by number.

Photodissociation region.

The nebular spectrum taken in the H1 region displays the [O i] 63  μm, the [O i] 146  μm and the [C ii] 158  μm lines that probably indicate the presence of a PDR around HD 148937. In PDRs, the ratio between the [O i] 63  μm and the [C ii] 158  μm lines is measured to be smaller than ten (Tielens & Hollenbach 1985), whilst in shock regions, it is larger. Given the ratio between these two lines, it is unlikely that their presence is due to shocks.

Neutral oxygen is found only in neutral zones. Therefore, the [O i] lines exclusively arise from the PDR (Malhotra et al. 2001) but neutral atomic carbon has an ionization potential lower than that of hydrogen and one thus expects to find C ii emission arising from both the H ii regions and the PDR.

The C/O abundance ratio can be estimated from the PDR line fluxes following a method described by Vamvatira-Nakou et al. (2013). This method was used to disentangle 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 (1)where we define with F[C ii] 158 being the total flux of the [C ii] 158  μm line, α a factor to be determined and ϵ the volume emissivity. Assuming that ⟨ C+ ⟩ / ⟨ N+ ⟩ = C/N and calculating the emissivities using the nebular/IRAF package with ne = 65 cm-3 and Te = 7000 K, we find (2)From N/O = 1.06, and from the line fluxes reported in Table 3, we obtain (3)using the measured F[C ii]/F[N ii] ratio. Assuming that there is a pressure equilibrium between the ionized gas region and the PDR, we have (4)that is used to define a locus of possible values in the diagram of Fig. 9. From these values (red lines in Fig. 9), we have , where by definition [C/O] = log (C/O)−log (C/O). Using the values of the line fluxes from Table 3 and (C/O) = 0.5 (Grevesse et al. 2010), we solve the equations to obtain α = 0.95 ± 0.03 and C/O = 1.42 ± 0.39. That means that almost all the flux of the [C ii] 158  μm line comes from ionized gas. We must emphasize, however, that the estimate of the flux for the [C ii] 158  μm line in the two areas observed by PACS must be taken with caution, given the high flux measured for this line in the off-source spectrum.

thumbnail Fig. 9

Temperature-density PDR diagnostic diagram. The grid of flux ratios F[O i] 63/F[O i] 146 versus was calculated by solving the level population equations for a range of temperatures and densities. The black and the red lines correspond to the H2 and the H1 regions, respectively. The quasi-vertical dashed lines indicate the pressure equilibrium constraint between the H ii region and the PDR, with the corresponding errors (dash-dot lines). The horizontal dashed line corresponds to the observed ratio with the errors (see Vamvatira-Nakou et al. 2013, for more details).

5.2.2. The H2 region

Electron density.

Among the different regions of the nebula studied by Leitherer & Chavarria (1987) and Dufour et al. (1988), three are close to the H2 area observed by Herschel/PACS. Assuming Te = 7000 K through the H2 region, as we did for the entire nebula, is in agreement with the measurements of Leitherer & Chavarria (1987) and Dufour et al. (1988) who estimated Te to be equal to 6800 K and 7600 K, respectively. We compute a ratio of 3.46 ± 0.19 between the [N ii] 122  μm and the [N ii] 205  μm lines. To determine the electron density inside this H2 region, we use the nebular package of the IRAF/STSDAS environment. From Te = 7000 K and the ratio of 3.46, we determine an electron density of ne = 125 ± 15 cm-3. This value is smaller than what we found from Spitzer data at position #7. The aperture size and a large gradient of electron density in this region of the nebula can explain such a difference.

Hα flux.

We use the Hα image and the same extinction value as reported in Sect. 5.2.1 to derive the Hα flux in the H2 region. We determine a reddened and a de-reddened flux of F0(Hα) = 3.32 × 10-11 erg cm-2 s-1 and F0(Hα) = 1.46 ± 0.6 × 10-10 erg cm-2 s-1, respectively. The values for the Hα flux in the H2 region are therefore larger than in the H1 region.

Nebular abundances.

As we did for the H1 region, we estimate the volume emissivities of the [N iii] 57  μm and the [O iii] 88  μm lines by assuming that ne = 125 cm-3 and Te = 7000 K. From these values and the fluxes reported in Table 3, we compute a N++/O++ ratio of 1.54 ± 0.06 for the H2 region. We also compute a N++/N+ ratio of 1.9, indicating that the N++/O++ is representative of the global N/O ratio. Therefore, we conclude that the N/O value is larger than that determined for the H1 region. Furthermore, we note that the N/O ratio estimated in the H2 region is similar to the value provided by Dufour et al. (1988) (N/O ~ 1.5, for their “A” and “B” regions), given the uncertainties.

The nitrogen content is also computed by assuming a case-B recombination with Te = 7000 K. The sum of the ionic abundances N+/H+ and N++/H+ gives a final value of N/H = 4.3 ± 1.3 × 10-4, by number. The N/H ratio seems similar in the H1 and the H2 regions, within the uncertainties. This means that the oxygen is responsible for the difference between the N/O ratios obtained in the H1 and H2 regions, confirming its possible depletion in the nebula (Dufour et al. 1988).

Photodissociation region.

As already detected in the H1 region, the [O i] 63  μm, [O i] 146  μm and [C ii] 158  μm lines are also visible in the H2 region. By following the same method, we determine the percentage of flux arising from the ionized region for the [C ii] 158  μm line as well as the C/O ratio.

In this context, we use the N/O ratio of 1.54 ± 0.06 that was measured in the H2 region. The temperature-density PDR diagnostic diagram (see the black lines in Fig. 9) provides us , obtained from the ratio between the [O i] 63  μm and the [O i] 146  μm lines of about 1.0 ± 0.1 (Table 3). By solving the equations as reported in Sect. 5.2.1, we obtain α = 0.88 ± 0.01 and C/O = 2.24 ± 0.10. The [C ii] lines measured from the H1 and the H2 regions are thus clearly dominated by the H ii region. Given the size of the Strömgren radius (~1.1°) surrounding the nebula NGC 6164/5, it therefore appears normal that the PDR is only a small fraction of this immense sphere. Although the C/O ratio must be taken with caution given the background spectrum, the C/O = 1.42 and 2.24 in the H1 and H2 regions, respectively, can also suggest an oxygen depletion that has been observed in the nebula (see, e.g. Dufour et al. 1988, for further details) and already mentioned before.

5.2.3. Mass of the ionized nebula

Vamvatira-Nakou et al. (2013) provided two equations to determine the ionized mass of a nebula: one using the Hα flux and one using the radio emission.

thumbnail Fig. 10

Best-fit CMFGEN spectrum of HD 148937 (red line) compared to FEROS spectrum taken on June 26th, 2005. The lines that are not modeled are diffuse interstellar bands (DIBs). Note the presence of an emission component in H and He i (narrow) lines, due to confined winds.

The Hα flux is measured from the Gemini Hα image, after having recalibrated it with the HST image (see Sect. 2.1.2). We measure a reddened flux of F0(Hα) = 5.15 ± 0.95 × 10-10 erg cm-2 s-1. This value is compatible, at 3σ, with the Hα flux estimated by Frew et al. (2013). After having corrected it for the extinction, we obtain a de-reddened Hα flux for the whole nebula of F0(Hα) = 2.12 ± 0.39 × 10-9 erg cm-2 s-1. Equations (B.7) and (B.16) provided by Vamvatira-Nakou et al. (2013) are dependent on the fraction of the volume of the nebula that is filled by ionized gas (the ϵ parameter) as well as the fraction of He in the form of He+. The former parameter is unknown for NGC 6164/5 whilst the latter has been estimated to 0.14 by Leitherer & Chavarria (1987). Considering Te = 7000 K, an angular radius of the nebula equal to 185″, and a distance of 1.3 kpc (Herbst & Havlen 1977), we estimate the mass of the ionized gas to be M. Therefore, for an arbitrary relatively small value of the filling factor, ϵ = 0.01, the ionized gas mass is 1.6 M. A more conservative estimation would be to consider that all the ionized material is located in the lobe regions. In this case, we consider a shell as first approximation of the 3D ionized nebula with an outer radius of 190″ (~1.2 pc at 1.3 kpc) and an inner radius of 142″ (~0.9 pc at 1.3 kpc). Under this assumption, we estimate a volume of 4.4 pc3 for the shell and of 7.4 pc3 for the sphere taking the whole structure into account, giving a fraction of the volume of the ionized gas of 0.6. This value would suggest a mass of the ionized gas of 12.4 M.

The radio emission used to determine the mass of the ionized gas is given by Milne & Aller (1982). These authors reported an integrated flux of 2.50 Jy at ν = 5 GHz and of 2.54 Jy at ν = 14.7 GHz. We thus estimate the mass of the ionized gas to be M with the same assumptions as before. Both values are in agreement within the errors on the different parameters.

To summarize, the exact determination of the ionized mass is difficult given the uncertainties on several parameters. We can certainly assume that the mass ejected by HD 148937 is larger than 1.6 M, which agrees with the results of Leitherer & Chavarria (1987).

6. Modeling of the central star

To perform the modeling of the central star, we employ the non-LTE CMFGEN atmosphere code (Hillier & Miller 1998). This code simultaneously solves the radiative transfer equation for a spherically symmetric wind in the co-moving frame and the rate equations under the constraints of radiative and statistical equilibrium. The velocity structure is constructed from a pseudo-photospheric structure connected to a β-velocity law of the form v = v(1−R/r)β where v is the terminal velocity. The density structure is computed from mass conservation. The photospheric structure is obtained after a few iterations of the hydrodynamical solution in which the radiative force computed from the level populations and atomic data is included.

Our final model includes H i, He i-ii, C ii-iv, N ii-v, O ii-vii, Ne ii-iv, Si iii-iv, Mg ii, S iii-v, Ar iii-v, Ca iii-iv, Fe ii-vii, and Ni ii-vi with the solar composition of Grevesse et al. (2010) unless otherwise stated. CMFGEN also uses the super-level approach to reduce the memory requirements. On average, we include about 1870 super levels for a total of 8050 levels. The final synthetic spectrum is obtained from a formal solution of the radiative transfer equation. A microturbulent velocity varying linearly with velocity from 10 km s-1 to 0.1 × v was used. We include X-ray emission in the wind since this can affect the ionization balance and the strength of key UV diagnostic lines. In practice, we adopt a temperature of three million degrees and we adjust the flux level so that the X-ray flux coming out of the atmosphere matches the observed LX/Lbol ratio equal to 10-6.1 (Nazé et al. 2008b).

The final spectrum is then convolved by a rotational profile estimated by the Fourier transform method (Simón-Díaz & Herrero 2007), and then convolved again by a Gaussian profile to mimick an isotropic macroturbulence. A description of the main diagnostic lines is provided by Mahy et al. (2015).

The best-fit CMFGEN spectrum (Fig. 10) gives an effective temperature of 40 000 ± 2000 K, and log g of 4.0 ± 0.1. The luminosity of HD 148937 is taken from Wade et al. (2012), that is, . The helium lines are pretty well reproduced as well as the other absorption lines except for the cores of the Balmer lines that seem affected by emission. Only the N iv 4058 Å and the C iv 5801–12 Å doublet are poorly fitted. The surface carbon abundance, by number, is estimated to 2.1 ± 0.5 × 10-4, the surface nitrogen abundance to 4.5 ± 0.8 × 10-4 and the surface oxygen abundance to 1.7 ± 0.3 × 10-4. HD 148937 is thus depleted in carbon and oxygen and overabundant in nitrogen. The N/O and C/O ratios determined for HD 148937 are 2.65 ± 1.14 and 1.24 ± 0.62, respectively. N/O is higher than in the nebula, as expected qualitatively. Our stellar parameters as well as our abundances agree with the values of Martins et al. (2015). The wind parameters are determined from the UV lines, the He ii 4686Å and Hα lines. We are, however, not able to fit the UV and the Hα lines with a single set of wind parameters (, v, clumping filling factor, and β). A possible cause can be the spectral variability of the object, but it is most probable that, because of the magnetic confinement, the spherical symmetry assumption, implicit to CMFGEN, is no longer valid. The fundamental parameters (stellar, wind, and abundances) of HD 148937 are listed in Table 4.

Table 4

Fundamental parameters of HD 148937.

7. Discussion

7.1. Evolutionary tracks

HD 148937 is a magnetic star with a 7.03-day variability that is associated to its rotational period. Wade et al. (2012) determined veq/vcrit = 0.12 and an obliquity of the magnetic field equal to 38°.

In the present study, we confirm the effective temperature and the gravity for HD 148937. Combined to the known luminosity of the star, these parameters provide a radius of 16.6 R, giving HD 148937 a rotational velocity of veq ≲ 120 km s-1 when accounting for the period of 7.03 days. The analysis of the optical spectrum of HD 148937 infers a projected rotational velocity of vsini ≲ 22 km s-1. We thus obtain an inclination of the stellar rotation axis of 11°, which is consistent, within the errors, with i ≤ 30° derived from the magnetic field analysis (Wade et al. 2012).

thumbnail Fig. 11

Current location of HD 148937 (red cross) in the HR diagram and locations at the time of nebular ejections, estimated from the N/O ratios of the H1 region in blue and in the H2 region in green. Evolutionary tracks in black are from Ekström et al. (2012) whilst the red ones are from Brott et al. (2011).

The stellar parameters of HD 148937 allow the determination of its position in the Hertzsprung-Russell (HR) diagram (Fig. 11). We compare this position with the evolutionary tracks of Brott et al. (2011) and of Ekström et al. (2012). The evolutionary tracks of Brott et al. (2011) take the transport of angular momentum by magnetism into account through the Tayler-Spruit dynamo formalism (Spruit 2002; Petrovic et al. 2005) but do not consider possible transport of chemical elements as a result of this dynamo process, whilst the tracks of Ekström et al. (2012) do not account for magnetism, either by a dynamo mechanism or by any strong fossil field.

The presence of a magnetic field can, however, reduce the rotational rate of the star as it was observationally detected by Townsend et al. (2010) in the case of σ Ori E. As models do not take this effect into account, we favor evolutionary tracks with low vsini to constrain the current status of HD 148937. In this context, the tracks of Brott et al. (2011), computed with an initial rotational velocity of 100 km s-1, yield an initial mass of 58 M to HD 148937. This value agrees with the initial mass of 57 M determined from the tracks of Ekström et al. (2012) computed with an initial rotational rate of 0.1.

The tracks of Brott et al. (2011) reproduce the value of veq/vcrit = 0.12 at the location of HD 148937 in the HR diagram. However, the N/O ratio that we find from the atmosphere models does not match that derived from the evolutionary tracks for a star with Teff = 40 000 K and (i.e., N/O = 0.15). In the same way, the evolutionary tracks of Ekström et al. (2012) fail at reproducing the N/O ratio, inferring values slightly too small in comparison to the observations, as well as a value too small for the current veq/vcrit. It is, however, worth noting that the inclusion of magnetism, under the assumption that the Tayler-Spruit dynamo is working, in the evolutionary models predicts an increase of the surface velocity and an increase of the mixing (Meynet & Maeder 2005), implying higher surface abundances for a given initial rotational rate. An agreement can then be found between the N/O ratios derived for the nebula and the stellar value by considering tracks with an initial rotational velocity of 300 km s-1. For the tracks of Ekström et al. (2012), the star would have an initial mass of 50 M and would have an age of 3.8 Myr; the material located in the H1 region with N/O = 1.06 would have been ejected by a 2.6-Myr star whilst the ejecta situated in the H2 region (where N/O = 1.54) would be produced by a 3.1-Myr star. For the tracks of Brott et al. (2011), computed with an initial rotational velocity of 300 km s-1, the star would have a mass of 57 M and currently have an age of 2.1 Myr, the H1 region would have been ejected when the star was 0.8 Myr old and the H2 region when the star was 1.6 Myr old. Therefore, the ejecta of the H1 region would have been ejected 1.2–1.3 Myr ago and those of the H2 region about 0.6 Myr ago. This suggests that the H2 region is younger than the H1 region even though it is further located in projection in the nebula. It thus turns out that the ejecta seem to have a complex morphology/kinematics. These values also infer an expansion velocity of about 2 km s-1which is not realistic for such an object if we consider that the material travels straight toward the lobes. All these values should however be considered as preliminary as all magnetic processes are not yet (fully) implemented in the models.

7.2. Evolutionary scenarios

Up to now, HD 148937 is the only magnetic O-type star known to be surrounded by a nebula. The exact formation process of this nebula raises many questions. The fact that HD 148937 is, so far, the most massive Galactic star with a detected magnetic field could play a role in the presence of this circumstellar nebula even though its magnetic field is not the most powerful detected among the massive star population. Moreover, no evidence exists that HD 148937 is a binary system. Neither Nazé et al. (2008b) from their multiwavelength spectroscopic survey, nor Sana et al. (2014) from their interferometric observations have detected the presence of a companion.

Therefore, two possible scenarios can be considered to explain the ejection of such a nebula: a giant eruption triggered by the stellar wind and the magnetic field, or a merger event between two massive stars in a binary configuration.

7.2.1. Giant eruption scenario

The ejecta constituting the nebula is known to be enriched, and has thus a stellar origin (see Leitherer & Chavarria 1987; Dufour et al. 1988; and the present study). We have roughly estimated a lower limit of 2 M for the mass of the ejecta. García-Segura et al. (1999) showed that the formation of planetary nebulae with such bipolar shapes is possible around magnetic stars through interactions of two succeeding time-independent stellar winds. Therefore, as a common origin for shaping bipolar nebulae, the critical rotation during a phase of strong mass loss could explain the presence of such a nebula. In this scenario, the stellar material would be first ejected in the equatorial plane. A second faster wind would then be ejected in the same plane, and would collide with the slower and denser wind and would then collimate along the polar direction. The nebula would thus be composed of two parts: an equatorial disk and a bipolar structure oriented perpendicularly.

To see whether this scenario can be applied here, we must analyze the orientation of HD 148937. Let us assume that the bipolar structure forming the nebula NGC 6164/5 is oriented in the polar direction. In this case, the equatorial ejection would correspond to the inner parts detected in the different images of Fig. 2, close to the central star, whilst the lobes would be created by the collimated material. Under this assumption, the poles/lobes would be younger than the equatorial ejecta, which is what we determine from the analysis of the abundances. However, such a configuration, where the lobes are in the polar direction, cannot be reconciled with the low inclination of about 15° suggested from the 7.03-day rotational period and from the nebula (see Sect. 3). It also suggests that the central star was close to its critical rotational velocity and should remain high even though the ejecta have taken out some momentum, which is also incompatible with the observations.

Let us assume now that the bright lobes are in the equatorial plane and that the inner zones of ejecta constitute the polar structure of the nebula. This agrees better with the low inclination but, in this case, the chemical abundances determined for the lobes should be less enriched than in the inner parts of the nebula, contrary to what is observed. Furthermore, the RVs of the ejecta, close by projection to the central star, should be approaching the observer and have thus a negative value, which is also not confirmed by the different studies of the kinematics.

thumbnail Fig. 12

Sketch of the helical morphology in three dimensions as viewed from the side. An Earth-based observer is at the left.

The scenarios involving a giant eruption associated with wind-wind interactions thus fails at putting together all the pieces of the puzzle and at characterizing the global properties (morphology, abundances, and kinematics) of the nebula. Furthermore, by assuming that the timescales of the ejections, determined from the evolutionary models, are correct, the expansion velocity of the nebula, estimated to 2 km s-1, would be too small to justify a giant eruption. As an alternative, let us consider the morphology presented by Carranza & Agüero (1986) proposed from their kinematic analysis of the nebula. The nebula would have been ejected through a giant eruption in the equatorial plane of HD 148937. Because of the misalignment between the magnetic and the rotational axes, the material would then follow a spiral path toward the poles of the magnetic axis, giving to the nebula a helical morphology (Fig. 12). The material located close to the central star by projection would have been ejected first, and thus would be less enriched than the ejecta found in the lobes, as observed. This also agrees with the low inclination of the rotational axis of HD 148937. This scenario would turn, however, the magnetic wind channeling backwards (from the equatorial plane toward the poles) rather than from the poles toward the equational plane, as in the usual confinement scenario (Babel & Montmerle 1997; ud-Doula et al. 2009, and references therein).

7.2.2. Binary merger scenario

The binary merger scenario was introduced by Langer (2012) to explain the nature of the nebula around HD 148937. In this interpretation, this author considered HD 148937 as a binary merger, that may shed material in the circumstellar medium either through a bipolar nebula or through a circumbinary disk (see de Mink et al. 2014). The ages of the lobes suggested by the evolutionary tracks tend to support a less violent process, such as a merging, even though these values must still be confirmed. The merging process could even generate a magnetic field (Tout et al. 2008) but the main issue in this scenario is the high rotational rate that would be obtained by merging. This rate would indeed be much higher than what has been derived for HD 148937 (veq ≤ 120 km s-1). One may, however, wonder at what rate the magnetic field would brake the rotation in such a disturbed event. If this rate is very high, it could have slowed down the star quickly, but this scenario needs to be simulated to assess its validity. Another question would be to know the chemical repartition in such an event and see whether it is possible to find enriched material far from the central merger, as it is observed for HD 148937. All these questions must still be investigated before such a scenario can be seen as physically possible and coherent considering the parameters and abundances of the star and of its ejecta obtained in the present analysis.

8. Conclusions

We have presented the analysis of the Herschel/PACS imaging and spectroscopic data of the nebula NGC 6164/5 surrounding the Of?p star HD 148937 together with Hα and mid-infrared images, as well as high-resolution optical spectra. The Hα image shows a bipolar or “8”-shaped ionized nebula whilst the infrared images describe a dust nebula smaller and closer to the central star. These images also exhibit a cavity in the dust material around HD 148937 that is supposed to be created by the strong stellar wind.

The far-infrared spectra taken in a H1 region close to the star and in a H2 region located further away, in the brightest part of NGC 6164, were analyzed. They indicate the presence of ionized material with N/O ratios of 1.06 and 1.54, respectively. This result, combined with the studies of Leitherer & Chavarria (1987) and Dufour et al. (1988), shows that the ejecta forming the brightest lobes in the Hα image are more enriched than the material located close to HD 148937. By coupling the analysis of the abundances with the kinematic analysis of the nebula, notably presented by Carranza & Agüero (1986), the most probable scenario capable of explaining the formation process of NGC 6164/5 would be a giant eruption ejecting material in the equatorial plane. This material, thanks to the presence of the magnetic field, would move towards the magnetic poles through spiral motions. This scenario would suggest two different inclinations, one for the rotational axis and another one for the magnetic axis. The helical motion would thus be triggered by a cumulative effect of the stellar wind and of the magnetic field of HD 148937. A stellar merger event can, however, not be excluded, but detailed modeling is required to see whether theoretical predictions match observations.

Acknowledgments

This research was supported by the PRODEX XMM and Herschel contract (Belspo), the Fonds National de la Recherche Scientifique (F.N.R.S.) and through the ARC grant for Concerted Research Actions, financed by the French Community of Belgium (Wallonia-Brussels Federation). We also thank Pr. G. Rauw and Dr. F. Martins for their helpful comments and Pr. D. J. Hillier for making his code CMFGEN available. We are grateful to the staffs of La Silla ESO Observatory and of Gemini Observatory for their technical support.

References

  1. Babel, J., & Montmerle, T. 1997, A&A, 323, 121 [NASA ADS] [Google Scholar]
  2. Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245 [NASA ADS] [CrossRef] [Google Scholar]
  4. Carranza, E., & Agüero, E. L. 1986, Ap&SS, 123, 59 [NASA ADS] [CrossRef] [Google Scholar]
  5. Catchpole, R. M., & Feast, M. W. 1970, The Observatory, 90, 136 [NASA ADS] [Google Scholar]
  6. de Mink, S. E., Sana, H., Langer, N., Izzard, R. G., & Schneider, F. R. N. 2014, ApJ, 782, 7 [NASA ADS] [CrossRef] [Google Scholar]
  7. Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton University Press) [Google Scholar]
  8. Dufour, R. J., Parker, R. A. R., & Henize, K. G. 1988, ApJ, 327, 859 [NASA ADS] [CrossRef] [Google Scholar]
  9. Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Frew, D. J., Bojičić, I. S., & Parker, Q. A. 2013, MNRAS, 431, 2 [NASA ADS] [CrossRef] [Google Scholar]
  11. García-Segura, G., Langer, N., Różyczka, M., & Franco, J. 1999, ApJ, 517, 767 [NASA ADS] [CrossRef] [Google Scholar]
  12. Grevesse, N., Asplund, M., Sauval, A. J., & Scott, P. 2010, Ap&SS, 328, 179 [NASA ADS] [CrossRef] [Google Scholar]
  13. Groenewegen, M. A. T., Waelkens, C., Barlow, M. J., et al. 2011, A&A, 526, A162 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Herbst, W., & Havlen, R. J. 1977, A&AS, 30, 279 [NASA ADS] [Google Scholar]
  15. Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef] [Google Scholar]
  16. Houck, J. R., Charmandaris, V., & Brandl, B. R. 2005, in The Initial Mass Function 50 Years Later, eds. E. Corbelli, F. Palla, & H. Zinnecker, Astrophys. Space Sci. Lib., 327, 527 [Google Scholar]
  17. Hubrig, S., Schöller, M., Schnerr, R. S., et al. 2008, A&A, 490, 793 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Langer, N. 2012, ARA&A, 50, 107 [NASA ADS] [CrossRef] [Google Scholar]
  19. Lebouteiller, V., Barry, D. J., Spoon, H. W. W., et al. 2011, ApJS, 196, 8 [NASA ADS] [CrossRef] [Google Scholar]
  20. Lebouteiller, V., Barry, D. J., Goes, C., et al. 2015, ApJS, 218, 21 [NASA ADS] [CrossRef] [Google Scholar]
  21. Leitherer, C., & Chavarria, C. 1987, A&A, 175, 208 [NASA ADS] [Google Scholar]
  22. Mahy, L., Rauw, G., De Becker, M., Eenens, P., & Flores, C. A. 2015, A&A, 577, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Maíz-Apellániz, J., Walborn, N. R., Galué, H. Á., & Wei, L. H. 2004, ApJS, 151, 103 [NASA ADS] [CrossRef] [Google Scholar]
  24. Malhotra, S., Kaufman, M. J., Hollenbach, D., et al. 2001, ApJ, 561, 766 [Google Scholar]
  25. Martins, F., & Plez, B. 2006, A&A, 457, 637 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Martins, F., Escolano, C., Wade, G. A., et al. 2012, A&A, 538, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Martins, F., Hervé, A., Bouret, J.-C., et al. 2015, A&A, 575, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Meynet, G., & Maeder, A. 2005, A&A, 429, 581 [Google Scholar]
  29. Milne, D. K., & Aller, L. H. 1982, A&AS, 50, 209 [NASA ADS] [Google Scholar]
  30. Nazé, Y., Walborn, N. R., & Martins, F. 2008a, Rev. Mexicana Astron. Astrofis., 44, 331 [Google Scholar]
  31. Nazé, Y., Walborn, N. R., Rauw, G., et al. 2008b, AJ, 135, 1946 [NASA ADS] [CrossRef] [Google Scholar]
  32. Nazé, Y., ud-Doula, A., Spano, M., et al. 2010, A&A, 520, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Nazé, Y., Zhekov, S. A., & Walborn, N. R. 2012, ApJ, 746, 142 [NASA ADS] [CrossRef] [Google Scholar]
  34. Nazé, Y., Petit, V., Rinbrand, M., et al. 2014, ApJS, 215, 10 [NASA ADS] [CrossRef] [Google Scholar]
  35. Ott, S. 2010, in Astronomical Data Analysis Software and Systems XIX, eds. Y. Mizumoto, K.-I. Morita, & M. Ohishi, ASP Conf. Ser., 434, 139 [Google Scholar]
  36. Parker, Q. A., Phillipps, S., Pierce, M. J., et al. 2005, MNRAS, 362, 689 [NASA ADS] [CrossRef] [Google Scholar]
  37. Peretto, N., & Fuller, G. A. 2009, A&A, 505, 405 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Petrovic, J., Langer, N., & van der Hucht, K. A. 2005, A&A, 435, 1013 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1 [CrossRef] [EDP Sciences] [Google Scholar]
  40. Pismis, P. 1974, Rev. Mexicana Astron. Astrofis., 1, 45 [NASA ADS] [Google Scholar]
  41. Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  42. Roussel, H. 2013, PASP, 125, 1126 [Google Scholar]
  43. Sana, H., Le Bouquin, J.-B., Lacour, S., et al. 2014, ApJS, 215, 15 [NASA ADS] [CrossRef] [Google Scholar]
  44. Scowen, P., Hester, J., Code, A., et al. 1995, ApJ, 450, 196 [NASA ADS] [CrossRef] [Google Scholar]
  45. Simón-Díaz, S., & Herrero, A. 2007, A&A, 468, 1063 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Spruit, H. C. 2002, A&A, 381, 923 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Stock, D. J., & Barlow, M. J. 2014, MNRAS, 441, 3065 [NASA ADS] [CrossRef] [Google Scholar]
  48. Tielens, A. G. G. M., & Hollenbach, D. 1985, ApJ, 291, 722 [NASA ADS] [CrossRef] [Google Scholar]
  49. Tout, C. A., Wickramasinghe, D. T., Liebert, J., Ferrario, L., & Pringle, J. E. 2008, MNRAS, 387, 897 [NASA ADS] [CrossRef] [Google Scholar]
  50. Townsend, R. H. D., Oksala, M. E., Cohen, D. H., Owocki, S. P., & ud-Doula, A. 2010, ApJ, 714, L318 [NASA ADS] [CrossRef] [Google Scholar]
  51. ud-Doula, A., Owocki, S. P., & Townsend, R. H. D. 2009, MNRAS, 392, 1022 [NASA ADS] [CrossRef] [Google Scholar]
  52. Vamvatira-Nakou, C., Hutsemékers, D., Royer, P., et al. 2013, A&A, 557, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Vamvatira-Nakou, C., Hutsemékers, D., Royer, P., et al. 2016, A&A, 588, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. van Marle, A. J., Meliani, Z., Keppens, R., & Decin, L. 2011, ApJ, 734, L26 [NASA ADS] [CrossRef] [Google Scholar]
  55. Wade, G. A., Grunhut, J., Gräfener, G., et al. 2012, MNRAS, 419, 2459 [NASA ADS] [CrossRef] [Google Scholar]
  56. Walborn, N. R. 1972, AJ, 77, 312 [NASA ADS] [CrossRef] [Google Scholar]
  57. Werner, M. W. 2005, Adv. Space Res., 36, 1048 [NASA ADS] [CrossRef] [Google Scholar]
  58. Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Emission line fluxes for each spaxel of the H1 region

Table A.1

Line fluxes in each spaxel.

Appendix B: Emission line fluxes for each spaxel of the H2 region

Table B.1

Line fluxes in each spaxel.

Appendix C: Emission line fluxes for each spaxel of the off-source region

Table C.1

Line fluxes in each spaxel.

All Tables

Table 1

Line fluxes measured on the Spitzer IRS spectra.

Table 2

Electron temperatures, electron densities, and abundances across the nebula.

Table 3

Line fluxes measured on the integrated spectra of the H1 and H2 regions of NGC 6164/5, after subtracting the off-source values.

Table 4

Fundamental parameters of HD 148937.

Table A.1

Line fluxes in each spaxel.

Table B.1

Line fluxes in each spaxel.

Table C.1

Line fluxes in each spaxel.

All Figures

thumbnail Fig. 1

Top: general picture of the nebula NGC 6164/5 and definition of the nomenclature used in the present paper. North is up and east is left. The small white circles represent the positions of the Spitzer IRS spectra whilst the green squares give the footprints of the Herschel/PACS field of view. The NW and SE lobes are indicated by gray arrows and the inner zones of ejecta are represented in red. Bottom: Hα image of NGC 6164/5 taken from the SuperCosmos Hα Survey (Parker et al. 2005). The white boxes represent the areas of the nebula analyzed by Leitherer & Chavarria (1987) whilst the red arrows represent the positions analyzed by Dufour et al. (1988).

In the text
thumbnail Fig. 2

Images of NGC 6164/5 around HD 148937. Top: the Gemini Hα + [N ii] image (left), the WISE 12  μm and 22  μm images (middle and right, respectively). Bottom: PACS images of the nebula at 70  μm, 100  μm and 160  μm from left to right, respectively. The white asterisk indicates the position of HD 148937. The dash-lined box provides the size of the Hα image because the sizes of these images are different. North is up and east is left.

In the text
thumbnail Fig. 3

Three-color image of PACS photometer field of view with 70  μm in blue, 100  μm in green and 160  μm in red. North is up and east is left. The box indicates the location of the Hα image (see Fig. 2). The dashed line represents the location of the wind-blown structure at 12′ from HD 148937 (seen on the Hα image).

In the text
thumbnail Fig. 4

Spitzer IRS spectrum taken at position #2 (see Fig. 1). The step at 19  μm is due to the different apertures used for the SH and LH spectra.

In the text
thumbnail Fig. 5

PACS spectral lines of the H1 region in the nebula around HD 148937, integrated over the 25 spaxels and background corrected, using the off-source region. The continuum has been removed for clarity.

In the text
thumbnail Fig. 6

PACS spectral lines of the H1 region (in solid line) and of the off-source region (in dashed line) integrated over the 25 spaxels. The continuum has been removed for clarity.

In the text
thumbnail Fig. 7

As for Fig. 5 but for the H2 region.

In the text
thumbnail Fig. 8

As for Fig. 6 but for the H2 region.

In the text
thumbnail Fig. 9

Temperature-density PDR diagnostic diagram. The grid of flux ratios F[O i] 63/F[O i] 146 versus was calculated by solving the level population equations for a range of temperatures and densities. The black and the red lines correspond to the H2 and the H1 regions, respectively. The quasi-vertical dashed lines indicate the pressure equilibrium constraint between the H ii region and the PDR, with the corresponding errors (dash-dot lines). The horizontal dashed line corresponds to the observed ratio with the errors (see Vamvatira-Nakou et al. 2013, for more details).

In the text
thumbnail Fig. 10

Best-fit CMFGEN spectrum of HD 148937 (red line) compared to FEROS spectrum taken on June 26th, 2005. The lines that are not modeled are diffuse interstellar bands (DIBs). Note the presence of an emission component in H and He i (narrow) lines, due to confined winds.

In the text
thumbnail Fig. 11

Current location of HD 148937 (red cross) in the HR diagram and locations at the time of nebular ejections, estimated from the N/O ratios of the H1 region in blue and in the H2 region in green. Evolutionary tracks in black are from Ekström et al. (2012) whilst the red ones are from Brott et al. (2011).

In the text
thumbnail Fig. 12

Sketch of the helical morphology in three dimensions as viewed from the side. An Earth-based observer is at the left.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.