A Morphokinematic Study of the Enigmatic Emission Nebula NGC 6164/5 Surrounding the Magnetic O-type Star HD 148937

Massive stars are not only the sources of enormous ultraviolet (UV) radiation in the universe, but also the factories of heavy elements as they evolve. In addition, they are good tracers of star-forming regions in the galaxies (Garcia et al. 2009) because their lifetimes are only about several million years. In such regions, star formation may be regulated by feedback from massive stars in both constructive (Elmegreen & Lada 1977) and destructive (Dale et al. 2012, 2013) ways. At the end of stellar evolution, massive stars leave compact objects that can be associated with high-energy transient events such as γ-ray bursts. Hence, many fields of astronomy require an understanding of the formation and evolution of massive stars. In particular, mass loss is a crucial parameter for understanding the evolution of massive stars (Smith 2014).

HD 148937 is a massive star of spectral type Of?p (Walborn 1972). One specific feature of this spectral class is the presence of periodic line profile variations and those of HD 148937 have the shortest period recorded up to now, only 7 days (Nazé et al. 2008). Like all other stars of this type, HD 148937 has been found to be strongly magnetic, with a dipolar field of polar intensity 1 kG (Wade et al. 2012). The spectral peculiarities are then explained by an oblique magnetic rotator scenario (Babel & Montmerle 1997). In this model, the stellar winds are channeled toward the equator by the strong field, where they collide. This confined material emits throughout the electromagnetic spectrum, with an enhanced and harder X-ray emission as well as numerous emission lines in the optical range. Because of the obliquity of the magnetic axis, the confined winds are seen under different angles throughout the rotation of the star, leading to the observed variations: the 7 days period of HD 148937 is thus interpreted as the rotation period of the star. However, the low amplitude of the variations, compared to other magnetic O stars, suggested a low inclination for the star (Nazé et al. 2010), which was then confirmed by the magnetic monitoring (Wade et al. 2012): HD 148937 is thus seen nearly pole-on. In addition, interferometric and spectroscopic measurements recently indicated that HD 148937 actually is a binary with two similarly massive components: a highly eccentric orbit and a long orbital period (Wade et al. 2019). Note that the inclination angle of the orbital plane of the system is also low. Finally, the stellar parameters of HD 148937 were well constrained by several previous studies (mass ∼50–60M_⊙; effective temperature ∼40,000 K; $\mathrm{log}g\sim 4.0$; and $v\sin i\,\geqslant 45$ km s⁻¹—Nazé et al. 2008, 2010; Wade et al. 2012; Martins et al. 2015). It is also worth noting that HD 148937 shows a nitrogen enrichment as high as seen in O-type supergiants (Martins et al. 2015).

HD 148937 is surrounded by several layers of nebulosities. Henize (1959) was the first to recognize that the nebulae NGC 6164 and NGC 6165 actually appear symmetrical with respect to the star. He then proposed the two features to belong to one single elongated structure of about $5^{\prime}$ diameter, at the time identified as a planetary nebula. In addition, a series of arcs delineate a larger shell of about $12^{\prime}$ radius centered on the star and even larger nebular features can be spotted at a distance of $64^{\prime}$ from the star (Westerlund 1961). The largest structure is interpreted as a Strömgren sphere linked to the intense UV radiation of the massive star (Fairall et al. 1985), the broken-arc shell being a bubble blown by the strong stellar wind of the massive star(s) and the NGC 6164/5 pair being recent ejecta (Leitherer & Chavarria 1987). While ejected material is commonly spotted around evolved massive stars (of Wolf–Rayet or luminous blue variable type), NGC 6164/5 is a rare example of nebulosities associated to a "normal," unevolved massive star, making its study crucial to better understanding the mechanisms of stellar feedback.

The kinematics of the NGC 6164/5 nebula has been studied several times in the past. First, the overall difference between the two parts was recognized: NGC 6164, to the northwest of the star, displays positive radial velocities (RVs), while NGC 6165, to the southeast of the star, shows negative RVs (Catchpole & Feast 1970). This is characteristic of material expanding from a center, so it was taken as an additional argument in favor of the NGC 6164/5 unification scheme. More detailed investigations, however, revealed the complexities of the kinematic structure, with velocities ranging from +117 to –70 km s⁻¹ in NGC 6164 and from –170 to –3 km s⁻¹ in NGC 6165 (Pismis 1974; Carranza & Agüero 1986). Even nearby areas could show very different velocities (Pismis 1974). In fact, the nebular line profiles often reveal several components, especially in the inner parts of the nebula, demonstrating the presence of several layers along the line of sight (Leitherer & Chavarria 1987). Some previous studies have suggested that the main lobes may be part of a spiral structure (Mahy et al. 2017), a bipolar structure (Pismis 1974; Leitherer & Chavarria 1987), or a helical structure (Carranza & Agüero 1986).

