© The Authors 2023

1 Introduction

Planetary nebulae (PNe) are astronomical objects representative of one of the latest stages of stellar evolution of low-and intermediate-mass stars (0.8–8.0 M_⊙). This phase occurs between the asymptotic giant branch (AGB) and white dwarf (WD) stellar phases. PNe are most known for their complex morphological structures, which are believed to arise during a brief (~10²–10⁴ yr) strong mass loss experienced in the postAGB phase (see the imaging catalogues by López et al. 2012; Sabin et al. 2014; Parker et al. 2016). Although the simplest morphologies are well explained with the interactive stellarwind model and its generalisation (Kwok et al. 1978; Balick 1987, respectively), more recent studies (e.g. Soker & Harpaz 1992; García-Segura et al. 2014, 2016) demonstrate that more complex morphologies, such as bipolar PNe, can naturally arise from a binary interaction via a common-envelope evolution (CEE; see Jones & Boffin 2017, for a review). However, the complexity of the different morphologies observed in PNe also indicates internal dynamical activity forming different structural components in the nebula such as jets, collimated outflows, rings, or knots indicative of different physical processes that are not very well constrained. The accurate determination of physical and chemical parameters of such PNe, along with detailed studies of their morpho-kinematic nature, are useful to understand their formation history (e.g. Kn 26, IPHASX J193718.6+202102, NGC 2392, NGC 40, and NGC 6543 by Guerrero et al. 2013, 2021; Sabin et al. 2021; Rodríguez-González et al. 2022; Clairmont et al. 2022, respectively).

M 1–16 (α_J2000 = 07^h37^m18^s.9, δ_J2000 = −09°38′47.″9), which was first discovered by Minkowski (1946), has a diameter of ≃3″ according to a study of optical images by Acker et al. (1982). This measurement was confirmed by Kwok (1985) using radio interferometry. An expansion velocity of the ionised gas of ~10 km s⁻¹ was measured in the [O III] λ5007 emission line by Sabbadin et al. (1986). Schwarz (1992) found that M 1–16 is a transition object, often called proto-PN, by means of optical, infrared (IR), and millimetric wavelength observations. They also found that M 1–16 exhibits multiple high-velocity outflows with shock velocities of around ~300 km s⁻¹ (Corradi & Schwarz 1993) and that the object is surrounded by a massive envelope of molecular gas. Aspin et al. (1993) revealed that M 1–16 also shows an extended emission in the near-IR (NIR) bands, mainly produced by vibrational excited molecular hydrogen (H₂), which appears as symmetrically opposed bipolar lobes (as compared with the optical image) with an extension of ~25″. They argued that these extended emissions of H₂ are related to collimated ‘jets’ with a maximum velocity of 28 km s⁻¹ with respect to the central star (in the plane of the sky). Aspin et al. (1993) also found that the effective temperature, T_eff, of the central star of PN (CSPN) is ~35 000 K, as determined by the ratio of the He I 2.058 μm and Br-γ, which is consistent with a young PN that has recently emerged from the proto-PN phase and is starting to ionise its surrounding gas. Also, high-resolution mapping in the CO J = 2−1 line of M 1–16 revealed a molecular envelope with a size of 50″, surrounding the ionised central 3″ nebula (Huggins et al. 2000). Huggins et al. (2000) argued that collimated outflows disrupted the central ionised nebula making cavities, with imprinted CO velocities up to 30 km s⁻¹, along the main axis.

M 1–16 shows different structural components in the inner regions, close to the central star, and in the outer regions as previously studied by Acker et al. (1982), Schwarz (1992), Aspin et al. (1993), and Corradi & Schwarz (1993). Lorenzo et al. (2021) also found extended and strong wings in the CO J = 9−8 band, indicative of fast collimated outflows. Guerrero et al. (2020) studied the jet nature of M 1–16 and found that the distance of the jets, with respect to the CSPN, is >0.5 pc, which is not expected in PNe with high-velocity jets.

In this paper, we present a deep analysis of the morpho-kinematic structure and chemical composition of the PN M 1–16 which reveal a different evolutionary path from those proposed in previous works. The observational data are described in Sect. 2, which includes deep imaging and high-resolution spectroscopy, and the results are presented in Sect. 3. The proposed origin of M 1–16 is presented in Sect. 4 and is supported by our morpho-kinematic model as well as the deep chemical analysis. Finally, the conclusions are presented in Sect. 5.

