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

HD 148937 is a peculiar massive star (Of?p) with a strong magnetic field (1 kG). The hourglass-shaped emission nebula NGC 6164/5 surrounds this star. This nebula is presumed to originate from episodic mass-loss events of the central O-type star, but the detailed formation mechanism is not yet well understood. Grasping its three-dimensional structure is essential to uncovering the origin of this nebula. Here we report the high-resolution multiobject spectroscopic observations of NGC 6164/5 using the GIRAFFE on the 8.2 m Very Large Telescope. Integrated intensity maps constructed from several spectral lines delineate well the overall shape of this nebula, such as the two bright lobes and the inner gas region. The position–velocity diagrams show that the two bright lobes are found to be redshifted and blueshifted, respectively, while the inner region has multiple layers. We consider a geometric model composed of a bilateral outflow harboring nitrogen-enriched knots and expanding inner shells. Its spectral features are then simulated by using a Monte Carlo radiative transfer technique for different sets of velocities. Some position–velocity diagrams from simulations are very similar to the observed ones. According to the model that best reproduces the observational data, the two bright lobes and the nitrogen-enriched knots are moving away from HD 148937 at about 120 km s−1. Their minimum kinematic age is estimated to be about 7500 yr. We discuss possible formation mechanisms of this nebula in the context of binary interaction.


INTRODUCTION
Massive stars are not only the sources of enormous ultraviolet 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 lifetime is 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(Dale et al. , 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 to understand 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 Lim et al.
those of HD 148937 have the shortest period recorded up to now, only 7 d (Nazé et al. 2008).As 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 towards 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 d 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 ∼ 40000K, log g ∼ 4.0, v sin i ≥ 45 km s −1 - Nazé et al. 2008Nazé et al. , 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 symmetrically with respect to the star.He then proposed the two features to belong to one single elongated structure of about 5 ′ diameter, at the time identified as a planetary nebula.In addition, a series of arcs delineate a larger shell of about 12 ′ radius centered on the star and even larger nebular features can be spotted at a distance of 64 ′ from the star (Westerlund 1961).The largest structure is interpreted as a Strömgren sphere linked to the intense ultraviolet (UV) radiation of the massive star (Fairall et al. 1985), the broken-arc shell as a bubble blown by the strong stellar wind of the massive star(s), and the NGC 6164/5 pair as recent ejecta (Leitherer & Chavarria-K. 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 understand 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 NW of the star, displays positive radial velocities (RVs) while NGC 6165, to the SE 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 −1 in NGC 6164 and from -170 to -3 km s −1 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-K. 1987).Some previous studies suggested that the main lobes may be part of either a spiral structure (Mahy et al. 2017), a bipolar structure (Pismis 1974;Leitherer & Chavarria-K. 1987), or a helical structure (Carranza & Agüero 1986).
The abundances of NGC 6164/5 were also studied in detail.Danks & Manfroid (1977) already 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-K. 1987).These non-solar 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 counter-intuitive 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 a N/O ratio about one dex above the solar abundance while the C/O ratio appeared three times the solar one.Larger values were found for the outer NW 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 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-K. 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 see also review by Langer 2012).However, it may be difficult to reach a rotation as slow as that of HD148937 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 blue-shifted velocities, which is not observed (Pismis 1974;Carranza & Agüero 1986;Leitherer & Chavarria-K. 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) has 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, a more in-depth investigation, with modern means, of  (Reid et al. 1991;Lasker 1994;Djorgovski et al. 1998).This optical image can be found in MAST: http://dx.doi.org/10.17909/T9QP4J.The fiber position numbers are labeled on the image.The positions of the fibers are relative to HD 148937 (R.A. = 16 h 33 m 52.s 387, decl.= −48 • 06 ′ 40.′′ 477).
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 threedimensional (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 positionvelocity (PV) maps in Section 3. We simulate a model based on the 3D distribution of the nebula in positionvelocity space and present a comparison of the observational data with the simulated one in Section 4. The formation of this nebula is discussed in Section 5, and our results are summarized in Section 6.

OBSERVATIONS AND DATA
The spectroscopic observations of NGC 6164/5 were performed with the European Southern Observatory 8.2m 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 highresolution spectrograph GIRAFFE1 , 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 1195s.Figure 1 displays the positions of the FLAMES fibers .The orderseparating filter HR15N was used to cover the spectral range of 6470 Å to 6790 Å.Some dome flat and ThAr lamp spectra were also obtained after the target observation for calibration.
The observational data were reduced by following the standard reduction procedures under the EsoReflex environment (Freudling et al. 2013).We extracted onedimensional 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.
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.

Integrated Intensity Maps
Integrated intensity maps from different emission lines are useful tools to investigate 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.
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 of 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.) with respect to RVs.Some emission components are centered in a narrow RV range between −18 km s −1 and −32 km s −1 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 ranges 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.