The abundances of NGC 6164/5 were also studied in detail. Danks & Manfroid (1977) have mentioned an enhanced He abundance, while the detailed study of Dufour et al. (1988) found an enrichment in N and He as well as a depletion in Ne and O (see also Leitherer & Chavarria 1987). These nonsolar abundances constitute an additional argument in favor of a nebula formed by an ejection event from HD 148937. It should be noted in this context that the star itself displays a clear nitrogen enrichment, as would be expected for such a scenario. However, in a counterintuitive way, the enrichment was found to be larger in the farther lobes than in the inner parts of the nebula (Dufour et al. 1988). This was confirmed in the Herschel data analysis of Mahy et al. (2017). The inner parts close to the star (structure "H1" in that paper) displayed an N/O ratio about 1 dex above the solar abundance, while the C/O ratio appeared three times the solar one. Larger values were found for the outer northwest lobe (structure "H2"), which also appeared denser. Assuming the nebular material was ejected by the star, its abundances would then reflect those at the surface of the star at the time of the ejection. Compared with the predictions of single-star evolution models, the observed nebular abundances then suggest an ejection about 0.6 Myr ago for the lobes and about 1 Myr ago for the inner parts.

Note that the extinction was often found to be rather uniform across the nebula (e.g., Leitherer & Chavarria 1987), although a dark cloud has been reported in the area (Peretto & Fuller 2009). Interpreting the various nebular features in a coherent scheme has proven to be difficult. Assuming a single event, one could think of ejection following a merger. This scenario is linked to the fact that HD 148937 is magnetic and to the theoretical proposal that magnetic fields could be generated in merging events (Tout et al. 2008; Ferrario et al. 2009; and see also the review by Langer 2012). However, it may be difficult to reach a rotation as slow as that of HD 148937 so soon after the merging. Besides, the detection of a companion by Wade et al. (2019) potentially implies that the history of the system is more complex, and it remains to be demonstrated that the initial triple system can keep the distant companion in the currently observed orbit after the merger.

Because of the bipolar appearance of the nebula, its axis, assumed to be similar to the stellar (rotation) or orbital axis, was often thought to be close to the plane of the sky, i.e., the star would be seen equator-on. In this case, the inner parts would correspond to slower material ejected first in the equatorial plane, while the distant lobes would correspond to a second event having faster ejecta collimated along the polar axis. This configuration would be similar to what is seen in the famous η Carinae homunculus (Smith 2008). However, the inclinations of the rotation and orbital axes were clearly determined to be low, i.e., the system is rather seen pole-on. In such a case, one could instead consider the inner parts to come from a polar ejection and the lobes from a subsequent equatorial ejection. It may be noted that Danks & Manfroid (1977) had already proposed two different directions of ejection for the inner parts and the lobes, with a possible pole-on configuration. However, the inner parts should then present the largest red- or blueshifted velocities, which is not observed (Pismis 1974; Carranza & Agüero 1986; Leitherer & Chavarria 1987). Besides, the multiple velocity components are not easily reconciled with a simple ejection event as envisaged in such scenarios.

Based on the above arguments, a more complex configuration needs to be anticipated. Carranza & Agüero (1986) have notably proposed a helical geometry, where the north and south ejection flows could get into the same line of sight, as well as different parts of the same flows. However, it is unclear how the material would flow along the polar axis away from the star. Indeed, magnetic channeling works in the other direction. Definitely, then, a more in-depth investigation, with modern means, of NGC 6164/5 is warranted, with the hope of clarifying the exact geometry of the ejected flows.

We perform the spectroscopic study of NGC 6164/5 with a high spectral resolution and a dense fiber configuration. The aim of this study is to uncover the three-dimensional (3D) geometry of the nebula and investigate its origin in the context of massive star evolution. In Section 2, the observational data we used are presented. The observational features of NGC 6164/5 are investigated using integrated intensity maps and position–velocity (PV) maps in Section 3. We simulate a model based on the 3D distribution of the nebula in PV space and present a comparison of the observational data with the simulated data in Section 4. The formation of this nebula is discussed in Section 5, and our results are summarized in Section 6.