2 Observational data

2.1 Narrow-band optical imaging

Narrow-band CCD images of M 1–16 were obtained on 2015 April 4 using the 2.1 m telescope at the Observatorio Astronómico Nacional at the Sierra de San Pedro Mártir (OAN-SPM), Mexico. Observations were carried out using the image mode of the Manchester Echelle Spectrograph (MES; Meaburn et al. 2003). The detector was a 2048 × 2048 pixel E2V CCD (13.5 μm) with a plate scale of 0.″351 pix⁻¹ (using a 2 × 2 pixel binning mode). The filters used to acquire the images were Hα+[N II] (λ_c = 6563 Å; Δλ = 90 Å) and [O III] (λ_c = 5007 Å; Δλ = 60 Å). Exposure times were 1800 s for both filters with a seeing of ≃2″.

In addition, narrow-band CCD direct images of M 1–16 were obtained on 2021 April 13 using the 0.84 m telescope (f/15) at the OAN-SPM. The CCD camera MEXMAN was used to obtain deep images of M 1–16 in the Hα (λ_c = 6565 Å, Δλ = 11 Å) and [N II] (λ_c = 6583 Å, Δλ = 10 Å) narrow bands. The detector was a 2024×2024 E2V CCD with a size of 15 μm pix⁻¹ (resulting in a nominal plate scale of 0.22″ pix⁻¹). Exposure times were 1200 s and 600 s in Hα and [N II], respectively, using a binning of 2×2 pixels. Seeing was around 1.9″ during the observations. The images were processed using standard IRAF¹ routines. The MEXMAN images provide a different view of M 1–16 in the Hα and [N II] emission lines separately, contrary to the images obtained with the MES instrument in which the observations are from the combined Hα+[N II] emission lines.

2.2 High-resolution optical spectroscopy

Long-slit high-resolution optical spectra were also acquired with the MES instrument. We set a slit width of 150 μm (1.9″) with a fixed length of 6.5′ Spectra were obtained during a five-day observational campaign that began on 2015 April 4. Observations were carried out using the 2 × 2 binning mode (plate scale of 0.351″ pix⁻¹). The spectra were centred around the Hα filter with Δλ = 90 Å to isolate the 87th order (0.1 Å pix⁻¹ spectral scale). The exposure time for all of the spectra was set to 1800 s. The seeing was ≃1″ as measured from a full width at half maximum (FWHM) fit taken from the field stars’ continuum along the spatial axis of the spectra.

All the spectra were processed with standard techniques for long-slit spectroscopy using IRAF. The wavelength calibration was performed using a ThAr arc lamp to an accuracy of ±1 km s⁻¹. The FWHM of the arc-lamp emission lines was measured to be ≃12 km s⁻¹.

2.3 Boller & Chivens long-slit low-resolution optical spectroscopy

Low-resolution optical spectra for the M 1–16 were acquired using the Boller & Chivens spectrograph mounted in the 2.1 m telescope at OAN-SPM, during one observing run on 2015 November 19. An E2V CCD with a 2048 × 2048 pixel array and plate scale of 1.8″ pix⁻¹ (1×2 binning mode) was used as a detector. A spectral resolution of ≃5.0Å (FWHM), as judged by the arc calibration lamp spectrum, was obtained by using the 400 lines mm⁻¹ grating in conjunction with the 2″ wide slit, covering a spectral range of 4100–7600 Å.

Different exposure times were used in order to enhance different structural components in M 1–16. We used the median of three 180 s exposure spectra in order to study the bright central part of the nebula and ten spectra (three of 300 s, three of 180 s, two of 1800 s, one of 900 s, and one of 60 s) were co-added, leading to a total of 6000 s, to study the faintest regions along the main axis.

The spectra were reduced using standard procedures for long-slit spectroscopy under IRAF software. The calibration was performed using the standard star Feige 34. The stellar component was extracted and separated from the surrounding nebular emission when possible.

3 Results

3.1 Morphology