PV diagrams
We investigated the velocity field of NGC 6164/65 in position-velocity (PV) space.Our multi-object 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(Lim et al. , 2019(Lim et al. , 2021)).We then   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 −1 (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 −1 and −27 km s −1 , respectively, i.e. they correspond to a symmetrically expanding feature.
Hα and He I λ6678 lines trace the distribution of hot gas filling the space between the two lobes, while [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 Hα and He I λ6678 lines ranges from −75 km s −1 to 50 km s −1 , equivalent to −52 km s −1 and 73 km s −1 relative to HD 148937, respectively.
The spatial resolution of our multi-object spectroscopy is about 0. ′ 5×0.′ 5 for the inner 6. ′ 0×6.′ 0 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 bestfit multiple Gaussian distributions.
Figure 6 displays the 3D PV diagrams.The RV distribution of 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 −1 , 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 a axisymmetric 3D structure.
The  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-K. 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 −1 relative to the central star.Such high-RVs may be close to 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 −1 ) as shown in Figure 7.Given the fact that the RVs of knots in the main lobes are about 30 km s −1 relative to the central star with an uncertainty of several km s −1 , the inclination angle can then be constrained in a range of 15 • to 22 • (see the dashed line in the figure ).In what follows, we thus adopted an inclination angle of 20 • .This is larger than that the 5 • angle inferred by Leitherer & Chavarria-K. (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 sky and the equator of the O-type star.

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 right cartoon of Figure 8).In this study, θ LOS was set to 20 • .
We considered that NGC 6164/5 is composed of three components: nitrogen-enriched knots (orange), extended envelopes (red) surrounding the nitrogenenriched knots and expanding thin (hollow) shells.The thickness of the knots and envelopes are characterized by an outer radius R o and inner radius R i = 0.7R o .Both opening angles of the knots and envelopes, θ 1 and θ 2 , were set to 15 • .The hollow shells occupy an inner region characterized by the inner and outer radii decreasing as a function of cos θ z , where θ z is measured from z-axis and ranges from θ 1 + θ 2 to 75 • .Therefore, the inner radii decrease from 0.4R o cos θ z to ∼ 0.10R o (close to the central star), and the outer ones decrease 0.5R o cos θ z to ∼ 0.13R o .
The knots display N [II] λ6584/Hα ratios higher than those of the envelopes.Their similar position 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 knots and envelopes were set to increase with the distance from the central star and reach v 1 (knots) and v 2 (envelopes) at R o .This study considered the situation that the two components are comoving (v 1 = v 2 ).Two different setups for their velocities at R o , 100 km s −1 and 120 km s −1 as shown in Figure 7, were applied to our simulations.On the other hand, the hollow shells are inflated from z-axis at a constant velocity v 3 .Three different setups for v 3 (80, 100, and 120 km s −1 ) were considered in our simulations.

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) 2 .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 6 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 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 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 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 (3) The Doppler factor ∆V of the collected photon, which is equivalent to the RV, is given by where v kine is the velocity at the location emitting the photon, which is computed by v 1 , v 2 , and v 3 .
In our simulation, we applied the projection effects to the the structure and velocities of the modelled 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) .(Reid et al. 1991;Lasker 1994;Djorgovski et al. 1998, right panel).The left panel shows the intensity map in the projected radius, considering the weight factors of 3, 1, and 0.1 for the nitrogen-enriched knots, envelopes, and hollow shells, respectively.Thus, the bright regions to the north and south are composed primarily of photons emitted from nitrogen-enriched knots.The color bar represents the level of normalized intensity.The Cartesian coordinates xp and yp correspond to R.A. (different sign) and Dec., respectively.The right panel displays the optical image for comparison with our model.The color map highlights the structure of NGC 6164/5.

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 comparison of the integrated intensity map of our model with the observational image from 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 aim, 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 −1 because multiple components along the line of sight include both backward-and forward-moving gas.This is reminiscent of what is found in observation (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 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 −1 as found in observations.On the other hand, the inner hollow shells cover different RVs, depending on v 3 = 80, 100, and 120 km s −1 .
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 ′ × 6 ′ ).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 RV were corrected by the stellar RV (−23 km s −1 ).The models adopting low velocities (v 1 = v 2 = 100 km s −1 , v 3 = 80 or 100 km s −1 ) tend to underestimate RVs of the nebula.On the other hand, the PV diagrams from the models assuming highvelocities agree best with the observed ones, particularly for v 1 = v 2 = v 3 = 120 km s −1 . 5. DISCUSSION Leitherer & Chavarria-K. (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 19th century, η Car has ejected a large amount of mass (10-35M ⊙ ) at velocities of about 650 km s −1 (Smith & Ferland 2007;Mehner et al. 2016).Some ejecta of the Homunculus nebula are moving at velocities of several hundred to several thousand km s −1 (Weis et al. 2001;Smith 2008).NGC 6164/5 may not originate from such a η 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 maximum 100 km s −1 .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 WR124 is composed of a bipolar structure and a torus (Zavala et al. 2022).The bipolar component has an expanding velocity of 88 km s −1 (Sirianni et al. 1998), while the velocities of the torus range from 30 km s −1 to 60 km s −1 (Zavala et al. 2022).It is often thought that WR 124 does not have a companion, however, several observational studies do not definitely rule out this possibility (Toalá et al. 2018;Moffat et al. 1982;Gvaramadze et al. 2010).Zavala et al. (2022) explained the formation of M 1-67 and WR124 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 years.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 to interact 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 2. ′ 7 away from HD 148937.Since the distance to the central star is about 1.1 kpc according to the offset-corrected Gaia parallax (Gaia Collaboration et al. 2023;Lindegren et al. 2021), 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 −1 , the travel time of the knots from the central star to their current positions is estimated to be about 7,500 years 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 7,500 years.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 have at least 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 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 the result of the deformation of the stellar surfaces by tidal forces on the primary star.While the tidal force can strongly affect the stellar surface 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.

SUMMARY
The magnetic O-type star HD 148937 is surrounded by the emission nebula NGC 6164/5.In this study, we 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 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 morpho-kinematics 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 −1 .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 further 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 −1 .This yields a kinematic age of 7500 years.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.

Figure 2 .
Figure2.One-dimensional spectrum of a given position on NGC 6164/5.In the upper panel, the red dashed line traces the continuum level from the third-degree polynomial fitting, while the lower panel shows the continuum-subtracted spectrum.Nebular spectral lines are labeled in the upper panel.

Figure 3 .
Figure 3. Integrated intensity maps from four spectral lines: Hα (first panel), [N II] λ6584 (second panel), He I λ6678 (third panel), and [S II] λ6731 (fourth panel).The signals in the RV range from −200 km s −1 to 200 km s −1 are integrated at a given fiber position.The integrated intensities are shown in a logarithmic scale.The boxes outlined by dashed lines represent a 6 ′ × 6 ′ region from HD 148937.These maps are smoothed by using a Gaussian smoothing technique, where the kernel sizes of 0. ′ 25 and 0. ′ 5 are adopted for the inner (6 ′ × 6 ′ ) and outer regions, respectively.A clear difference can be observed between permitted (Hα and He I λ6678) and forbidden ([N II] λ6584 and [S II] λ6731) lines.

Figure 4 .
Figure 4. Distributions of log[N II] λ6584/Hα (left panel) and log[S II] λ6730 / [S II] λ6716 (right panel) with respect to RVs.The integrated flux of individual lines were obtained from the best-fit Gaussian distributions.
constructed the PV diagrams by the sum of intensities along R.A. and decl., respectively.
[N II] λ6584/Hα line ratios of the individual components are shown by different sizes of 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 envelope of the two lobes and the inner hot gas tend to have lower [N II] λ6584/Hα line ratios.

Figure 5 .
Figure 5. Position-velocity diagrams obtained from Hα (upper panel), [N II] λ6584 (middle panel), and He I λ6678 (lower panel).The blue sphere represents the central O-type star HD 148937.The contour shows the integrated intensities in ADU.The coordinates are relative to HD 148937 as shown in Figure 1

Figure 6 .
Figure 6.PV diagrams (left panels) and [N II] λ6584/ Hα line ratios (right panels) of each gas component.In the left panels, the RVs of individual gas components are shown by different colors.The symbol sizes and colors are proportional to the line ratios [N II] λ6584/ Hα in the right panels.The star symbols represent HD 148937.

Figure 7 .
Figure7.Expected velocities of ejecta with respect to different inclination angles.We considered ejection velocities of 80 (dotted line), 100 (dashed line), and 120 (solid line) km s −1 .The observed RV, measured on the main lobes, is shown by a dashed line.

Figure 8 .
Figure 8. Schematic illustration of our bilateral model in a cartesian coordinate.The x-axis corresponds to the rotational axis, and the equator of the central star is placed along the y-z plane.The left cartoon displays the structure and velocity vectors.The nebula has three components, nitrogen-enriched knots (orange), envelopes (red), and expanding hollow shells (blue).The velocities of the knots and envelopes (v1 and v2) are proportional to the distance from the central star (green dots).The hollow shells are also expanding from z-axis at v3.The right cartoon shows the inclination of this system from the line of sight or the plane of sky.The inclination angle (θLOS) of 20 • is adopted in this model (see the main text for detail).
of sight k defined by kp = (sin θ p cos ϕ p , sin θ p sin ϕ p , cos θ p ), (1) where θ p and ϕ p are the azimuthal and polar angles of the line of sight.The basis vectors of the projected plane xp and ŷp are given by xp = (sin ϕ p , − cos ϕ p , 0), (2) ŷp = (− cos ϕ p cos θ p , − sin ϕ p cos θ p , sin θ p ), ϕ p is fixed at 0 • because the geometry is symmetric to z−axis.θ 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.

Figure 9 .
Figure 9. Integrated intensity map of our bilateral model (left panel) and the Digitized Sky Survey image(Reid et al. 1991;Lasker 1994; Djorgovski et al. 1998, right panel).The left panel shows the intensity map in the projected radius, considering the weight factors of 3, 1, and 0.1 for the nitrogen-enriched knots, envelopes, and hollow shells, respectively.Thus, the bright regions to the north and south are composed primarily of photons emitted from nitrogen-enriched knots.The color bar represents the level of normalized intensity.The Cartesian coordinates xp and yp correspond to R.A. (different sign) and Dec., respectively.The right panel displays the optical image for comparison with our model.The color map highlights the structure of NGC 6164/5.

Figure 10 .Figure 11 .
Figure10.Maps of mean RVs (top panels) and RVs of random spots (bottom panels) for different setups, calculated from the Doppler shifts of the expanding substructures.In the bottom panels, nitrogen-enriched knots, envelopes, and hollow shells are shown by triangle, circle, and plus symbols, respectively.