The spectroscopic observations of NGC 6164/5 were performed with the European Southern Observatory 8.2 m Very Large Telescope (UT2) on Cerro Paranal in Chile on 2019 March 24. We took a total of 220 spectra over the nebula using the intermediate-to-high-resolution spectrograph GIRAFFE, ⁹ whose fibers are fed by FLAMES (Pasquini et al. 2002). A total of three frames were obtained for the same field. The exposure time of each frame was set to 1195 s. Figure 1 displays the positions of the FLAMES fibers. The order-separating filter HR15N was used to cover the spectral range of 6470–6790 Å. Some dome flat and ThAr lamp spectra were also obtained after the target observation for calibration.

Zoom In Zoom Out Reset image size

The observational data were reduced by following the standard reduction procedures under the EsoReflex environment (Freudling et al. 2013). We extracted one-dimensional (1D) spectra using the IRAF/SPECRED package. Figure 2 exhibits one of the extracted spectra. Polynomial fitting was performed on the continua of all spectra (see the red dashed line in the upper panel), and the continua were subtracted by the best-fit curves to obtain usable normalized spectra. Our spectra show the usual nebular emission lines in this domain, e.g., [N ii] λ λ 6548, 6584, Hα, He i λ6678, and [S ii] λ λ6719, 6731.

Zoom In Zoom Out Reset image size

We fit the emission lines Hα, [N ii] λ6584, and [S ii] λ λ 6716, 6730 with Gaussian profiles using the MIDAS command DEBLEND/LINE. Since most lines have several velocity components along the line of sight, multiple Gaussian profiles were used to fit such complex line profiles. The RVs were taken from the centers of the best-fit Gaussians and then converted to heliocentric RVs. We also obtained the full width at half maximum and fluxes of the lines.

3.1. Integrated Intensity Maps

Integrated intensity maps from different emission lines are useful tools for investigating the spatial distribution of the gas. We created integrated intensity maps using a Gaussian smoothing technique for the positions and counts (ADU) from given fiber positions. Figure 3 displays the integrated intensity maps from the four different emission lines Hα, [N ii] λ6584, He i λ6678, and [S ii] λ6731.

Zoom In Zoom Out Reset image size

The integrated intensity map of Hα line reproduces well the overall structure seen in optical images. The diffuse gas is detected mostly over NGC 6164/5. The region extending out of the nebula (outside the white box delineated by a dashed line in Figure 3) may be part of the filamentary halo heated by shocks and photoionization from the central O-type star (Fairall et al. 1985). The northern and southern lobes contain bright knots. The He i λ6678 line traces a similar structure for the nebula. However, this line was not detected in the halo.

The [N ii] λ6584 map highlights the presence of the northern and southern lobes, while the line strength of the diffuse gas is not as strong as those seen in the Hα map. The [S ii] λ6731 map traces almost the same features as those in the [N ii] λ6584 map, and the halo gas was also detected. Given the fact that the formation of these two forbidden lines is sensitive to electron density, the two main lobes have different electron densities (or density structures) from the other regions.

The spectral lines of NGC 6164/5 are affected by emission from the filamentary halo (Fairall et al. 1985), as seen in the integrated intensity map for Hα lines. Figure 4 displays the line ratios (log[Nii] λ6584/Hα and log[S ii] λ6730/[S ii]) with respect to RVs. Some emission components are centered in a narrow RV range between −18 and −32 km s⁻¹ around the system velocity of HD 148937, and they show moderate spreads in the two line ratios (−1.3 to 0.0 for log[N ii] λ6584/Hα and −0.4 to 0.3 for log[S ii] λ6730/ [S ii] λ6716). Most of these emissions may originate from the filamentary halo. Meanwhile, the other components associated with NGC 6164/5 show a much broader range of RVs and of line ratios. In order to probe the kinematics of NGC 6164/5, it is necessary to minimize the contribution of the halo component. The emission components associated with the halo were discarded from further analysis.

Zoom In Zoom Out Reset image size

3.2. PV Diagrams