Our deep images, Hα+[N II]λ6584 Å, Hα, [N II]λ6584 Å, and [O ΙII]λ5007 Å, obtained with MES in image mode and with the MEXMAN camera, show the complex seemingly multi-polar morphology of M 1–16. Initially, we identified two pairs of lobes or outflows coming out from the central bright region of M 1–16 (Fig. 1). One pair of large lobes (LL), previously discovered by Schwarz (1992) and Aspin et al. (1993), is seen with a position angle (PA) on the plane of the sky of ≃ −31°, and they are labelled in the Fig. 1 as LL-SE and LL-NW (the SE and NW lobe, respectively). We measured their angular extension to be ≃47″ and ≃51″, with respect to the bright central region or main body (MB), for the LL-SE and LL-NW, respectively.

A similar value of 90″ for the full extent of the LL was measured by Schwarz (1992). A pair of Middle-sized Lobes (ML) are seen close to the MB, whose size in the sky is ≃24″ and ≃15″ for the ML located to the SE (ML-SE) and NW (ML-NW), respectively, with a PA≃ −41°. The ML-SE correspond to the loop seen by Aspin et al. (1993) in the NIR which primarily traces the emission of the vibrationally excited H₂.

Finally, we used the unsharp-masking technique to enhance the structural components of the MB in M 1–16, allowing us to identify new features (see the right image at the top panel of Fig. 1). Our images show a smaller lobe (SL-SE) inside the MB, which has a size of ~8″ and a PA≃ −39°. We did not observe any NW counterpart in our images. In this figure, we also found a central ellipse (CE) ≃ 9 × 7″ with its major axis at PA = −55°.

Figure 2 shows the image ratio [N II]/[O III] of M 1–16. The [O III] nebular emission is slightly stronger in the central part of M 1–16 (in white in the figure), faint in the ML (mostly ML-SE), and barely visible in the large lobes (see the middle image in the top panel of Fig. 1). The [N II] nebular emission, which is a known tracer of low-ionisation structures (Corradi et al. 1996), is concentrated in the lobes showing different structural components such as over-densities and knots (darker components in Fig. 2).

Fig. 1 Narrow band optical images of M 1–16 a seen in different emission lines. Top panel: M 1–16 as seen in Ηα+[Ν II] (left), [O III] (middle), and unsharp masked Ηα+[Ν II] (right) images obtained with MES in image mode. Bottom panel: Separated Hα (left) and [Ν II] (right) images obtained with the MEXMAN camera. Arbitrary contrast scales were set to all images to enhance the different features seen in the nebula. North is up and east is to the left in all of the images.

Fig. 2 Grey-scale image of the [Ν II]/[O III] of M 1–16 where dark values indicate a dominant [Ν II] emission line over [Ο III].

Fig. 3 M 1–16 as seen in Ηα+[Ν II] unsharp-masked images taken with MES in image mode. The MES and the B&Ch slits are overlaid in blue dashed lines and red boxes, respectively. The CS region is shown as a green box.

Fig. 4 Grey-scale image and contour plot of the [N II] emission line PV maps B, D, H, and I of the central region of M 1–16 using the spectral unsharp masking technique (see Appendix B). Both the image and the contour plot are in logarithmic representation with an arbitrary colourbar and levels, respectively.

3.2 Kinematics

Position-velocity (PV) diagrams corresponding to the slit paths in Fig. 3, were obtained and they are shown in Figs. A.1 and A.2, labelled from A to I. The position angle (PA) for slits A, B, and C correspond to −30°, −25°, and −40°, respectively. The other slits have PA = 60^° (D to G) and PA = 90^° (H and I). Spectra corresponding to slits passing through the central star (A, B, C, D, and H) were used to measure the systemic velocity of the shell, which was estimated as $V_{\text{sys}}^{LSR}=48\pm 1{\text{ km s}}^{-1}$. A visual inspection of PV maps A, B, and C (Figs. A.1 and A.2) shows that they are dominated by the central emission, which exhibits a widening at the position of the central star (also shown in PV maps D and H). The [N II] PV maps A, B, and C, also show an evident protrusion leaving the central region towards the SE, along the slit, and that is blue-shifted. This protrusion continues up to 20″ (or even more in PV map B), increasing its radial velocity to ≃120 km s⁻¹. This feature is not visible in Hα (just a marginal emission in B). PV map C also shows a fainter antisymmetrical protrusion towards the NW. The PV maps and measurements of the observed radial velocities are with respect to the aforementioned systemic velocity, $V_{\text{sys}}^{LSR}$.