We investigated the velocity field of NGC 6164/65 in PV space. Our multiobject spectroscopic data contain four-dimensional information, R.A., decl., velocity, and intensity. However, the fibers were not placed in a regular grid. A Delauney triangulation technique was applied to interpolate the intensities and velocities of the emission lines (Hα, [N ii] λ6584, and He i λ 6678) into data cubes composed of 90 × 90 × 90 regular volume cells (Lim et al. 2018, 2019, 2021). We then constructed the PV diagrams by the sum of intensities along R.A. and decl., respectively.

Figure 5 displays the PV diagrams of NGC 6164/5. We also plotted a dot to mark the position of the central O star HD 148937 in the diagrams. This star is known to be a binary (Wade et al. 2019) and its systemic RV is about −23 km s⁻¹ (Wade et al. 2019). In the PV diagrams, there are two prominent gas components with respect to HD 148937, corresponding to the northern and southern lobes. The knots in the northern and southern lobes are moving away from the central star at RVs of about 33 km s⁻¹ and −27 km s⁻¹, respectively, i.e., they correspond to a symmetrically expanding feature.

Zoom In Zoom Out Reset image size

The Hα and He i λ6678 lines trace the distribution of hot gas filling the space between the two lobes, while the [N ii] λ6584 line is too weak to probe the velocity field of the inner hot gas. The RVs of the inner hot gas from the Hα and He i λ6678 lines range from −75 to 50 km s⁻¹, equivalent to −52 and 73 km s⁻¹ relative to HD 148937, respectively.

The spatial resolution of our multiobject spectroscopy is about 05 × 05 for the inner 60 × 60 region, which is higher than in previous spectroscopic observations. Despite this, some small structures can be blurred, due to the limited resolution and uneven grid of fibers. In particular, velocity components with weak line intensities may be more affected by such effects in the construction of PV diagrams. We thus investigated the individual velocity components decomposed from the best-fit multiple Gaussian distributions.

Figure 6 displays the 3D PV diagrams. The RV distributions of the individual gas components are consistent with that shown in Figure 5. Also, we confirmed that there are high-velocity gas components between the two main lobes. Their velocities relative to the central O-type star exceed 100 km s⁻¹, which is faster than the velocities of the two main lobes. The PV diagrams in Figures 5 and 6 show an almost symmetric velocity field, indicating that the nebula surrounding HD 148937 may have an axisymmetric 3D structure.

Zoom In Zoom Out Reset image size

The [N ii] λ6584/Hα line ratios of the individual components are shown by the different sizes of the dots in the right panels of Figure 6. The majority of points with high [N ii] λ6584/Hα line ratios (red dots) are found in the knots of each lobe, although some are close to the central region, as seen in Figure 3. The envelopes of the two lobes and the inner hot gas tend to have lower [N ii] λ6584/Hα line ratios.

4.1. Geometry

The morphology of NGC 6164/5 has been explained by three different geometric models invoking a helical structure (Carranza & Agüero 1986), a spiral structure (Mahy et al. 2017), and a bipolar structure (Pismis 1974; Leitherer & Chavarria 1987). In this study, we revisited the geometric model of the bipolar structure. Since HD 148937 has a pole-on configuration (Nazé et al. 2010; Wade et al. 2012), we will consider that NGC 6164/5 was probably ejected from the equator of the star. Hereafter, the bipolar structure is thus referred to as the bilateral structure.

If the main lobes from the plane of the sky were not inclined from the plane of the sky, its RVs would not be detected by spectroscopy. This inclination angle is one important parameter in setting up a geometric model. In order to infer the inclination angle, it is necessary to constrain the velocity of the nebula in the actual (physical) 3D space, from the known RVs. The diagrams shown in Figures 5 and 6 indicate that some parts of the Hα-emitting nebula have RVs exceeding 100 km s⁻¹ relative to the central star. Such high RVs may be close to the actual velocities if the associated nebular region is moving nearly parallel to the line of sight, but may be very different from the actual ones if inclined on the line of sight.

We thus computed the expected RVs of the main lobes with respect to different inclination angles for three ejection velocities (80, 100, and 120 km s⁻¹), as shown in Figure 7. Given the fact that the RVs of the knots in the main lobes are about 30 km s⁻¹ relative to the central star with an uncertainty of several km s⁻¹, the inclination angle can then be constrained in a range of 15°–22° (see the dashed line in the figure). In what follows, we thus adopt an inclination angle of 20°. This is larger than that the 5° angle inferred by Leitherer & Chavarria (1987), but it remains compatible with the rotational inclination angle (≤30°) measured from the line of sight (Wade et al. 2012), i.e., the angle between the plane of the sky and the equator of the O-type star.