On the other hand, PV map E shows a deformed ellipse with an expansion velocity of υ_exp ≃ 30 km s⁻¹ centred at ≃ −80 km s⁻¹ with respect to the systemic velocity. Additionally, PV maps F and G show double flat structures, spanning radial velocities of 140 and 250 ± 1 km s⁻¹, respectively, with the first one being centred at −80 and the last one at 120± 1 km s⁻¹. The PV map I marginally shows emission at the systemic velocity spanning 20″, as well as two faint knots at +50 and −50 km s⁻¹. PV maps at Hα A, B, C, D, and H are dominated by the central emission, also showing the He IIλ6560 Å emission line at the centre, whereas E, F, G, and I are a faint version of the corresponding PV maps in [N II]. Finally, a spectral logarithmic unsharp masking technique (see Appendix B) was applied on the [N II] PV maps A, D, H, and I showing that, at the centre, the spectral line is divided into some components in velocity, spanning radial velocities from −40 to +20 ± 1 km s⁻¹ (see contour-grey PV maps in Fig. 4).

3.3 Chemical and physical conditions

We identified and extracted five nebular zones along the major axis of the nebula. In addition, we extracted the area corresponding to the central star. The regions are identified from R1 to R5 and CS, and are shown in Fig. 3 as red boxes. We obtained the median of the spectra extracted from the inner regions, R1 and R2, to study the bright central part of M 1–16 and we added the multiple spectra taken for R3, R4, and R5 to study its faintest regions.

The calculations of the physical conditions, abundances, and respective uncertainties were performed with the code PYNEB (Luridiana et al. 2015), which is a tool dedicated to nebular analysis. The atomic data used to determine the ionic abundances are those described in Morisset et al. (2020).

Line fluxes for each extracted region were measured with the ‘splot’ task by fitting a Gaussian function to each line. The error budget included read-out noise and photon noise produced by both emission lines and the sky at the line position. However, uncertainties arising from wavelength calibration, flux calibration, or emission line fitting were not considered. To calculate the uncertainties on the line fluxes, we followed the independent and random errors approach (Taylor 1997). The error propagation (for all the parameters) was based on a Monte Carlo approach similarly to the work by Sabin et al. (2022) for the PN PC 22.

From the observed line intensities and their associated errors, PYNEB calculates the logarithmic extinction c(Hβ), the electronic density N_e, and temperature T_e in an iterative way, starting as a first step with a theoretical ratio of Hα/Hß = 2.85, T_e = 10⁴ K, and N_e = 10³ cm⁻³, assuming a case B recombination (Osterbrock & Ferland 2006). After convergence, the final values of c(Hβ), N_e, and T_e were obtained. The T_e was calculated using the [O III](λ5007 Å+λ4959 Å)/λ4963 Å and [N II](λ16548 Å+λ16583 Å)/λ15755 Å ratios, for high- and low-excitation regions, respectively. In the case of N_e, we used the [S II]λ6716 Å/λ6731 Å, [Ar IV]λ4711 Å/λ4740 Å, and [Cl III]λ5537 Å/λ5518 Å ratios, when available.

The first (top) half of Table 1 shows the intensity of the emission lines obtained for the different regions studied here. The second (bottom) half shows the extinction value, c(Hβ), and the physical conditions, electronic density (N_e), and temperature (T_e), obtained for the different regions. The emission lines were dereddened with the estimated value of c(Hβ) and using the extinction law of Cardelli et al. (1989). The values of N_e and T_e could only be estimated simultaneously for CS and regions R1, R2, R3, and R4, only. The interpretation of these results is shown in Sect. 4.2.

4 The nature of M 1–16

4.1 Morpho-kinematic structure

Morphology and kinematics of M 1–16 were modelled using SHAPEX, a software that specialises in morpho-kinematic reconstructions of gaseous nebulae (Steffen et al. 2011). Our analysis determined that the structures that provide the best simultaneous fit to both the imaging and spectra data are: (1) two large lobes, (2) one medium semi-lobe, (3) one (apparently compact) small semi-lobe, (4) one central ellipsoidal structure, (5) two distinct bipolar blobs, and (6) two knots in the northern lobe (see Fig. 5).