Zoom In Zoom Out Reset image size

4.2. Setup

Figure 8 displays the schematic illustration of our bilateral model. The x-axis corresponds to the rotational axis of HD 148937, and the y–z plane is thus an equatorial plane. The rotational axis is inclined from the line of sight by θ_LOS (see the right cartoon of Figure 8). In this study, θ_LOS was set to 20°.

Zoom In Zoom Out Reset image size

We considered that NGC 6164/5 is composed of three components: nitrogen-enriched knots (orange), extended envelopes (red) surrounding the nitrogen-enriched knots, and expanding thin (hollow) shells. The thickness of the knots and envelopes is characterized by an outer radius R_o and an inner radius R_i = 0.7R_o . Both opening angles of the knots and envelopes, θ₁ and θ₂, were set to 15°. The hollow shells occupy an inner region characterized by the inner and outer radii decreasing as a function of $\cos {\theta }_{z}$, where θ_z is measured from the z-axis and ranges from θ₁ + θ₂ to 75°. Therefore, the inner radii decrease from $0.4{R}_{o}\cos {\theta }_{z}$ to ∼0.10R_o (close to the central star), and the outer ones decrease from $0.5{R}_{o}\cos {\theta }_{z}$ to ∼0.13R_o .

The knots display N [ii] λ6584/Hα ratios higher than those of the envelopes. Their similar positions may indicate that the knots have almost the same expanding velocities as the envelopes or slightly higher velocities. We therefore postulated that the knots and envelopes may have different origins.

The velocities of the knots and envelopes were set to increase with the distance from the central star and reach v₁ (knots) and v₂ (envelopes) at R_o . This study considered the situation that the two components are comoving (v₁ = v₂). Two different setups for their velocities at R_o , 100 and 120 km s⁻¹, as shown in Figure 7, were applied to our simulations. On the other hand, the hollow shells are inflated from the z-axis at a constant velocity v₃. Three different setups for v₃ (80, 100, and 120 km s⁻¹) were considered in our simulations.

4.3. Simulation

We developed a new 3D Monte Carlo radiative transfer simulation to investigate our spectral cube data, based on the Monte Carlo techniques used in Chang & Lee (2020). ¹⁰ In this simulation, we set the kinematics and geometry of the emission nebula as described in Section 4.2. In order to compare our model with observational data, our simulations uniformly generate a total of 10⁶ photon packets within the model nebulae. Each photon packet carries a velocity v_kine and a position r _i = (x_i , y_i , z_i ) at the emitting location. Furthermore, these components have different intensities for Hα and [N ii] λ6584 lines (see Figure 3), and therefore the weight factors of the photon packets, representing the relative emissivity in a unit volume, of the nitrogen-enriched knots, envelopes, and hollow shells were set to 3, 1, and 0.1, respectively. After that, we collected photons along the line of sight. Since the scope of our study is to explore the kinematic structure of NGC 6164/5 in comparison with the observed data, detailed physical parameters, such as electron densities and the number of photons from the central star, were not considered in these simulations.

The model nebula in 3D coordinates is projected onto the plane corresponding to the plane of the sky. The number of photons at any given point was summed along the line of sight to construct the integrated intensity map. Given that the unit vector of the line of sight $\hat{{\boldsymbol{k}}}$ is defined by

$$\begin{eqnarray}&&\hat{{{\boldsymbol{k}}}_{p}}=(\sin {\theta }_{p}\cos {\phi }_{p},\sin {\theta }_{p}\sin {\phi }_{p},\cos {\theta }_{p}),\end{eqnarray} \tag{ 1 }$$

where θ_p and ϕ_p are the azimuthal and polar angles of the line of sight, the basis vectors of the projected plane $\hat{{{\boldsymbol{x}}}_{p}}$ and $\hat{{{\boldsymbol{y}}}_{p}}$ are given by

$$\begin{eqnarray}\begin{array}{rcl}\hat{{{\boldsymbol{x}}}_{p}} & = & (\sin {\phi }_{p},-\cos {\phi }_{p},0),\\ \hat{{{\boldsymbol{y}}}_{p}} & = & (-\cos {\phi }_{p}\cos {\theta }_{p},-\sin {\phi }_{p}\cos {\theta }_{p},\sin {\theta }_{p}),\end{array}\end{eqnarray} \tag{ 2 }$$

where ϕ_p is fixed at 0° because the geometry is symmetric to the z-axis and θ_p is 20° because it corresponds to θ_LOS.

The positions and RVs in the projected plane are needed to construct the PV diagram as shown in Figure 6. The coordinates in the projected plane are newly defined using the projected coordinates x_p and y_p , which correspond to R.A. and decl., respectively, given by

$$\begin{eqnarray}&&{x}_{p}={{\boldsymbol{r}}}_{{\boldsymbol{i}}}\cdot \hat{{{\boldsymbol{x}}}_{p}},\quad {y}_{p}={{\boldsymbol{r}}}_{{\boldsymbol{i}}}\cdot \hat{{{\boldsymbol{y}}}_{p}}.\end{eqnarray} \tag{ 3 }$$

The Doppler factor ΔV of the collected photon, which is equivalent to the RV, is given by

$$\begin{eqnarray}&&{\rm{\Delta }}V=-{{\boldsymbol{v}}}_{\mathrm{kine}}\cdot \hat{{\boldsymbol{k}}},\end{eqnarray} \tag{ 4 }$$

where v _kine is the velocity at the location emitting the photon, which is computed by v₁, v₂, and v₃.

In our simulation, we applied the projection effects to the structure and velocities of the modeled nebula. The following section will discuss the comparison with the simulated and observational data via the integrated intensity map (Figure 9) and the PV diagram (Figure 10).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.4. Comparison with the Observational Data

The integrated intensity map was rotated by 45° clockwise in the x_p − y_p plane to reproduce the tilt of NGC 6164/5, as seen in Figure 1. Figure 9 shows a comparison of the integrated intensity map of our model with the observational image from the Digitized Sky Survey (Reid et al. 1991; Lasker 1994; Djorgovski et al. 1998). The observational image shows intensity variation across the nebula because of thin dust lanes. Such dark regions are often found in H ii regions and are not very bright in infrared passbands. Their presence can only be recognized with bright background illumination. In our simulations, we did not consider the presence of thin dust lanes in front of NGC 6164/5.

We first investigated the global velocity trends. To this end, we averaged along the line of sight the Doppler shifts of the photons (upper panels of Figure 10). The northern and southern lobes have, on average, positive and negative RVs, respectively. Meanwhile, the mean velocities of the hollow shells are nearly 0 km s⁻¹, because multiple components along the line of sight include both backward- and forward-moving gas. This is reminiscent of what is found in the observations (Figure 5).

In order to probe the local variations of the RVs, the Doppler shifts were calculated at about 200 arbitrary spots within the 3D structure. The lower panels of Figure 10 exhibit the positions of such spots projected onto the sky, with different symbols depending on their associated substructure (knots, lobes, or shell) and different colors depending on their RVs (red/blue for forward-/backward-moving gas). As can be seen, there are multiple components with different RVs along the line of sight in the inner hollow shells and envelopes, as in the observations.

Figure 11 shows the PV diagrams taken from the randomly sampled data in Figure 10. In the simulated PV diagrams, photon-emitting components can be traced, as shown by different symbols. In all simulations, the nitrogen-enriched knots are found at around ±30 km s⁻¹, as found in the observations. On the other hand, the inner hollow shells cover different RVs, depending on v₃ = 80, 100, and 120 km s⁻¹.

Zoom In Zoom Out Reset image size

We directly compared the simulated PV diagrams with the observed ones (black dots) in Figure 11. The model nebulae occupy a region covering a range from −1 to +1 for the x_p - and y_p -axes (arbitrary units). Our spectroscopic observation covers the overall structure of NGC 6164/5, which corresponds to a (true) size of ($6^{\prime} \times 6^{\prime}$). For comparison, the observed PV diagrams were scaled by dividing the position relative to the star by 3 in arcminutes. Compared to Figure 6, note also that the observed RVs were corrected by the stellar RV (−23 km s⁻¹). The models adopting low velocities (v₁ = v₂ = 100 km s⁻¹, v₃ = 80 or 100 km s⁻¹) tend to underestimate the RVs of the nebula. On the other hand, the PV diagrams from the models assuming high velocities agree best with the observed ones, particularly for v₁ = v₂ = v₃ = 120 km s⁻¹.