To determine the size, shape, and velocity law of each morphological structure, a visual fitting process was performed. The velocity law of a structure establishes the relationship between deprojected velocity and size or distance in a 3D vectorial basis. We used homologous velocity laws in which V ∝ r. Figures 6 and 7 display the comparison between these structures and the observations. The two blobs correspond exactly with the tips of the bipolar lobes, so we do not mention them afterwards. A summary of the main parameters of each structure, including their kinematic age and relationship with features labelled in Fig. 1, can be found in Table 2.

It is noteworthy that both knots were modelled using the same velocity law as the large lobes and placed close to the surface. As a result, the deprojected polar velocity of the large lobes and the deprojected velocity of the knots from the central star exhibit good agreement. Based on this model, it can be inferred that the three structures share the same kinematic age, suggesting that knots are part of or were impacted by the expansion of the bipolar outflow that created the large lobes. However, we also modelled the knots as small cylinders (1″ in diameter and length), with asymmetric density and an internal expansion velocity law (that is, with respect to the centre of each knot) of V_exp ≃ 100 km s⁻¹ × r. Our goal was to model bow-shock-like structures in order to reproduce the morpho-kinematical structures seen on the PV maps. The small size of the cylinders takes into account the spatially non-resolved bow shocks.

By considering the kinematic ages of all structures, it appears that the mass-loss history of this object may have progressed as follows: (1) the first collimated ejection of the SL stream (which may have been bipolar), (2) the second collimated ejection of the ML flow (which may have also been bipolar), (3) the bipolar collimated ejection of the LL, which interact with the LK and RK nodes, and finally (4) the ejection of the CE flow. The sequence of events leading to the observed structure may have been caused by a binary or multiple system, as similar scenarios have been suggested in the past (see Jones & Boffin 2017, for a review). The unusual morpho-kinematic structures of the SL and ML could result from the partial ionisation of bipolar structures that are not fully exposed to radiation, possibly due to the obstructing effect of a warped disk surrounding the central star. A careful analysis of the unsharp masking processed image (see the top-right panel of Fig. 1) reveals the presence of additional, faint lobes and ejections, some of which may be related to the semi-lobes. The rapid expansion of the knots, spreading at radial velocities of more than 200 km s⁻¹ with the tips reaching close to 250 km s⁻¹, could be explained by the effect of high-velocity ejected clumps colliding with the shell or the bipolar outflow interacting with previously ejected mass-loaded knots. This could also account for the slightly elevated temperature in [Ν II] relative to other regions.

High-velocity outflows are a common feature in PNe, especially in those with bipolar or point-symmetric shapes. These outflows are thought to be driven by collimated winds from the CSPN or its binary companion, and they interact with the slow-moving circumstellar envelope ejected during the AGB phase. One example of such a PN is M 3–38, which has a bipolar morphology and a high-velocity component defined as a jet (Rechy-García et al. 2022). This PN has been proposed to be the result of a CEE, where the primary star engulfed a main-sequence companion with a 17 day orbital period. Then, the CEE phase may produce a bipolar proto-PN with fast outflows and a circumbinary disk (García-Segura et al. 2021). A similar scenario has been suggested for other PNe with high-velocity outflows, such as MyCnl8, KjPn8, and also M 1–16 (Guerrero et al. 2020, and references therein). However, the origin and evolution of these outflows are still poorly understood and require further observations and modelling.

It is crucial to note that, from the comparison between observations and the model depicted in Fig. 7, it can be inferred that the nature of the core of M 1–16 is clumpy. While it is true that these clumps may generally adhere to the model’s pattern, their non-uniformity suggests a discontinuous ejection of material.

Fig. 5 SHAPEX reconstruction of M 1–16 based on all PV diagrams. The top panel shows the fill model and the bottom panel shows the inner-central part of the model.

Fig. 6 Comparison of the observed PV maps in the [O III] emission line for B, E, F, and G slits (top panel) and the SHAPEX synthetic spectra obtained from the adopted model (bottom panel).

Fig. 7 Similar to Fig. 6 but for A, D, and H slits of the central part of M 1–16.

4.2 Chemical and physical parameters