Leitherer & Chavarria (1987) considered that HD 148937 has an equator-on configuration, and they modeled the mass ejection along the rotational axis of HD 148937, which they considered as nearly perpendicular to the line of sight. However, both spectroscopic and magnetic monitoring observations support a pole-on configuration (Nazé et al. 2008; Wade et al. 2012). Based on these more recent observations, a new simulation of the geometry of NGC 6164/5 was performed in this study. Our simulation reproduced well the observed velocity field, particularly for the kinematics of the inner region. It is now useful to compare the properties of NGC 6164/5 with those of mass ejection events from other massive stars in order to understand the possible formation mechanism of the emission nebula.

In the nineteenth century, η Car has ejected a large amount of mass (10–35 M_⊙) at velocities of about 650 km s⁻¹ (Smith & Ferland 2007; Mehner et al. 2016). Some ejecta of the Homunculus nebula are moving at velocities of several hundred to several thousand kilometers per second (Weis et al. 2001; Smith 2008). NGC 6164/5 may not originate from such an η Car-like giant eruption in terms of mass and velocity. Indeed, HD 148937 has lost some 2 M_⊙ (Mahy et al. 2017) and the ionized gas is receding away at a maximum of 100 km s⁻¹. In addition, there is a significant difference between the Homunculus nebula and NGC 6164/5. The great eruption of η Car occurred along the polar axis.

The emission nebula M 1–67 surrounding the Wolf–Rayet star WR 124 is composed of a bipolar structure and a torus (Zavala et al. 2022). The bipolar component has an expanding velocity of 88 km s⁻¹ (Sirianni et al. 1998), while the velocities of the torus range from 30 to 60 km s⁻¹ (Zavala et al. 2022). It is often thought that WR 124 does not have a companion, but several observational studies do not definitely rule out this possibility (Moffat et al. 1982; Gvaramadze et al. 2010; Toalá et al. 2018). Zavala et al. (2022) explained the formation of M 1–67 and WR 124 in the context of the common envelope scenario (Paczyński 1967).

According to Wade et al. (2019), HD 148937 has a massive O-type companion in a very eccentric orbit (e = 0.75) with an orbital period of 26 yr (or e = 0.58 with an orbital period of 18 yr). In that study, the inclination of the orbital plane from the plane of the sky was deduced to be about 40°. If the plane where NGC 6164/5 is expanding is nearly parallel to the orbital plane, the origin of the nebula could be associated with a binary interaction. Unlike the M 1–67 system, NGC 6164/5 may not be formed through the common envelope scenario, because the eccentric orbit and the long orbital period prohibit the two stars from interacting closely as in a common envelope scenario. However, other binary interactions due to tidal force could possibly occur when the primary and secondary stars come close together at periastron. This periastron passage can therefore result in mass ejection events for highly eccentric binaries (Regös et al. 2005; Moreno et al. 2011). Such a process could also explain why ejection occurs in a specific direction of the equatorial plane, rather than over the full equator.

The knots in the main lobes are about 27 away from HD 148937. Since the distance to the central star is about 1.1 kpc, according to the offset-corrected Gaia parallax (Lindegren et al. 2021; Gaia Collaboration et al. 2023), the projected distance of the knots is about 0.86 pc. If we consider the inclination angle of 20°, the actual distance between the central star and the knots is about 0.92 pc. With a velocity of 120 km s⁻¹, the travel time of the knots from the central star to their current positions is estimated to be about 7500 yr, assuming that the knots have rapidly accelerated and reached their final velocity on a short timescale. If a gradual acceleration of the knot is considered, the kinematic age would be somewhat longer than 7500 yr. This kinematic age is much smaller than 0.6–1.3 Myr, the age estimated from the comparison of the N/O abundances with evolutionary models (Mahy et al. 2017).

Here, we have shown that we can reproduce the appearance and kinematics of the nebula by assuming that all gas components were ejected at almost the same velocity. If the nitrogen-enriched knots were ejected at the same time at slower velocities, these knots would not reach their current position. The overall morphology of NGC 6164/5 would then be different from the current one. Therefore, these knots should at least have velocities similar to the other nebular components with less nitrogen enrichment if ejected simultaneously.

The puzzling difference in nitrogen abundance between the different features may then be understood if their origins are different. For example, knots and other features could be associated with the primary and companion stars, respectively. Given the fact that the nitrogen abundance of HD 148937 is as high as those of O-type supergiants (Martins et al. 2015), the primary star (49 M_⊙—Wade et al. 2019) may be evolving into the supergiant stage, while its companion star (34 M_⊙) may still be on the main-sequence stage.

η Car and several eccentric binaries are often assumed to have undergone eruptions around periastron passages (van Genderen & Sterken 2007). The brightenings seen in light curves were interpreted as resulting from the deformation of the stellar surfaces by tidal forces on the primary star. While the tidal force can strongly affect the stellar surfaces of binary pairs, it is still not fully understood how the binary interaction can/could drive a mass ejection. Nevertheless, we propose the origin of NGC 6164/5 to be similar. In order to confirm this idea, it is necessary to monitor the light curve of HD 148937 on a long timescale and spectroscopically find the signature of other mass ejection events around periastron passages. In addition, theoretical approaches utilizing smoothed particle hydrodynamics will give us a better understanding of mass ejection processes in detail.

The magnetic O-type star HD 148937 is surrounded by the emission nebula NGC 6164/5. In this study, we have investigated the kinematic and morphological properties of the nebula with high-resolution spectroscopy.

We measured six strong emission lines ([N ii] λ λ 6548, 6584, Hα, He i λ6678, and [S ii] λ λ 6719, 6731) in our data. The integrated intensity maps constructed from the Hα and He i λ6678 lines show the overall structure of the nebula, while the two forbidden lines [N ii] λ6584 and [S ii] λ6731 only trace the bright lobes.

The morphokinematics of NGC 6164/5 was investigated in PV diagrams. It confirms that the bright lobes are receding away from HD 148937 at RVs of about 30 km s⁻¹. The bright knots in these two lobes show high [N ii] λ6584/Hα line ratios, indicative of nitrogen enrichment. In addition, some nitrogen-enriched spots were also detected closer to HD 148937.

On the other hand, the gas surrounding the main lobes and in the inner regions tends to have lower [N II] λ6584/Hα ratios. In those regions, multiple components at different RVs are detected, indicating the presence of multiple layers along the line of sight.

In order to interpret the observed PV diagrams, we constructed a geometric model to reproduce the spectroscopic results. This model assumes a hollow expanding shell close to the star and lobes farther away. The envelopes of the lobes host nitrogen-enriched knots. Overall, the modeled nebula is expanding along the equatorial plane, which is tilted by 20° from the plane of the sky. To match the data, all features should expand with ∼120 km s⁻¹. This yields a kinematic age of 7500 yr. A single ejection velocity thus seems sufficient to reproduce the observed features, but abundances differ between them. This could however be explained in a binary interaction model, with ejection occurring at periastron from both stars. Further observational and theoretical studies are required to understand the mass ejection process of this system in detail.

The authors thank the anonymous referee for constructive comments and suggestions. This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia archive website is https://archives.esac.esa.int/gaia. The Digitized Sky Surveys were produced at the Space Telescope Science Institute under U.S. Government grant NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope on Palomar Mountain and the UK Schmidt Telescope. The plates were processed into the present compressed digital form with the permission of these institutions. The Second Palomar Observatory Sky Survey (POSS-II) was made by the California Institute of Technology with funds from the National Science Foundation, the National Geographic Society, the Sloan Foundation, the Samuel Oschin Foundation, and the Eastman Kodak Corporation. The Oschin Schmidt Telescope is operated by the California Institute of Technology and Palomar Observatory. The UK Schmidt Telescope was operated by the Royal Observatory Edinburgh, with funding from the UK Science and Engineering Research Council (later the UK Particle Physics and Astronomy Research Council), until 1988 June, and thereafter by the Anglo-Australian Observatory. The blue plates of the southern Sky Atlas and its Equatorial Extension (together known as SERC-J), as well as the Equatorial Red (ER), and the Second Epoch [red] Survey (SES), were all taken with UK Schmidt. This work was supported by a National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT; grant No. 2022R1C1C2004102). Y.N. acknowledges support from the Fonds National de la Recherche Scientifique (Belgium), the European Space Agency (ESA), and the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Programme.

Facility: VLT (FLAMES/GIRAFFE) - .

Software: EsoReflex (Freudling et al. 2013), IRAF (Tody 1986, 1993), MIDAS (Banse et al. 1992), NumPy (Harris et al. 2020).