Table 1 shows the variation of the logarithmic extinction across the nebula. The most internal, and slightly more obscured regions’ CS, shows c(Hβ) = 0.88 with a mean value of ≃0.74 (taking into account the central region R2), whereas external regions R3, R4, and R5 have a similar mean value of ≃0.70.

The electronic density calculated with the low-ionisation species, [S II], shows fluctuations along the central region CS, RI, R2, and R3 regions, with N_e being a factor ~4 lower in R1, R2, and R3 compared to the CS. This implies that the collisional excitation rate is less important in R1, R2, and R3 with respect to the CS. The estimation of N_e with the high-ionisation species, that is [Cl III] and [Ar IV], could only be performed in the internal region CS.

According to Fig. 3, we are currently measuring the N_e from the SL-SE and the MB in R2 and the CS, respectively. The value of N_e([S II]) = 2300 cm⁻³ of the CS is in agreement with the value obtained by Corradi & Schwarz (1993) as measured in the central region of M 1–16, and much higher than the one estimated by Schwarz (1992). The N_e = 610 cm⁻³ for the LL-NW, as measured in R3, however, is higher than the value obtained by Corradi & Schwarz (1993) of 200 cm⁻³ (for that low regime; see Sect. 3.2). The T_e, estimated from [N II] and [O III] for low- and high-ionisation regions, respectively, does not show fluctuation and a mean value of ~14 kK is obtained.

Estimated ionic and elemental abundances are given in Tables 3 and 4, respectively. To estimate the elemental abundances, we used the ionisation correction factors (ICF) of Delgado-Inglada et al. (2014). When compared to solar abundances derived by Lodders (2010; see Table 4), M 1–16 is clearly N and Cl overabundant with (N/H) ≃ 3.53 (N/H)_⊙ and (Cl/H) ≃ 4.44 (Cl/H)_⊙, respectively. On the contrary, M 1–16 shows a sub-solar abundance in O, S, and Ar with (O/H) - 0.60 (O/H)_⊙, (S/H) ≃ 0.098 (S/H)_⊙, and (Ar/H) ≃ 0.45 (Ar/H)_⊙, respectively. In the case of He, the variation is only 1.82 with respect to the solar abundance.

The elemental abundance results indicate that M 1–16 is a Type-I PN because of the log(N/O) = −0.103 and He/H = 0.153 obtained for the CS (defined as log(N/O) > −0.3 and He/H≥0.125; Stanghellini & Haywood 2010). We also estimated the CSPN mass by using method three described in Maciel et al. (2010), yielding ~0.74 M_⊙. According to Miller Bertolami (2016), such a CSPN mass indicates that M 1–16 evolved from a progenitor star with a mass between 3 and 4 M_⊙ for Z = 0.02 or 2.5 M_⊙ for Z = 0.001.

Schwarz (1992) and Aspin et al. (1993) argued that M 1–16 is a young object, the so-called proto-PN, and the latter estimated a T_eff = 35 000 K by means of the ratio between He I 2.058 μm and Brγ. However, the presence of high excitation ions, such as He IIλ4686 Å, [O ΙII]λ5007 Å, and [Ar V]λ6435Å, in the optical spectra, indicate a T_eff of at least 50 000–60 000 K. We used the 3MdB (Morisset et al. 2015), a database of photionisation models made with the CLOUDY code (v17.01; Ferland et al. 2017), to estimate the T_eff of the CSPN of M 1–16. We looked for models made with the theoretical atmosphere models from Rauch (2003) and the ones that simultaneously reproduce the log(He ΙIλ4686 Å/H/β), log([O III] λ5007 Å/[N II] λ6583 Å), log([Ar V] λ7006 Å/[Ar IV] λ4740 Å), and log([O III] λ5007 Å/[O I] λ6300 Å). A total of 13 models were identified to fit the observed ionisation stage of M 1–16, corresponding to a T_eff ≃140 000 K and a log(L/L_⊙) ≃ 2.3. Therefore, according to the theoretically evolutionary post-AGB tracks of Miller Bertolami (2016), this sets a progenitor mass between 2.0–3.0 M_⊙ and a mass for the CSPN between 0.618–0.713 M_⊙ depending on the metallicity (comparable with the estimated value for the CSPN mass using the log(N/O)).

5 Conclusions

From the kinematic study, we derived a systemic velocity for M 1–16 of $V_{\text{sys}}^{LSR}=48\pm 1{\text{ km s}}^{-1}$ and a deprojected expansion velocity of ≃30 km s⁻¹ for the brightest lobe, ML-SE, as obtained from the PV maps presented here. A detailed morpho-kinematic model was build using the computational model SHAPE fitting our data and improving our interpretation. The model is composed by two pairs of lobes, corresponding to SL and ML, and one medium semi-lobe, one central ellipsoidal structure, and two pairs of knots located to the N at the tip of the LL-NW. According to the model, M 1–16 first ejected the SL stream with a deprojected velocity of 125 km s⁻¹ (i = 10°) followed by a second collimated ejection of the ML flow at 140 km s⁻¹ (i = 40°). Then, the bipolar collimated ejection of the LL occurred and interacted with the LK and RK nodes. Finally, the ejection of CE flow occurred. A binary system or multiple system at the centre of M 1–16 is proposed to explain the sequence of mass-loss events in this nebula, as similar scenarios have been observed in other cases (see Jones & Boffin 2017, and references therein). The unusual morpho-kinematic structures of the SL and ML could result from the partial ionisation of bipolar structures that are not fully exposed to radiation, possibly due to the obstructing effect of a warped disk surrounding the central star.

The chemical analysis along with the kinematical age obtained from the morpho-kinematic model demonstrated that M 1–16 may be an evolved PN instead of a proto-PN as previously suggested. The presence of high ionisation lines, such as [O II], [Cl III], and [He II], and the chemical abundances indicate that the CSPN have a mass of ~0.74 M_⊙ and evolved from a progenitor star with a mass between 2.5 and 4.0 M_⊙ (depending on the metallicity); a T_eff ≃ 140 000 K and a log(L/L_⊙) ≃ 2.3 were estimated by using the 3MdB photoionisation models in conjunction with theoretical post-AGB evolutionary models indicating a CSPN mass in the range of 0.618–0.713 M_⊙ and a progenitor mass between 2.0–3.0 M_⊙.

Acknowledgements

M.A.G.M. acknowledges support from the ACIISI, Gobierno de Canarias and the European Regional Development Fund (ERDF) under grant with reference PROID2020010051 as well as from the State Research Agency (AEI) of the Spanish Ministry of Science and Innovation (MICINN) under grant PID2020-115758GB-I00. M.A.G.M. also acknowledges the postdoctoral fellowship granted by the Instituto de Astronomía of the Universidad Nacional Autónoma de México. This work was supported by UNAM-PAPIIT IN106720 grant. L.S. acknowledges support from grant UNAM-PAPIIT IN110122. S.Z. acknowledges support from the ITE-UNAM agreement 1500-479-3-V-04, and from O. Velázquez-DISC/ITE. This work is based upon observations carried out at the Observatorio Astronómico Nacional on the Sierra San Pedro Mártir (OAN-SPM), Baja California, Mexico. We thank the daytime and night support staff at the OAN-SPM for facilitating and helping obtain our observations, in particular to the telescope operators Mr. Gustavo Melgoza (also known as ‘Tiky’) and Mr. Felipe Montalvo. This research made use of APLpy, an open-source plotting package for Python (Robitaille & Bressert 2012). We thank the referee for they comments that improved the manuscript.

Appendix A Position-velocity maps

Fig. A.1 PV maps of the [Ν II] emission line for the selected slits labelled in the bottom right corner of each panel.

Fig. A.2 PV maps of the Hα emission line for the selected slits labelled in the bottom right corner of each panel.

Position-velocity maps were extracted from all the spectra. Figure A.1 corresponds to emission in [Ν II], whereas Fig. A.2 is for Hα. Some notable features are apparent in the panels corresponding to the slits crossing the central star (A, B, C, D, and H): (a) In Fig. A.1, some artefacts can be observed, such as a ‘big bubble’ surrounding all the emission, as well as two horizontal lines above and below the continuum of the central star. These known artefacts were briefly mentioned by Meaburn et al. (2003) and examples can be found in the database of López et al. (2012). (b) These artefacts are also present in the corresponding panels in Fig. A.2, where the presence of the He IIλ6560Å emission line is also noticeable. Moreover, it is important to note that the emissions observed in panels F and G from both Figs. A.1 and A.2 are authentic, and they should not be confused with the artefacts previously mentioned, as they correspond to gaseous knots with low brightness.

Appendix B The spectral unsharp masking technique

Unsharp masking is a technique used in image processing to enhance the perception of edge details in an image. It works by subtracting a blurred version of the image from the original image, accentuating the edges and fine details present in the original but not in the blurred version. The result is a sharper, more detailed image. The technique is often used in digital photography and image editing software to improve the clarity and sharpness of photographs. In signal processing, an unsharp mask is equivalent to a filter that amplifies the high-frequency components of a signal.

The spectral unsharp masking (SUM) technique is presented here as a variant. In astronomical image processing, a blurred image is typically obtained by applying a Gaussian convolution. The main parameter of this convolution is the standard deviation or σ value, which is related to the ‘seeing’ of the observations. A 1-σ convolved image typically has no effect, as the seeing already accounts for any image blur. Commonly, a 5-σ convolution or higher is used to obtain the blurred image. We shall call this variable σ_spatial.

In the case of spectra, one axis corresponds to the slit (spatial axis), and the other corresponds to the spectral axis (usually wavelength, but in our case, velocity). The spectral resolution is measured by fitting Gaussian curves to comparison arc lines. This measurement considers all types of broadenings, including those caused by local effects and instrumental errors, with natural broadening being dominant. Natural broadening also has a Gaussian profile, and we usually report the FWHM as the main parameter to describe it. FWHM is related to σ as $\text{FWHM}=2\sqrt{2\text{ ln }2}\sigma$. We shall call this variable σ_spectral.

Since the values of σ_spatial and σ_spectral are not related to each other, and are typically different in terms of pixels, we cannot use a two-dimensional circular Gaussian, as is done in imaging. Instead, we must use a two-dimensional elliptical Gaussian. Therefore, if we want to apply a 5-σ blur to a PV map, we must multiply the number five by the corresponding value (in pixels) to both σ_spatial and σ_spatial separately. These values are the major and minor axes of the Gaussian ellipse used to obtain the blurred image when applied to the corresponding axes. The other steps for getting the unsharp masking are the same as for imaging.

References

All Tables

All Figures

Fig. 1 Narrow band optical images of M 1–16 a seen in different emission lines. Top panel: M 1–16 as seen in Ηα+[Ν II] (left), [O III] (middle), and unsharp masked Ηα+[Ν II] (right) images obtained with MES in image mode. Bottom panel: Separated Hα (left) and [Ν II] (right) images obtained with the MEXMAN camera. Arbitrary contrast scales were set to all images to enhance the different features seen in the nebula. North is up and east is to the left in all of the images.

In the text

	Fig. 2 Grey-scale image of the [Ν II]/[O III] of M 1–16 where dark values indicate a dominant [Ν II] emission line over [Ο III].
In the text

	Fig. 3 M 1–16 as seen in Ηα+[Ν II] unsharp-masked images taken with MES in image mode. The MES and the B&Ch slits are overlaid in blue dashed lines and red boxes, respectively. The CS region is shown as a green box.
In the text

	Fig. 4 Grey-scale image and contour plot of the [N II] emission line PV maps B, D, H, and I of the central region of M 1–16 using the spectral unsharp masking technique (see Appendix B). Both the image and the contour plot are in logarithmic representation with an arbitrary colourbar and levels, respectively.
In the text

	Fig. 5 SHAPEX reconstruction of M 1–16 based on all PV diagrams. The top panel shows the fill model and the bottom panel shows the inner-central part of the model.
In the text

	Fig. 6 Comparison of the observed PV maps in the [O III] emission line for B, E, F, and G slits (top panel) and the SHAPEX synthetic spectra obtained from the adopted model (bottom panel).
In the text

	Fig. 7 Similar to Fig. 6 but for A, D, and H slits of the central part of M 1–16.
In the text

	Fig. A.1 PV maps of the [Ν II] emission line for the selected slits labelled in the bottom right corner of each panel.
In the text

	Fig. A.2 PV maps of the Hα emission line for the selected slits labelled in the bottom right corner of each panel.
In the text