© The Authors 2026

1 Introduction

The transition from asymptotic giant branch (AGB) stars surrounded by their circumstellar envelopes (CSEs) to a white dwarf (WD) and their planetary nebulae (PNe) is a brief and spectacular metamorphosis experienced by low-to-intermediate mass (∼0.8−8 M_⊙) stars, and it is still poorly understood, particularly regarding the physical mechanisms driving the transformation (see Balick & Frank 2002). Far from simply being round fossils of AGB quasi-spherical stars and the isotropically expanding envelopes of AGB stars (AGB CSEs), ∼80% of PNe depart from spherical symmetry and display a wide range of morphologies organised around (at least) one axis, whose origin we have been unable to grasp fully despite thorough efforts in the past decades (e.g. Balick & Frank 2002; Jones & Boffin 2017, Stanghellini et al. 2016). Similarly, the shaping of nebulae around intermediate post-AGB stars, not yet able to ionise the surrounding nebula, termed pre-planetary nebulae (pPNe), have proved to be challenging as well (e.g. Sahai & Trauger 1998; Ueta et al. 2000; Sahai et al. 2007; Lagadec et al. 2011; Sahai et al. 2011; Szczerba et al. 2012). It has become increasingly clear that a source of angular momentum provided from a binary or sub-stellar companion is a key ingredient to this puzzle, which almost certainly involves the launching of jets likely powered by accretion discs in rotation. However, the specifics of how this happens remain stubbornly elusive, with hydrody-namic/magnetohydrodynamic (HD/MHD) models struggling to match the vast zoo of PNe/pPNe shapes from mostly uncertain initial conditions (e.g. Mastrodemos & Morris 1999; Icke 2003; García-Segura et al. 2005; Huarte-Espinosa et al. 2012; Soker 2015).

In this context, morpho-kinematic modelling of PNe/pPNe has proven to be a valuable tool for investigating their 3D spatial distribution and velocity field, which includes their orientation with respect to the plane of the sky and their kinematic (distancedependent) age (e.g. Solf & Ulrich 1985; Santander-Garcia et al. 2004; Akashi & Soker 2013; Akashi & Soker 2018). In this respect, the SHAPE software (Steffen et al. 2011) has become a standard in the PNe community, as it allows relatively easy modelling of these nebulae with custom geometries and a few parameters. Over the years, it has been used to successfully describe the present-day ejecta of a number of PNe displaying the most intricate morphology and kinematics (e.g. García-Díaz et al. 2012; Santander-García et al. 2012; Clyne et al. 2015).

Following the release of SHAPE, our group developed shapemol (Santander-García et al. 2015), a plug-in that provides precise calculations of excitation and radiative transfer of rotational lines of ¹²CO and ¹³COusing the Large Velocity Gradient (LVG) approximation, as established by Castor (1970). This tool generates synthetic spectral profiles and maps to match observational data. Simultaneously analysing low- and high-excitation lines allows for the determination of gas excitation conditions, more specifically, the density, temperature, and molecular abundance of each species relative to the total number of gas particles. From these data we can estimate the mass distribution in AGB envelopes and PNe, offering insights into molecular chemistry conditions and contributing to the refinement of formation models (e.g. Santander-García et al. 2012; Santander-García et al. 2017; Doan et al. 2017; Díaz-Luis et al. 2019).

As millimetre-range radio telescopes and interferometers improve in bandwidth and sensitivity, a plethora of spectral lines from species other than ¹²CO and ¹³CO are detected in observations of pPNe and PNe, allowing for a more complete description of the physico-chemical properties of the studied objects. Tracers such as HCN or HCO⁺, for instance, can reveal dense and hot pockets or shocked regions within the molecular envelopes of pPNe and PNe, while analysis of the spatial distribution of the HNC/HCN abundance ratio can be used as a diagnostic of ultraviolet radiation (Viti et al. 2002; Jørgensen et al. 2004; Bublitz et al. 2019; James et al. 2020).

Bearing in mind these new possibilities to gain insight into pPNe and PNe formation, we have expanded shapemol to include ten new molecular species and implemented a code to improve the direct matching of interferometric observations. Using these enhancements, we present new findings from the concurrent analysis of numerous IRAM-30 m and HER-SCHEL/HiFi spectra along with NOEMA interferometric maps focusing on different species within the pPN M1-92.

The pPN M1-92, also known as Minkowski’s Footprint (Herbig 1975), is one of the most iconic objects of its class. First discovered by Minkowski ca. 1946 (Minkowski 1946), it consists of a bipolar reflection nebula around a B-type post-AGB central star, about 12″×6″ in size, with the major length following a conspicuous symmetry axis oriented at a PA ∼ 310° (Herbig 1975; Solf 1994; Bujarrabal et al. 1998b; Ueta et al. 2007). In the optical, the reflection nebula depicts a bi-lobed structure divided by a thick equatorial waist that partially obscures the south-east lobe. This orientation agrees with the expansion velocities measured in emission lines, which show approaching and receding movements in the north-west and south-east lobes, respectively (Herbig 1975; Solf 1994). These line emissions arise from a pair of compact knots placed along the symmetry axis, about 2″.5-3″.0 away from the centre, roughly in the middle of each lobe (Herbig 1975; Solf 1994; Bujarrabal et al. 1998b). They are detected in optical forbidden lines and vibrationally excited H₂ (see Davis et al. 2005, and references quoted before), indicating the presence of active shocks.

The nebular composition is dominated by the molecular gas, which has been studied through observations of low-J transitions of ¹²CO and ¹³CO, including interferometric maps at 075-175 spatial resolutions (Bujarrabal et al. 1997, 1998a; Alcolea et al. 2007). The molecular gas is located in the equatorial structure dividing the two lobes, in the walls of these two formations, and at the two polar ends of the nebula (the polar tips). In contrast to the optical images, because of the negligible dust opacity at millimetre-wavelengths, the emission from molecular gas also displays a remarkable equatorial mirror symmetry, i.e. both lobes are very similar to each other. A linear velocity gradient of 7.4 km s⁻¹ per arc-sec. was estimated, with expansion deprojected velocities up to 70km s⁻¹, which for a distance of 2.5 kpc translates into a kinematic age of 1200 a. This velocity gradient seems to hold in all directions, suggesting that the nebula was accelerated to its present velocities in a single brief (<100 a) event (Alcolea et al. 2007). In addition to CO and H₂, the nebula has also been detected in other molecular species such as OH, SO, SO₂, and SiO, clearly indicative of O-rich composition, but also in CS, HCN, and HCO⁺, providing evidence of a rich and diverse chemistry (Seaquist et al. 1991; Alcolea et al. 2022).

As stated before, in this paper, we fit the observational results of the four most abundant isotopologues of CO, HCO⁺, H¹³CO⁺, HCN, and H¹³CN from M1-92 while putting into practice the new capabilities of the upgraded shapemol code. Details on these observations are given in Sect. 2. In Sect. 3, we describe the improvements in both SHAPE and shapemol developed for this work. Our new model for M1-92 is presented in Sect. 4. Implications of this model on physical and chemical characteristics are discussed in Sect. 5. Finally, the main results of this paper are summarised in Sect. 6.

A very preliminary work was published in Masa et al. (2024), where a first approximation of the model and only some of the observational data were presented. Here we provide the final version of the model in addition to the complete work, from the new capabilities of SHAPE and shapemol to full analysis of the nebula’s properties, and include all the observational data used.

2 Observational data

The data presented in this paper were obtained with the IRAM 30 m-MRT, the IRAM NOEMA interferometer, and the HiFi instrument onboard the Herschel Space Observatory (30 m, NOEMA, and HiFi hereafter).

2.1 Single-dish data

The 30 m data were acquired through various projects, namely 050-15, 160-15, and 047-16. They comprise the J =1−0 and J=2−1 lines of the four most abundant CO isotopologues, ¹²CO, ¹³CO, C¹⁷O, and C¹⁸O, and the J =1−0, J=2−1, and J=3−2 lines of HCN, H¹³CN, HCO⁺, and H¹³CO⁺. Data were obtained with the Eight MIxer Receiver (EMIR) receiver in the 3, 2, and 1.3 mm bands, using as a backend the fast Fourier transform units (FFTs) in its wide configuration, which provides a native spectral resolution of 200 kHz. This results in a spectral resolution of 0.7-0.22 km s⁻¹, although final spectra were resampled to a uniform resolution of 3.25 km s⁻¹ for increasing their S/N ratio, also facilitating the comparison with previous NOEMA observations and model results. The observations were performed in the wobbler switching mode (WSW), alternating between the ON and OFF positions at a frequency of 0.5 Hz. The OFF positions are symmetrically taken 100″ away in the azimuth axis to ensure a proper subtraction of the sky and instrumental contributions. This method provides very flat and stable baselines. Observations are automatically calibrated in T_* units, which includes every 20 mins estimates of the receiver noise and also corrections for the perturbations introduced by the atmosphere. Atmospheric absorption was corrected from the sky emissivity obtained in these calibration measurements and a model for the atmosphere at the 30 m site, near Pico Veleta (Sierra Nevada, Spain), at an altitude of 2850 m. Pointing was checked by the observation of quasars and other radio-continuum sources close in the sky to our target, about every two hours. After applying these corrections, we can assure that the pointing of the telescope is better than 2″. Relative calibration between the different runs were performed by the observation of strong spectral line emitters: specifically, the pPNe CRL 2688 and the young PN NGC 7027. For CRL 2688 we adopted integrated area values of 218, 170, 56.5, 151, 324, and 300 K km s⁻¹ (in T^*_a scale) for ¹²CO J =1−0, ¹²CO and ¹³CO J=2−1 and HCN J =1−0, J=2−1, and J =3−2, respectively. For NGC 7027, we adopted integrated area values of 246, 238, 13.6, 9.8, 27.6, 21.5, 26.6, 49.6, and 42.0 K km s⁻¹(in T^*_a scale) for ¹²CO J =1−0, ¹²CO and ¹³CO J=2−1 HCN J =1−0, J=2−1, and J=3−2, and HCO⁺ J =1−0, J=2−1, and J =3−2, respectively. After this final recalibration, we adopt an absolute flux uncertainty of ±5-10% in the 3 and 2 mm bands, and of ±10-15% at 1.3 mm. After eliminating bad and noisy scans, data from both horizontally and vertically polarised receivers were averaged. Zero- or first-degree polynomial baselines were also removed. Finally, all data have been rescaled into T_mb (K) using the main-beam efficiencies given in the EMIR user manual: 1.20 and 1.55 for CO J =1−0 and J=2−1 respectively, and 1.15, 1.40, and 1.70 for HCO⁺ and HCN J =1−0, J=2−1, and J=3−2 respectively.

The HiFi/HSO data for the J=5−4 and J=9−8 lines of ¹²CO have been taken from the programme GT1_dteyssie_1 (PI: D. Teyssier), which have already been presented in the work of Lorenzo et al. (2021). Those for the J=7−6 one have been taken from the project OT1_vbujarra_4 (PI: V. Bujarrabal). In this latter case, data calibration and reduction have been performed as in the case of the HIFIStars programme (PI: V. Bujarrabal), fully described in the work by Alcolea et al. (2013).

2.2 Interferometric data

NOEMA interferometric data, comprising the maps of the J=2−1 line of ¹³CO, C¹⁷O, C¹⁸O, HCN, H¹³CN, HCO⁺, and H¹³CO⁺, were obtained in project W17BJ. Data for CO isotopologues were acquired in the 1.3 mm band (NOEMA band 3) while those for HCN, HCO⁺, and their ¹³C substitutions were acquired using the 2 mm receiver (NOEMA band 2). The observations were performed in the most extended A-configuration, together with the more compact C-configuration, resulting in final spatial resolutions of about 0″.7. The PolyFix correlator was configured to provide an 8+8 GHz full bandwidth coverage at a resolution of 200 kHz. Final maps were also resampled to a spectral resolution of 3.25 km s⁻¹. Observations happened with 9 antennas in March and April 2018 for two spectral setups, in A and C configurations: at the LO frequency 231.022 GHz (setup 1) five tracks were obtained, including one in C configuration, and at the LO frequency 168.022 GHz (setup 2) four tracks were observed, one in C configuration. In all the tracks a bright quasar was observed to perform the RF bandpass calibration. Two close calibrators (1901+319 and 1932+204) were observed every 25 min to determine instrumental phase and amplitude gain changes in time, and also to ensure a proper pointing and focus during the observations (every 25 min or 50 min). The main absolute flux calibrator in all the tracks was MWC349, with an adopted flux of 1.86 Jy at 219.6 GHz and 1.64 Jy at 178.8 GHz (with S.I. 0.6). The flux calibration accuracy was found to be better than 10% in both bands, considerably better at 2mm. Data calibration was made in CLIC by using the standard calibration procedures.

Whole band datasets were first mapped and then cleaned and inspected for line absorption/emission. After flagging channels affected by spectral signatures, data were channel-averaged in the uv-plane and imaged and cleaned once more for obtaining maps of the continuum emission at the different observed frequencies (156 204.7-164 310.9, 171 690.6-179 797.3, 219 188.8227 303.8, and 234 930.9-242 795.8 MHz). Pure (line-free) continuum uv-data were subtracted from the original datasets for obtaining pure (continuum subtracted) spectral line emission uv-data. Final datasets represented in LSRK-velocity (assumed systemic velocity of −0.5 km s⁻¹ in the LSRK frame) channels were obtained after correction for the corresponding rest frequencies (obtained from the JPL line catalogue) and resampling at a spectral resolution of 3.25 km s⁻¹ . This resample increases the signal-to-noise ratio and, given the velocity gradient in this source, translates into a spatial resolution along the line of sight direction of about 0″5, which is comparable to the spatial resolution of the interferometric maps. In this paper, we show channel maps for a bandwidth of 140 km s⁻¹ that includes the whole detected emission. Image synthesis was performed using MAPPING within the GILDAS software¹. For the dirty mapping, we adopted a uv-plane sampling of 7.5 m (half the diameter of the antennas), a pixel size of 0″.1×0″.1, appropriate for the final attained spatial resolution, and a natural weighting scheme. Cleaning was performed using the Hogbom method, a gain of 0.15, and support polygons to restrict the area looking for clean components. We checked the convergence of the cleaning procedure by inspecting the cumulative cleaned flux. Finally, we also compared the total fluxes from the maps with those obtained with the 30 m, concluding that there is no flux filtered out in the interferometric data as both measurements are compatible within the expected absolute calibration uncertainties. In this paper, we only present the ¹³CO, C¹⁷O, C¹⁸O, HCO⁺ and HCN J=2−1 maps, with a resolution of around 0″.50×0″.75, and a position angle (PA) of the beam of 76° for CO isotopologues and 86° for the other two species. These maps can be found in Figs. B.1-B.5, together with their moment zero maps and integrated emission spectra in Fig. B.6 and a tabulation of integrated fluxes from both IRAM 30 m and NOEMA observations in Table A.2.

3 Software

3.1 shapemol update

SHAPE is a software code for the numerical simulation of the emission from gas nebulae (Steffen et al. 2011). The original version of SHAPE is not intended for the modelling of the emission of molecular lines; this is done by the accompanying code shapemol (see Santander-Garcia et al. 2015). In the version of shapemol presented there, only the treatment of ¹²CO and ¹³CO was implemented. We introduce a new version of shapemol to work along with the newest version of SHAPE², which includes more molecular species and an updated treatment of the two most abundant CO isotopologues. The molecules chosen for this update are species relevant in the radio astronomy field but also with simple structure (linear molecules), where the vibrational modes can be ignored when calculating the radiative transfer at the studied frequencies (in the millimetre and sub-millimetre domains).

As stated previously, shapemol is a complementary code for SHAPE that computes synthetic line profiles and maps for the molecular line emission of a numerical nebula model. In summary, shapemol solves the statistical equilibrium population of a given molecular species using the large velocity gradient (LVG) approximation formalism (see Castor 1970), transforming the problem into a local one: the excitation properties of any molecular species at a given point does not depend on the properties of the rest of the nebula as a result of the different Doppler shifts. This calculation requires the simultaneous computation of several level populations that depend on both the collisional and radiative excitation transitions. Under the LVG formalism, since this is transformed into a local problem, the only relevant parameters are the density n, temperature T, relative abundance X of the considered species, and the logarithmic velocity gradient ε = dlog(V)/dlog(r) = (r/V)(dV/dr), all referred to the particular point of the nebula for which the computations are being performed. The coupling of any given point with the rest of the nebula is controlled by the ratio of the macroscopic velocity to the local velocity dispersion: the larger the value, the more decoupled the different nebular points are. As a consequence, e determines the effective escape probability of any given photon and so the importance of radiative trapping. LVG approximation applies to cases of an increasing radial velocity, which means an ε ≥ 0. This version of shapemol has only been tested in these cases, and therefore, we recommend using it only for monotonously increasing radial velocities.

Although the LVG formalism was developed for large velocity gradients, i.e: large values of e, due to its simplicity, it is also applied for low e values where, even though the expected accuracy in these cases is lower, its low computational cost and still very reasonable results, make it a very useful tool for a modelling process (see Bujarrabal 1994a,b; Bujarrabal & Alcolea 2013).

shapemol uses values of the absorption and emission coefficients, k_v and j_v respectively, of a set of rotational lines of a given molecular species. These values have been pre-computed for an ample range of densities, temperatures, relative abundances and values of e, and stored in a structured set of tables where shapemol seeks the most appropriate values for a given nebular point (as a function of its local n, T, X, and e values), and transfers them to the main code of SHAPE for the production of synthetic spectral profiles and maps to be compared to real observations. Similar tables as the original ones are now available for C¹⁷O, C¹⁸O, HCN, HCO⁺, HNC, SiO, CS, N₂H⁺, H¹³CN, and H¹³CO⁺. The tables for ¹²CO and ¹³CO introduced in the original work (Santander-García et al. 2015) have also been extended to match the new ranges of temperature, density and abundance covered by the new species (see Table A.1).

For the pre-computation of the molecular excitation (the k_v and j_v tables), we only considered collisional rates with H₂. Other collisional partners, such as He and electrons, as well as more complex molecular structures, such as the hyperfine structure, are not taken into account. Neglecting the hyperfine structure in shapemol may have some impact in the J =1−0 lines of HCN and N₂H⁺, resulting in slightly more opaque lines and narrower spectra. We note that in some sources, such as the one presented in this paper, this effect is not relevant due to the line width being dominated by the local turbulence and the source’s kinematics. The collisional coefficient data are from the Lambda Database³ when they reached the temperature range we aimed for and from other more specific works otherwise. In the cases where separate tables for collisions with ortho-H₂ and para-H₂ are available, an ortho/para abundance ratio of three has been assumed. When not stated otherwise, Einstein coefficients and level multiplicities have also been taken from the Lambda Database, while the molecular constants used were provided by the Cologne Database for Molecular Spectroscopy (CDMS)⁴ and the National Institute of Standards and Technology (NIST)⁵ databases. (For further details on shapemol’s characteristics, see Santander-García et al. 2015.)

For calculations of these pre-computed tables, the same nonlocal thermodynamic equilibrium (LTE) model based on the LVG approximation has been followed, and the general format has been maintained. Therefore, the absorption and emission coefficients have been calculated for ground-state rotational transitions from J =1−0 up to J =17-16, for each combination of temperature, density and relative abundance, in the cases of e being 0.2, 1 and 3. Temperature ranges from 5 to 1000 K in steps of 5 K and density from 10² to 10¹² cm⁻³ in multiplicative factors of $\sqrt[4]{10}$ for all CO isotopologues, $\sqrt[8]{10}$ for SiO and CS, and $\sqrt[16]{10}$ for HCN, HNC, H¹³CN, HCO⁺, H¹³CO⁺, and N₂H⁺.

The relative abundance range has been chosen for each molecule so that the density column calculated for a radius of 10¹⁷cm covers those typically found in different objects, from AGB shells to planetary nebulae. The step chosen also has a multiplicative factor, this time of $\sqrt[NT-1]{5.3{10}^4}$ for CO isotopologues, SiO, CS, HCO⁺, and H¹³CO⁺, and $\sqrt[NT-1]{{10}^7}$ for HCN, HNC, H¹³CN, and N₂H⁺, where NT is the number of different abundance values and so the number of tables provided. A total of 52 tables for each species are available, except for HCN, H¹³CN, HNC, and N₂H⁺, whose typical range of values is about 2-3 orders of magnitude larger than for other species, and in which case 76 tables have been calculated, and SiO, calculated for a slightly higher abundance range with a total of 56 tables. This difference in number of tables aims to keep a similar abundance gap between tables. A summary of the characteristics and parameter coverage is shown in Table A.1.

Contrary to abundance, where different ranges require different iteration steps, the differences between multiplicative factors in density for each species are purely due to converging purposes. When it was needed, smaller density steps were taken to ensure smoother transitions between iterations and guarantee the full computation of each table. This is complemented by another strategy implemented in the temperature iterations, where the number of energy levels considered for calculation is reduced depending on the temperature considered, as lower temperatures will only significantly populate the lower levels, generating calculation noise from upper-level populations and increasing the computational cost.

The accuracy of these calculations has been checked for the most abundant isotopologues, assuming an analogue behaviour for their less common counterparts. These checks where done by comparing the opacities obtained during the computation of the k_v and j_v tables with those provided by RADEX⁶ (van der Tak et al. 2007), as well as through the result of radiative transfer expected in the simple analytic case of local thermodynamic equilibrium.

Differences with RADEX, which provides non-LTE molecular radiative transfer comparable to the LVG approximation in the case of a linear velocity gradient, are in the 5-8% range when the same collisional coefficients are used, which is the case of CO isotopologues, and between 5 and 70% otherwise. These discrepancies can be explained by the collisional coefficients data themselves, as differences of up to a factor 2 can be found between the different datasets used by RADEX and this work (for analysis of both datasets for CS see Denis-Alpizar et al. 2018). These comparisons were made only when collisional data were available for both RADEX and shapemol. Temperature extrapolations or hyperfine structure were not considered when checking these results, as the huge discrepancies in those cases make the comparison irrelevant. Details on the differences between RADEX and shapemol, as well as collisional coefficients data sources, can be found in Table 1.

However, the most significant test of the final SHAPE+shapemol results is the comparison of its calculations in a case where we can analytically compute the solution. The test was run in a homogeneous thermalised sphere in which the kinetic temperature matches the excitation temperature of the molecule in those conditions while expanding with a linear velocity law. With Planck’s law and the LVG approximation, the final intensity emitted at any point on the sphere can be obtained. In this case, the error is between 0.05 and 10%, generally decreasing as the opacity increases. This is the expected result, since as the opacity increases the lines tend to thermalise, converging to the kinetic temperature regardless of the approximation. Details on the errors for each molecule and opacity under the conditions required for the test are shown in Fig. 1.

These new tables that expand and update shapemol’s capabilities can be downloaded via SHAPE’s website⁷ or the available electronic material of this work⁸. For checking purposes, tables of ¹²CO and ¹³CO applied to LTE conditions are also available for the same ranges of temperature, density and relative abundance.

3.2 SHAPE

Over the course of the expansion of shapemol, the SHAPE software has undergone a few upgrades in the frame of this work. Apart from fixing small errors and bugs, version BETA.2.0.0 has two new features that are useful for the simultaneous modelling of emissions from different transitions or molecular species. We briefly describe them below.

A new treatment of the cosmic microwave background (CMB) has been implemented in the Spectrum tab of the Render module. As in previous versions, by default, a black-body with a temperature of 2.73 K is assumed for the background, its emission being taken into account at each location in the first step of the radiative transfer calculation, and subtracted at the last step. In the new version, however, a background with a custom, user-defined temperature can be used instead, and the final subtraction can be disabled. This allows the modelling of sources against extended backgrounds with temperatures higher than the CMB.

SHAPE now also includes a new scripting module that enables automation of the workflow. Employing its own simple scripting language, for which some examples are provided, the user can sequentially render a nebula in different transitions and save the resulting spectral profiles or map data cubes, for example, or explore the parameter space by automatically rendering a grid of models with different values of the desired parameters, and saving the results for later inspection. This new feature does not require shapemol to work.

Fig. 1 Evolution of the error as opacity decreases, calculated as the discrepancy between shapemol and the theoretical calculations with respect to the theoretical value. Only the opacity conditions where the molecule had thermalised were used.

3.3 Model-data comparison

With every rendering, SHAPE computes the emission along the line of sight through a ray-tracing algorithm. This produces a synthetic map with the same spatial resolution as the 3D model. This map then gets convolved with the beam provided for each transition and observing mode, creating the final dataset that is compared to the observational data. In the case of single-dish, the data is already ready for comparison, while for interferometric maps we added an extra step that subjects the raw data from SHAPE to an ‘interferometrisation’ process with the same conditions as the observation. For this process to be reliable, we must first SHAPE-model an area wide enough so there is no significant signal at the borders. This image should then be enlarged by a factor of two to avoid aliasing problems in the next step, and it is subsequently Fourier transformed. From these pseudovisibilities, we only retained those matching the observational data uv-coverage, which were then processed with imaging synthesis techniques, similarly to the observed interferometric data. This procedure ensures a direct comparison of model and real data in the image space even in the event of overesolution and lost flux problems (which are not the case here). This process is done through GILDAS’s MAPPING software, after a conversion from SHAPE’s outcome file in TXT format into a FITS format file. Both scripts (‘interferometrisation’ and TXT to FITS conversion) are available via contacting the first author.

4 M1-92 model

Guided by the overall physical and chemical characteristics described in previous studies (Solf 1994; Bujarrabal et al. 1998b,a; Davis et al. 2005; Alcolea et al. 2007; Sánchez Contreras et al. 2008; Alcolea et al. 2022), we built a model of the pre-planetary nebula located 2.5 kpc away and oriented 40° with respect to the plane of the sky. Our de-projected model is 16″ long with a cylindrical symmetry along this axis and a bipolar structure, with a radius of 2″.25 in the equator and 0″.5 at both ends of the nebula. Thanks to SHAPE’s modelling features, an accurate position of the object, both in distance and orientation, as well as precise sizes, is easily reproduced.

Our model (see Fig. 2 for model building and Table 2 for their physical properties) is made of six structures with revolution symmetry around the main axis and mirror symmetry with respect to the equatorial plane. They are divided into two categories: the outer ones, made of cold and dense gas with a shell-like shape and uniform physical properties, and the inner ones, made of warmer gas, with variable physical properties.

The largest component is the main shell, a prolate spheroid with flattened ends, which contains most of the inner (smaller) structures. This spheroid has been emptied by subtracting a similar smaller shape, leaving only a shell that will be filled with the equatorial structure, creating two empty lobes. Similar to all outer structures in our model, this shell is made of cold and dense gas, with constant values across all points, and low turbulence. Exact values of physical properties can be checked in Table 2.

A complementary structure to this is the ring. The ring is set in the equatorial area of the main shell, with 3″.5 of length and the same thickness (0″.75) and physical properties as the shell. Its only role in the model is to provide a defined region for those molecular species located only in the central part of the nebula when they reach the outer parts of its equator, but not extending along its shell.

Finally, the last outer structures are the outer parts of the tips. With very similar properties as the shell, these structures at each end of the nebula represent the continuation of that main spheroid when it squeezes into a structure 0″.5 in radius.

The largest inner structure is the central cylinder. It is also 3″.5 long and located in the nebula’s centre, with a 1″.5 radius and concave ends that complete the central shape of the empty lobes. This structure is associated with warmer but also denser gas than the outer shell, with both temperature and density laws ensuring the continuity of physical variables across the interface.

Finally, we have the two types of inner structures made of much warmer but much less dense gas. Those are the blobs and the inner part of the tips. The blobs are represented by a sphere of 0″.5 radius inside each lobe travelling along the main axis of the nebula. These are the only structures needing a much higher value for turbulence than the rest of the nebula to reproduce the velocity dispersion observed. The inner part of the tips has similar physical properties to the blobs, but this time with turbulence more similar to the rest of the nebula, although still a bit larger. The full composition of the model and each structure’s sizes are shown in Fig. 2.

For each of these structures, we set a turbulence, temperature, density, and velocity law, plus a different abundance for each molecular species present in them, in order to reproduce the observed line profiles and maps. All the structures, except the spheres, are given the same velocity law of a radial expansion proportional to the distance as $V=\frac{64.75}{7.4}$ km s⁻¹arcsec⁻¹ × r(arcsec), consistent with that detected in the nebula by previous studies. The spheres’ velocity is defined by a step function in order to match the map observations, with the 0″.7 closest to and the 0″.3 farthest from the centre of the nebula expanding radially (from the nebula’s centre) at 55 and at 15 km s⁻¹ respectively. This is the equivalent to a front-like shock distribution, which stays consistent with the extreme velocity dispersion found across different species. The remaining physical variables set in the final model are shown in Table 2, and the final relative abundances of each species used in each structure are shown in Table 3. A more detailed description of the model and its structures can be found in the SHAPE file of the model, provided among other additional material of this work.

The 3D model was run inside a grid defined by a field of view of 40″ × 40″ and a resolution of 127 × 127 cells for line profiles and 255 × 255 cells for maps. Each rendering is done by filling the model with one species at a time and for a specific transition. In the case of single-dish line profiles, the results are convolved with a circular gaussian beam 3% bigger than the nominal one of the telescope to account for tracking errors. For interferometric maps, we follow the procedure described in Section 3.3, which ensures the same clean-beam geometry as in the observations. In every case, the render mode used was HD, with a CMB set at 2.73 K and being subtracted. The spectral resolution used was 1000 spectral bands dividing a velocity range of ±70 km s⁻¹. We set the number of channels for interferometric maps to 35, through a range of ±55.25 km s⁻¹ (i.e. a velocity spacing of 3.25 km s⁻¹). We also used the newly added scripting module to compute all the transitions sequentially, allowing for a myriad of combinations of different physical variables to be computed semi-automatically. This way we could explore a wide range of values for each parameter in a reasonable time, to finally find the best-fit model. The script used, together with the SHAPE file of the model, can be found among the examples provided in the webpage.

We considered the model to be successfully built once it was able to reproduce all observations across all transitions and species using the same values for its physical properties. This means that when we add a new species to the model, the only adjustments allowed are the values of its relative abundance in each structure, with no changes in any other variable. All transitions, both from maps and line profiles, undergo a visual inspection by eye to ensure they are fit simultaneously. Previously, especially during model building, a mathematical inspection is done by calculating the emission ratio between model and observations, both for single-dish and interferometric maps. However, the final decision must be a compromise between each line’s emission distribution, across all lines of the same molecular species for the relative abundance, and across all observational data for the common physical properties, making the human eye the most powerful tool to determine the fit. This method is widely used and accepted within the research community, particularly regarding morpho-kinematical modelling of complex nebulae (e.g. Solf & Ulrich 1985; Santander-García et al. 2004, 2008; Huckvale et al. 2013), and it has lead to solid results, such as the almost-perfect alignment between the orbital plane and the equatorial waist of every post-common-envelope planetary nebulae where both orbital parameters and nebular morphology have been studied, the latter relying on the human eye as fitting tool (see Fig. 8 in Hillwig et al. 2016, and references therein).

Fig. 2 Our best-fit model in a wired representation showing the observer’s point of view and its angle with respect to the plane of the sky. Each structure is shown in a different colour and separately: green is for the outer shell, purple for the ring, red for the central cylinder, bright yellow for the spheres, and light yellow for both the outer and inner parts of the tips. The shell and ring are shown together, as the ring is only a limited region of the shell meant to define the locations of the species, but its physical properties are the same as in the shell. Values are provided for the main dimensions of each structure. Pictures of independent structures are not to the same scale. Note that 3.6×1016 cm subtends 1″.0 at the assumed distance of 2500pc.

4.1 CO model

Among our observational data, CO is the molecule that can trace the overall structure of the nebula defined by the cold and dense gas, so it is be the main source of information to build our model. We have single-dish data of five ¹²CO rotational transitions: J =1−0, J=2−1, J=5−4, J=7−6 and J=9−8. For ¹³CO, C¹⁷O and C¹⁸O we also obtained line profiles of the first two transitions, J =1−0 and J=2−1, as well as interferometric maps of the J=2−1 transition from each of them. From the interferometric data, particularly those obtained for the ¹³CO transitions, with the brightest detected emission, we obtained an accurate description of the sizes and geometry of the outer structures of the nebula. The largest and densest structures are best studied from the low-J mapped transitions, and hence this is the case for the outer shell and the central area in the nebula. In contrast, we could best determine the physical conditions from the emission of the observed ¹²CO transitions due to their larger variety of involved energies. Thus, we were able to detect emission from all structures, including the blobs, where, for instance, no CO emission has been detected in any low-J transition so far and hence for which the excitation conditions have not been properly studied until now.

From the modelling analysis, we deduced that the central cylinder has the most complex physical properties among the structures that compose it. Its temperature law is defined as an inverse function, decreasing as the radius increases until it reaches the temperature of the outer shell. A similar tendency is also set on the density, but this time with a linear decrease between a higher density centre and the outer shell. As a result, our model has a very small volume with both high density and temperature, which will make it relevant across all transitions covered by our observational data. This property distribution was chosen as a compromise to best fit the data in all observed transitions, even though it produces a peak too narrow in higher lines and excessive brightness in intermediate ones.

To achieve an accurate reproduction of the emission from all ¹²CO transitions, we needed to introduce a double-layered structure of the tips. For the emission from these structures to be clearly detectable in the J=9−8 line but absent in J=5−4 and lower, this double layer is a must. The hotter area detected in higher lines must be contained in a smaller and less dense part, while a denser but colder shell covers it. The latter ensures that while it has no relevant contribution to the lower line profiles, its presence is still clear on the maps.

The only structures for which not much information is provided by any of the CO observations are the spheres inside the lobes, introduced to reproduce the observed blobs. While their importance for the highest ¹²CO line observed is clear, they are barely noticeable in the rest of the lines, contributing only by slightly broadening the emission further along in the velocity axis due to its high turbulence. In addition, most of their velocity range overlaps with that of the tips and even the central parts of the nebula, which results in an infinite number of possible combinations to obtain the final intensity observed in the line profiles. Lacking any further information, a proper assessment of its physical conditions as well as a final definite description of the model was not possible until the HCN and HCO⁺ maps, both in J=2−1, were evaluated.

Figure 3 shows the comparison of ¹²CO line profiles from observational data (in black) and model simulations (in green). We can see the same results for the line profiles of ¹³CO, C¹⁷O and C¹⁸O in Figure 4. In all of these lines, the central peak is always the less well-fitted part of the model. Not only in intensity, the distribution of the peak differs slightly from the observations in all lines, especially in the highest-excitation ones, where the central peak of the observational data is displaced towards positive velocities. A look into the ¹³CO position-velocity (P-V) diagram in Fig. 5 shows us that the source of this inconsistency is found in the central area of the nebula.

Thanks to the self-similar (i.e. Hubble-like) expansion of the nebula, where the velocity variable is proportional to radial position, we can use these P-V diagrams to show different cuts of the nebular structure. In Figure 5 we use this property to compare the observational data and the simulated map. A clear asymmetry in the central area shows that our model lacks a bright spot in the nebula’s equator, slightly displaced towards positive velocities. Although the physical properties of the central cylinder in our model attempt to reproduce the behaviour of this spot across all lines, it is clearly not possible with the limited resolution we have of this area. As a consequence, our model produces slightly higher peaks in lower lines and too narrow peaks in higher-excitation lines. A similar pattern is found in the C¹⁷O and C¹⁸O maps as seen in Figure 6.

While we can find other asymmetries in the observational maps, most of them are also reproduced by the model after going through the ‘interferometrisation’, including thickness variations between analogue points of the shell despite its perfect symmetry, or the instrumental emission and absorption artefacts. Also, none of them are as strong as this central spot, which means this is, without any doubt, a real feature in the nebula with different physical variables than those given to the central cylinder.

As a result, the intensity of the central area becomes the main source of disagreement between both datasets, also in the line profiles. Concerning the ¹²CO line profiles, we see how the central peak of the model is overestimated in the lower transitions up to J=5−4. For higher energy transitions, however, the situation is the opposite. Even though this difference in intensity in the very centre is within uncertainties given the noise margin, it is important to notice that the peak reproduced by our model is much narrower than that observed. This shows that an additional unmodelled structure must be present in the nebula covering a wider range of central velocities, either due to its size or to the turbulence of its gas, and that it is slightly redshifted.

In all, our model simulates very accurately the observational data of CO isotopologues. Line profiles and interferometric maps helped build a robust model of the structures that dominate the emission in these lines. However, we can improve the determination of the characteristics of the nebula by also modelling the emission from other molecular species.

Fig. 3 Comparison of line profiles obtained from single-dish IRAM-30 m and HSO/HIFI observations (black) and model reproduction (green) on 12CO lines in Tmb (K) vs. LSR velocity (km s−1).

Fig. 4 Comparison of line profiles obtained from single-dish IRAM-30 m observations (black) and model reproduction (green) on 13CO, C17O, and C18O lines in Tmb (K) vs. LSR velocity (km s−1).

Fig. 5 Position-velocity diagrams of observational data from IRAM-NOEMA (left), the model (centre), and residuals (observation - model; right) from the 13CO J=2-1 map. Contours are drawn at 20 mJy beam−1 intervals. The structures that dominate the emission are annotated in the observational data panel.

Fig. 6 Position-velocity diagrams of observational data from IRAM-NOEMA (left), the model (centre), and residuals (observation - model; right) from the C17O (top panels) and C18O (bottom panels) J=2-1 maps. Contours are drawn at 5 and 2 mJy beam−1 intervals respectively.

Fig. 7 Position-velocity diagrams of observational data from IRAM-NOEMA (left), the model (centre), and residue (observation-model; right) from the HCO+ J=2-1 map. Contours are drawn at 10 mJy beam−1 intervals. The structures that dominate the emission are annotated in the observational data panel.

Fig. 8 Comparison of line profiles obtained from IRAM-30 M observations (black) and model reproduction (green) fot HCO+ lines in Tmb (K) vs. LSR velocity (km s−1).

4.2 HCO⁺

As less dense areas with higher temperatures are invisible in the lower transitions of CO maps, a proper description of the blobs’ velocity law and position is impossible without the help of HCO⁺ and HCN maps. The best adjustment for these turbulent areas, while keeping a relatively simple description, is a constant value of 18 km s⁻¹ for the turbulence, but centred on different radial velocities. The 0″.3 of the sphere farther away from the centre of the nebula are set at 15 km s⁻¹ of radial expansion velocity, while the closer part is set at 55 km s⁻¹. The accuracy of this frontlike velocity distribution, as seen by its coverage on the synthetic map, implies the presence of an ongoing shock on the fore parts between the ejected gas and previously existing material.

Besides the blobs, we can see from the maps that HCO⁺ emission is bright on the tips and central area of the nebula, but the shell also seems to have a relevant HCO⁺ abundance, particularly intense at its equator. This distribution prompted the addition of the ring structure, which only has a role different from that of the shell in the abundance of this molecule. As seen in Table 3, to properly reproduce the observations, we needed to introduce different relative abundances for the different structures, contrary to CO species, where the same abundance is used across the model. The synthetic map resulting from our model, as well as the comparison to the observational data, can be found in Figure 7.

Once again, in the data we find a bright spot slightly displaced towards positive velocities in the central part of the nebula not reproduced by our model, with a clear effect also on the line profiles as seen in Figure 8. The attempt at compensating for this emission makes the adjustment of the central peak less precise than those of the side peaks, just as it happens with ¹²CO.

The model was also applied to H¹³CO⁺ for which the same abundance distribution across structures was deduced, just applying a constant factor with respect to the values of HCO⁺. The line profiles can be seen in Figure 9, where we note how the central peak is even more dominated by the bright spot than in HCO⁺ and that it is displaced towards positive velocities.

Besides the slight inaccuracy of the central area of the nebula in most lines and maps, the overall modelling can be considered highly successful, as the physical variables remained the same from CO analysis, with only relative abundance needing to be determined for each species and structure to reproduce the map and all six line profiles. Therefore, the physical conditions used for the modelling are robust enough to be considered representative of those of the real nebula.

Fig. 9 Comparison of line profiles obtained from IRAM-30 M observations (black) and model reproduction (green) for H13CO+ lines in Tmb (K) vs. LSR velocity (km s−1).

Fig. 10 Position-velocity diagrams of observational data from IRAM-NOEMA (left), the model (centre), and residue (observation - model; right) from the HCN J=2-1 map. Contours are drawn at 5 mJy beam−1 intervals.

4.3 HCN

The same process was followed for the HCN data. Compared to HCO⁺, the HCN is less extended across the nebula, with specific bright areas in the tips, blobs and central region, and less velocity dispersion overall. In HCN, the central bright region is located around the previously mentioned bright spot, and not found at all in the rest of the equator. As for HCO⁺, the relative abundance used in the modelling needs to be different in each structure to properly reproduce the observational data, being non-zero just for the inner structures. The interferometric J=2−1 map shows each of these three structures very clearly, as seen in Figure 10, where we show P-V diagrams for data, model, and data-model residuals.

However, the model predictions for the HCN line profiles do not reproduce the data as accurately as for other species. In particular, as seen in Figure 11, there are discrepancies in reproducing the side peaks across the three lines, an issue that was not found in any of the species modelled before. As this is an issue that only affects this species, and given the differences in the covered regions seen on maps, we expect this to be more likely due to the spatial distribution of HCN and the particular conditions linked to its emission, instead of related to wrongly adopted physical conditions for the overall nebula model.

Once again, the model was applied to its ¹³C isotopic substitution, H¹³CN, but since the related observational data were dominated by noise, we have just opted for the same abundance ratio determined for HCO⁺. In Figure 12, we see that the model results are compatible with the data, considering the large data uncertainties. Only the J =1−0 transition signal slightly exceeds the noise, but considering the overestimation in its ¹²C counterpart, the real emission would most likely be below the noise level.

Fig. 11 Comparison of line profiles obtained from IRAM-30 m observations (black) and model reproduction (green) on HCN lines in Tmb (K) vs. LSR velocity (km s−1).

5 Discussion

5.1 Derived masses, kinetic energy, and scalar momentum

From our 3D, radiative-transfer-enabled morpho-kinematic model, we derived the primary values of the nebula and their distribution across the object. After integrating all the model’s cells, we obtained the following total values:

• Total mass: 0.8±0.2 M_⊙.

• Total kinetic energy: (6.5±1.3)×10⁴⁵ erg.

• Total scalar linear momentum: (4.1±0.8)×10³⁹ g·cm·s⁻¹.

We note that these values were calculated as an integration over all cells composing our model. Therefore, these results are not derived from any assumed abundance for a particular species but on the best-fit values of the physical parameters of the nebula computed over the 3D model simultaneously for every transition studied. These values agree with previous studies (Bujarrabal et al. 1998a; Alcolea et al. 2007; Li et al. 2024; Khouri et al. 2025) within a reasonable margin. The highest discrepancy is found with this last study and the value they provide for the total mass of the nebula. Khouri et al. (2025) find a total mass for M1-92 of 0.35 M_⊙, assuming a distance of 2.6 kpc, which would scale to 0.32 M_⊙ at our adopted distance. They compute this mass through an assumed temperature of 10 K and a C¹⁸O abundance with respect to hydrogen nuclei of 8.5×10⁻⁷. Although the assumed temperature is close to the one we obtain for the structures that dominate the C¹⁸O emission, their assumed abundance differs quite a bit from ours, 6×10⁻⁷ with respect to the total number of collisional particles, the majority of them assumed to be molecular hydrogen. This means their assumed abundance is ∼2.8 times higher than ours, which is close to the ratio found between the total masses of both studies after taking into account the corrections for total gas mass. Given that Khouri et al. (2025) is a survey-type study, where the same standard abundance and physical conditions are assumed for different objects, we consider our approach through detailed model fitting of observations more accurate.

A smaller but still significant difference in the total mass retrieved is found when comparing with the results from Li et al. (2024), where a total mass of 1.1 M_⊙ is reported. This study also presents the estimated mass per structure: 1.02 M_⊙ for the torus, which is the equivalent to our central cylinder and ring, 0.072 M_⊙ in the lobe walls (our shell) and 2.7 10⁻⁴ M_⊙ for each of the blobs, with no modelling of the tips. Even though the value for the total mass of the nebula is very compatible with our results, we see how the mass distribution is a strong source of disagreement when comparing these values to those presented in Table 4. Li et al. (2024) locates most of the mass in the central area, describing a high density torus and very narrow lobe walls. This distribution is valid for observations in visible wavelengths, but it is in striking contrast with radio observations, which show that these lobe walls must be considerably thicker. This is probably due to the (visible) scattered light being mostly sensitive to the internal side of these structures. We note that the total mass derived in the model by Li et al. (2024) also depends on dust properties and the gas-to-dust ratio assumed (200 for all components of the nebula), which, given the mixed chemistry found in M1-92 might not be of standard characteristics. It is important to highlight that in all of these studies, the total mass estimated is a very significant portion of the initial mass of the star, resulting in no relevant changes in the conclusions derived for the death of its central star.

We can also visualise the distribution of these values across the nebula and its components. We provide the distribution of these three values by structure in Table 4, together with the time that radiation pressure would take to provide these linear momenta. As shown in some of the previously cited studies, the overall nebula presents a momentum excess that requires a mechanism more powerful than radiation pressure to be accounted for. We note that, as for the contribution of this momentum excess per structure, the outer and denser regions are the main contributors to this excess, in contrast with the warmer structures along the main axis, which do not show such an excess. However, given the strong collimation of these later ejections, we do not assume radiation pressure to be the process behind their formation. Nevertheless, this exercise shows that the launching mechanism that happened more recently, responsible for the axial components, should be significantly less powerful than the one that formed the bulk of the nebula. It is also worth noting that the mass contained by the spheres in our model is comparable to that estimated for the H₂ emission inside the lobes by Bujarrabal et al. (1998b).

We can also plot these values as a function of latitude. In Fig. 13, we can see the total amount of scalar linear momentum and kinetic energy with each structure’s contribution in each latitude range. This figure clearly shows how material along the equator represents just a fraction of momentum and kinetic energy, whereas these quantities increase towards the polar direction and then drop at the very ends, where less material is detected. Finally, in Fig. 14, we show the distribution of temperature, density, scalar linear momentum, kinetic energy, pressure, and Mach number across the plane containing the main axis of the nebula in our model. As axis symmetry and mirror symmetry with respect to the equatorial plane are kept in our model, a single quadrant is representative of the entire nebula.

We provide these total quantities and their distribution as a detailed constraint for HD/MHD models where this plotting is common (see García-Segura et al. 2005). Unfortunately, there are no such models of M1-92 published currently, only the total values reported in Bujarrabal et al. (2001) are used when comparing with model results in García-Segura et al. (2005) and Blackman & Lucchini (2014). As our results for these values are fairly compatible with the ones in Bujarrabal et al. (2001) we conclude that our work aligns with the results from these studies, where a common envelope ejection, or similar sudden event, is needed to reproduce these momentum and energy values. Therefore, we also avoid providing the reader with mass-loss rates. Since no precise ejection timescale is estimated for this source, we cannot present a value meaningful enough for its evolutionary process.

Fig. 12 Comparison of line profiles obtained from IRAM-30 m observations (black) and model reproduction (green) on H13CN lines in Tmb (K) vs. LSR velocity (km s−1).

Fig. 13 Distribution of scalar linear momentum (left) and kinetic energy (right) across the nebula as a function of latitude angle for the different structures.

5.2 Abundances

The full final description of relative abundances per species and structure can be found in Table 3. The initial values of relative abundance with respect to collisional particles of the CO species assumed from the bibliography worked relatively well throughout the modelling process, thus requiring no significant changes. As for HCO⁺ and HCN, we find a dichotomy between outer and inner structures. For outer ones, we also concluded values within the expected range in the typical CSEs of AGBs (see Pulliam et al. 2011; Schöier et al. 2013). In contrast, the inner structures, especially the blobs and inner part of tips, present relative abundance values on the very upper limit or above. The strong presence of these molecules on the most dynamically active regions is yet another indicator of shocks (Viti et al. 2002; Jørgensen et al. 2004; James et al. 2020).

With this information on relative abundances with respect to collisional particles from our model, we derived the abundance ratios between the isotopologues of the studied species. We also estimated the error for the three main ratios as the variation needed for the difference to be significant in our model, while keeping the same physical variables.

• ¹²CO/¹³CO = 30±3;

• C¹⁷O/C¹⁸O = 1.6±0.1

• H¹²CO⁺/H¹³CO⁺ = 10±1.

Note that we can reasonably assume the same atomic isotopic ratio for ¹²C/¹³C and ¹⁷O/¹⁸O as those obtained through molecular species, since no chemical mechanism favouring certain isotopes (such as fractionation) is likely to be operating in M1-92 (see Section 5.3).The C¹⁷O/C¹⁸O ratio matches the numbers obtained by line ratios once the frequency and Einstein coefficient corrections are applied to C¹⁷O and C¹⁸O, as well as those estimated by previous studies (Alcolea et al. 2022; Khouri et al. 2025).

However, we find a strong discrepancy in the ¹²C/¹³ C ratio depending on the molecule used as a probe (CO or HCO⁺, see above). While both 30 and 10 are reasonable values for a mixed chemistry environment and not far from that previously estimated in the bibliography (see Bujarrabal et al. 1997), they are mutually incompatible in our model if applied to all the species. In Fig. 15 we show the paradox of adopting either of these two ¹²C/¹³C ratios in the adjustment of the different molecules by plotting some line profiles on the species used for probing this ratio: the optically thinnest line of ¹²CO (J=9−8), both lines of ¹³CO (J=2−1 and J =1−0), the best reproduced line of H¹²CO⁺(J =1−0), and the less noisy lines of H¹³CO⁺ ( J=2−1 and J =1−0). As can be seen in this figure, the adjustment of the CO lines fails miserably if we use the ¹²C/¹³ C ratio obtained from HCO⁺ and vice versa. We note that the difference between both values (a factor of 3) is much larger than the uncertainty due to degeneracy between abundance and density.

A ¹²C/¹³C ratio of ten was also applied to H¹³CN, as the molecule traces an area very similar to HCO⁺. One can see that it is also compatible with the observational data, as seen in Fig. 11, where the model prediction would get hidden by the noise, with the exception of the J =1−0 line already discussed in Section 4.3.

Fig. 14 Distribution across the nebula of different physical values. In the left figure, the temperature (top right), density (top left), density of scalar linear momentum (bottom right), and density of kinetic energy (bottom left) are shown. In the right figure, the pressure (top left), Mach number (top right), scalar velocity (bottom left), and turbulence (bottom right) are shown. The axes are given in physical length units (centimetres). Note that 3.6×1016 cm subtends 1″.0 at the assumed distance of 2500 pc.

Fig. 15 Comparison of model reproduction between 12C/13C ratios of 30 (blue) and 10 (magenta). The vertical axes represent Tmb (K), while the horizontal axes are LSR velocity (km s−1). Note how 12CO and 13CO lines (top panels) can only be successfully reproduced with a 12C/13C ratio of 30, while H12CO+ and H13CO+ lines (bottom panels) require a value of ten.

5.3 Physical implications

Putting all of these results together, we can gain some insights into the physical circumstances that gave birth to this nebula. With the ¹⁷O/¹⁸O abundance ratio of 1.6 obtained from CO isotopologues, we can derive an initial mass for the central postAGB star. According to most nucleosynthesis models (De Nutte et al. 2017; Karakas 2014; Cristallo et al. 2011; Stancliffe et al. 2004), this implies a stellar mass of around 1.7M_⊙ at the beginning of the main sequence if we assume a solar metallicity, which despite not being measured, it is a reasonable assumption given the location of M1-92 in the galactic disc. This means that the nebular mass as obtained through our model would be around 70% of the total mass loss expected throughout the entire AGB phase: for a 1.7 M_⊙ star we expect a white dwarf of 0.57±0.2 M_⊙ (Catalán et al. 2008; Cummings et al. 2018). This initial mass should have also been enough for the central star to become C-rich at the end of the AGB phase, as it would have experienced enough third-dredge-up events to do so (Groenewegen et al. 1995; Pardo & Cernicharo 2007; Zhang et al. 2013; Marigo et al. 2020). However, despite finding molecular species typical of C-rich environments, such as HCN, HNC or CS, the nebula shows a mixed chemistry with molecular species only found in O-rich environments, such as OH masers, SO, SO₂, and H₂O ices (see Alcolea et al. 2022; Seaquist et al. 1991; Eiroa et al. 1983). This result is compatible with previous studies such as Khouri et al. 2025, where, despite the significant difference in C¹⁸O abundance assumed in their work, the C¹⁷O/C¹⁸O ratio, which is the value that nucleosynthesis models agree the most, is virtually the same.

The amount of mass ejected, 0.79 M_⊙, relative to the initial mass, 1.7 M_⊙, together with the expansion law, suggests that the formation of the nebula was triggered by a sudden event that also interrupted the AGB evolution and nucleosynthesis (and third-dredge-up events), as this would also explain its non-C-rich chemistry. Most likely, as proposed by previous studies (Alcolea et al. 2007), this event was a common envelope ejection, as these provide the energetics and a preferential axis for the ejection (Douchin et al. 2015; Jones & Boffin 2017). We notice that this kind of pPNe, with a speculated sudden formation, is relatively common. Other examples are IRAS 16342-3814, Hen 3-1475 or OH231.8+4.2 (see Sahai et al. 2017; Khouri et al. 2025; Alcolea et al. 2001, respectively.)

In addition, the warm, turbulent gas from the inner structures along the symmetry axis, tracing shocked gas, indicates the presence of later ejections, most likely related to polar jets, which are believed to indicate a binary system at its heart (Soker 2025). If indeed we are witnessing the ejection of a common envelope, the fact that the total mass of the nebula amounts to the majority of the envelope of the AGB progenitor star would better align with the companion being a post-AGB star instead of a main sequence star (thus making the system a so-called doubledegenerate) and therefore with the observed nebula being the product of the second common envelope ejection experienced by the system, as these seem to be substantially more massive than their single-degenerate counterparts (Santander-Garcia et al. 2022).

Since it seems clear that its physical properties differ quite significantly from the rest of the equatorial area, we can speculate about the nature of the central bright spot and whether it is also a fresh new ejection or has a different origin from all the other structures. However, we are limited to the spatial resolution of our maps, from which it is impossible to determine the actual position and size of this emission. Similar asymmetries have been found in other pPNe (see Castro-Carrizo et al. 2002), without an obvious explanation for their presence. Spatial resolution is certainly the most limiting aspect of our study, since the combined data is enough to derive the existence of a doublelayered structure in the nebula, but we lack the information to eproduce a more realistic and accurate geometry in our model, with the failure in reproducing the HCN lines and the nature of the central spot as main consequences.

The main challenge we face when it comes to putting these results in context comes from the values found for the ¹²C/¹³C ratio. Although a ratio of 30 is in line with the interrupted AGB evolution via common envelope ejection hypothesis, as it is an intermediate ratio between those expected at the beginning and end of a thermally-pulsing AGB phase of a 1.7 M_⊙ star (see F.R.U.I.T.Y. database⁹ from Cristallo et al. 2016), it is clearly different from the alternate value of 10 derived from HCO⁺ and compatible with HCN. Both results seem to be very strong in our model, and the only way of reconciling them is by assuming different isotopic ratios for different structures and setting the ¹³CO relative abundance 3 times higher only in the inner, warm structures. While this overestimates the ¹³CO emission in the equatorial area, as our central cylinder is bigger than the central bright spot, it does not affect any other aspect of the ¹³CO data we have. Therefore, we cannot conclude whether this is a molecular difference for the entire nebula or an isotopic ratio difference for carbon in certain structures, particularly affecting all the chemistry happening in them. In any case, we can ensure that this is not an opacity issue, as we base our results on the model properties, plus the smaller ratio would have been favoured for the most opaque lines, i.e. those of CO instead of HCO⁺. We can also discard isotopic fractionation as no part of our nebula has physical conditions with low enough temperature for it to happen, especially the ones traced by HCO⁺ and HCN (see Roueff et al. 2015). In the same way, photodissociation is also not an option as the structure is delimited by the ballistic dynamics of the gas, with no sign of previously ejected material around it (also discussed by Li et al. 2024). Therefore, we theorise that structural difference, rather than molecular, is more likely to be the origin of these different ¹²C/¹³C values. If we assume all the central structures, including the bright spot, come from more recent ejections, this would mean the central system kept evolving, including a decrease in the ¹²C/¹³C isotopic ratio in the most recently ejected material. We do not dare to propose an exact mechanism for it, but very low (≤5) ¹²C/¹³C ratios are also found in other nebulae of explosion-like origin (see Kaminski et al. 2017, where possible mechanisms are discussed). We note that a similar result was found in the pPN CRL 618 by (Pardo & Cernicharo 2007), where a change in ¹²C/¹³C ratio is reported, with a value of ∼40 in the halo decreasing to 15 in the core. This result is supported by other studies using different tracers for the isotopic ratio (see Wyrowski et al. 2003; Wesson et al. 2010; Lee et al. 2013). According to Lee et al. (2013) this change occurred in a timespan of ∼450 years, a timescale very similar to that in M1-92 assuming the blobs as the main source of this disagreement in our study.

Fig. 16 Cuts of the nebula performed perpendicularly to its symmetry axis taken at an intermediate position of the lobes for inclinations of 30°, 40°, and 50° with respect the plane of the sky once the projected velocities had been converted into positional points along the line of sight for the 13CO J=2-1 map. In red is the fit around the data, showing the elongation or lack of.

5.4 Model limitations and sources of uncertainty

All of our best-fit physical values are supported by the strong alignment of the general characteristics of the nebula (total mass, kinetic energy, morphology, etc.) with previous studies (Solf 1994; Alcolea et al. 2007; Alcolea et al. 2022). Among those, the inclination of the symmetry axis is the parameter that shows the largest discrepancy among studies, as discussed by Li et al. (2024). In this case, our inclination of 40° with respect to the plane of the sky (50° respect to the line of sight), results in the best-fit model in which revolution symmetry along the main axis is conserved, and a single radial gradient is enough for describing the velocity field. In particular, this angle is obtained by looking at the revolution symmetry of the lobes. Provided a velocity gradient to convert velocity information along the line of sight into positional information, cuts of the data at the correct angle will result in circular (revolution symmetry) images of the nebula. In Fig. 16 we show the result of these cuts at different angles, showing elongations at other angles and a circular fit at the angle used for this modelling. We assume, at most, an error of ±5 degrees, given the data noise, but further than that, revolution symmetry could not be assumed for the nebula. However, this is still an important source of uncertainty, as the real size will be affected. In the same way, the distance to the nebula, or the relative abundance of CO with respect to the colliding particles can change our results.

However, besides the uncertainties intrinsic to the observations, such as spatial resolution, axis inclination or distance to the nebula, the modelling method still imposes further limitations. The most important limitation is the very large number of variables in this modelling process and their interdependencies, which makes it virtually impossible for us to give exact error values for each of the estimated physical variables. However, this limitation is mitigated by the extended and varied amount of data used for this work, which heavily constrains our model. To check the extent of this robustness, an exploration of the parametric space has been carried for those pairs of physical variables that can produce the most degeneracy in the model reproduction. These pairs are Abundance-Density and Density-Temperature. In addition, each of these parameters has been explored individually while keeping the rest of the best-fit values fixed. These checks are performed on the ¹²CO data, which are the ones used to define the best-fit values of the physical variables of the model, and to which the rest of abundances are referred. This exploration is also performed globally in our model, which means that the parameters are modified by the same percentage of their bestfit value in all structures at once, otherwise, the large number of combinations for each parameter and structure would make this checking process unfeasible. Due to the large number of parameters in the model, and their complex interdependences (e.g. mass and inclination with density and abundance), providing exact uncertainties for all values is far from trivial.

We find the density to have a ±10% deviation margin with respect to its best-fit value before the model reproduction becomes incompatible with the observational results, including the observational data calibration errors discussed in Section 2. However, this uncertainty becomes larger when taking into account the degeneracy with the abundance parameter, which can compensate for a change in density up to a 20%. We note that this is not the only source of uncertainty for the density, where the volume of the structure and therefore the spatial resolution of maps also play a role. We estimate an uncertainty of 30% of the ratio of the beam size to each structure size, which adds an extra 15% for the largest structures and 30% for the smallest to the density uncertainty. All of these errors are included in the uncertainty provided in Table 2.

For temperature the individual margin is of ±8%. We note that temperature is very constrained in the structures where ¹²CO is optically thick, discarding a possible degeneracy of the temperature-density pairing. It is true that optically thin structures could present certain degeneracy; however, even the higher lines show a strong presence of optically thick structures ( J=9−8 line has a significant contribution of the central cylinder to the line profile), making the exploration of individual structures difficult to evaluate and overall inconsistent with the checking method used.

The ±20% error on the abundance-density combination translates also into the ¹²CO abundance and therefore to all other abundance values. This error, together with the individual errors accounting also for data calibration, are presented in Table 3. We highlight that this degeneracy on the abundance-density combination does not affect the abundance ratios since a change of density would require proportional changes in all other abundances. In the same way, the data calibration does not have any effect on the isotopic ratio values as observations (both singledish and interferometric maps) of same molecule isotopologues where observed simultaneously. In the cases where the bandwidth was not enough to observe them at the same time, a third transition from another species common to both setups acts as calibrator. This way, the relative calibration is the same for those transitions, allowing the ratios to be as accurate as possible by minimising the relative error.

Regarding velocity, the general velocity law is also adopted from previous studies, and does not impose further uncertainty than that of the spectral resolution (3.25 km s⁻¹). This does create a degeneracy in the possible velocity law of slower parts, i.e. central area of the nebula. The only structures whose velocity law has been determined throughout the modelling process are the spheres, where, due to their small size, the spatial resolution introduces an error in the determination of the shock-front limit of ±0.″1, as inferred from the emission distribution in the different velocity channels of the maps.

All of these modelling errors and instrument contributions are also taken into account when presenting the error margins of the derived masses, kinetic energy and scalar momentum. However, not all sources of errors affect these final values. For instance, the error on the total density value derived from geometrical uncertainties cancels out with volume uncertainties when computing the total mass.

6 Conclusions

The new observational data presented in this paper provides the most complete and up to date view of M1-92 and consequently motivated our detailed study and analysis of each structure of the nebula. We have presented a significant upgrade of the shapemol code together with minor updates on the SHAPE software as a whole. This upgrade provides a huge improvement to the capabilities of this modelling tool, expanding the previous range of physical conditions for ¹²CO and ¹³CO and adding ten new molecular species. The use of radiative transfer and LVG approximation provide an accurate reproduction of the observational data in many physical conditions, including optically thick cases. This sets SHAPE to the standard of top radio telescopes and interferometers in the world, which are currently detecting and mapping not just CO but many other molecules of great impact in the study of astrophysical objects.

We brought this new update to the real case of M1-92, making a complete morpho-kinematical study of the nebula. This is the most complete study of this object to date, collecting 23 line profiles and five interferometric maps to be reproduced simultaneously and giving strong constraints to the physical variables. As a result, we derived for the first time the existence of a two-layer structure along the preferential axis of the nebula together with its turbulent and shock-front features. In addition, we were able to study the chemical characteristics of the nebula, finding over-abundances of certain molecules and a discrepancy in the ¹²C/¹³C isotopic ratio depending on the molecule used as a probe. Lastly, we present our results as a global distribution of the key physical variables across the nebula and as a function of latitude, resembling the products of HD/MHD models, hoping that in this way, comparison with these models will be more immediate.

All of these findings shed light on the origin of this nebula and, by extension, of the formation of other pPNe of this kind. Our analysis shows that the central star had its thermally pulsing AGB phase interrupted by an event that triggered the main nebular ejection and that the ejection happened at once. However, there are still many features that raise new questions about the formation and evolution of the nebula, including the apparent ¹²C/¹³C ratio change during the post-AGB evolution and its link with the later collimated ejections as well as the nature of the central bright spot. We hope to answer these questions in future works, provided we obtain more spatially detailed observations.

Data availability

The electronic material produced in this study, i.e. shapemol tables and SHAPE model and script, can be found in https://doi.org/10.5281/zenodo.19473007.

Acknowledgements

This work is part of the coordinated research project NEBULAe WEB, grants PID2019-105203GB-C21 and PID2019-105203GB-C22, funded by MICIU/AEI/10.13039/501100011033, and coordinated project MESON, grants PID2023-146056NB-C21 and PID2023-146056NB-C22, funded by MICIU/AEI/10.13039/501100011033 and ERDF/EU. E.M. also thanks support from the IGN (MITMS, Spain) fellowship programme. This paper has made use of data obtained with the IRAM 30 m-MRT and NOEMA telescopes, and HSO/HIFI. The observations were carried out under project number W17BJ with the IRAM NOEMA Interferometer, projects number 050-15, 160-15 and 047-16 with the IRAM 30 m telescope and projects GT1_dteyssie_1 and OT1_vbujarra_4 with the HSO/HIFI telescope. IRAM, the Institut de radioastronomie millimétrique, is an international research institute funded by the French Centre National de la Recherche Scientifique (CNRS), the German Max-Planck Gesellschaft (MPG) and the Spanish Instituto Geográfico Nacional (IGN). HSO, the Herschel Space Observatory, is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. HIFI, the Heterodyne Instrument for the Far Infrared, has been designed and built by a consortium of institutes and university departments from across Europe, Canada, and the United States under the leadership of SRON (Netherlands Institute for Space Research), Groningen, The Netherlands, and with major contributions from Germany, France, and the US. Consortium members are: Canada: CSA, U.Waterloo; France: CESR, LAB, LERMA, IRAM; Germany: KOSMA, MPIfR, MPS; Ireland, NUI Maynooth; Italy: ASI, IFSI-INAF, Osservatorio Astrofisico di Arcetri-INAF; The Netherlands: SRON, TUD; Poland: CAMK, CBK; Spain: Observatorio Astronmico Nacional (IGN), Centro de Astrobiologa (CSIC-INTA). Sweden: Chalmers University of Technology-MC2, R SS & GARD; Onsala Space Observatory; Swedish National Space Board, Stockholm University-Stockholm Observatory; Switzerland: ETH Zurich, FHNW; USA: Caltech, J.P.L., NHSC. This work has made use of GILDAS software for data reduction and analysis (https://www.iram.fr/IRAMFR/GILDAS). We also acknowledge the work by Ismael Domínguez Salamanca during his Master’s Final Project.

References

Appendix A Additional tabular material

Appendix B Interferometric maps

Fig. B.1 Interferometric maps obtained with IRAM NOEMA for the 13CO J=2-1 transition for the different 3.25 km s−1 -wide velocity channels. The central LSRK velocity (in km s−1 ) for each channel is shown in the upper left corner of the corresponding panel. Note that these maps have been rotated 49° clockwise so the symmetry axis of the nebula is in the vertical direction. The Field of view is 6″.5×15″.0. Clean-beam size is 0″.51×07 78, with a PA of 76°, and is displayed in the inset at the lower-left corner. Contours are drawn at 20 mJy beam−1 intervals. The false-colour intensity scale displayed in Jy beam−1 is also indicated right of the last panel.

Fig. B.2 Same as Fig. B.1 for the C17O J=2-1 transition and 2.5 mJy beam−1 intervals for contours. Clean-beam size is 0″.51×0″.76, with a PA of 76°.

Fig. B.3 Same as Fig. B.1 for the C18O J=2-1 transition and 2.5 mJy beam−1 intervals for contours. Clean-beam size is 0″.51×0″.78, with a PA of 76°.

Fig. B.4 Same as Fig. B.1 for the HCO+ J=2-1 transition and 10 mJy beam−1 intervals for contours. Clean-beam size is 0″.53×0″.76, with a PA of 86°.

Fig. B.5 Same as Fig. B.1 for the HCN J=2-1 transition and 5 mJy beam−1 intervals for contours. Clean-beam size is 0″.53×0″.77, with a PA of 86°.

Fig. B.6 Moment zero maps (top panels) and integrated emission profiles (bottom panels) for the lines observed with the NOEMA interferometer: from left to right, 13CO J=2-1, C17O J=2-1, C18O J=2-1, HCO+ J=2-1, and HCN J=2-1. For the moment zero maps, the intensity is given in Jy beam−1. Contours are drawn at 10% to 90% by intervals of 10% of the maximum in the respective map. The contour step, 10% of maximum, in Jy beam−1, is given at the bottom of the panel. The colour scale is fixed between −10% and 100% of the maximum emission. Offsets are given along the major symmetry axis of the nebula (the maps have been rotated by 49° so the revolution axis of the nebula, at PA 311° is in the vertical direction). Integrated emission profile fluxes are given in Jy vs. LSR velocities in km s−1. To keep the same flux scale in all five panels, spectra other than that of 13CO J=2-1 have been multiplied by the factor indicated in the top right corner.

All Tables

All Figures

	Fig. 1 Evolution of the error as opacity decreases, calculated as the discrepancy between shapemol and the theoretical calculations with respect to the theoretical value. Only the opacity conditions where the molecule had thermalised were used.
In the text

Fig. 2 Our best-fit model in a wired representation showing the observer’s point of view and its angle with respect to the plane of the sky. Each structure is shown in a different colour and separately: green is for the outer shell, purple for the ring, red for the central cylinder, bright yellow for the spheres, and light yellow for both the outer and inner parts of the tips. The shell and ring are shown together, as the ring is only a limited region of the shell meant to define the locations of the species, but its physical properties are the same as in the shell. Values are provided for the main dimensions of each structure. Pictures of independent structures are not to the same scale. Note that 3.6×1016 cm subtends 1″.0 at the assumed distance of 2500pc.

In the text

	Fig. 3 Comparison of line profiles obtained from single-dish IRAM-30 m and HSO/HIFI observations (black) and model reproduction (green) on 12CO lines in Tmb (K) vs. LSR velocity (km s−1).
In the text

	Fig. 4 Comparison of line profiles obtained from single-dish IRAM-30 m observations (black) and model reproduction (green) on 13CO, C17O, and C18O lines in Tmb (K) vs. LSR velocity (km s−1).
In the text

	Fig. 5 Position-velocity diagrams of observational data from IRAM-NOEMA (left), the model (centre), and residuals (observation - model; right) from the 13CO J=2-1 map. Contours are drawn at 20 mJy beam−1 intervals. The structures that dominate the emission are annotated in the observational data panel.
In the text

	Fig. 6 Position-velocity diagrams of observational data from IRAM-NOEMA (left), the model (centre), and residuals (observation - model; right) from the C17O (top panels) and C18O (bottom panels) J=2-1 maps. Contours are drawn at 5 and 2 mJy beam−1 intervals respectively.
In the text

	Fig. 7 Position-velocity diagrams of observational data from IRAM-NOEMA (left), the model (centre), and residue (observation-model; right) from the HCO+ J=2-1 map. Contours are drawn at 10 mJy beam−1 intervals. The structures that dominate the emission are annotated in the observational data panel.
In the text

	Fig. 8 Comparison of line profiles obtained from IRAM-30 M observations (black) and model reproduction (green) fot HCO+ lines in Tmb (K) vs. LSR velocity (km s−1).
In the text

	Fig. 9 Comparison of line profiles obtained from IRAM-30 M observations (black) and model reproduction (green) for H13CO+ lines in Tmb (K) vs. LSR velocity (km s−1).
In the text

	Fig. 10 Position-velocity diagrams of observational data from IRAM-NOEMA (left), the model (centre), and residue (observation - model; right) from the HCN J=2-1 map. Contours are drawn at 5 mJy beam−1 intervals.
In the text

	Fig. 11 Comparison of line profiles obtained from IRAM-30 m observations (black) and model reproduction (green) on HCN lines in Tmb (K) vs. LSR velocity (km s−1).
In the text

	Fig. 12 Comparison of line profiles obtained from IRAM-30 m observations (black) and model reproduction (green) on H13CN lines in Tmb (K) vs. LSR velocity (km s−1).
In the text

	Fig. 13 Distribution of scalar linear momentum (left) and kinetic energy (right) across the nebula as a function of latitude angle for the different structures.
In the text

Fig. 14 Distribution across the nebula of different physical values. In the left figure, the temperature (top right), density (top left), density of scalar linear momentum (bottom right), and density of kinetic energy (bottom left) are shown. In the right figure, the pressure (top left), Mach number (top right), scalar velocity (bottom left), and turbulence (bottom right) are shown. The axes are given in physical length units (centimetres). Note that 3.6×1016 cm subtends 1″.0 at the assumed distance of 2500 pc.

In the text

	Fig. 15 Comparison of model reproduction between 12C/13C ratios of 30 (blue) and 10 (magenta). The vertical axes represent Tmb (K), while the horizontal axes are LSR velocity (km s−1). Note how 12CO and 13CO lines (top panels) can only be successfully reproduced with a 12C/13C ratio of 30, while H12CO+ and H13CO+ lines (bottom panels) require a value of ten.
In the text

	Fig. 16 Cuts of the nebula performed perpendicularly to its symmetry axis taken at an intermediate position of the lobes for inclinations of 30°, 40°, and 50° with respect the plane of the sky once the projected velocities had been converted into positional points along the line of sight for the 13CO J=2-1 map. In red is the fit around the data, showing the elongation or lack of.
In the text

Fig. B.1 Interferometric maps obtained with IRAM NOEMA for the 13CO J=2-1 transition for the different 3.25 km s−1 -wide velocity channels. The central LSRK velocity (in km s−1 ) for each channel is shown in the upper left corner of the corresponding panel. Note that these maps have been rotated 49° clockwise so the symmetry axis of the nebula is in the vertical direction. The Field of view is 6″.5×15″.0. Clean-beam size is 0″.51×07 78, with a PA of 76°, and is displayed in the inset at the lower-left corner. Contours are drawn at 20 mJy beam−1 intervals. The false-colour intensity scale displayed in Jy beam−1 is also indicated right of the last panel.

In the text

	Fig. B.2 Same as Fig. B.1 for the C17O J=2-1 transition and 2.5 mJy beam−1 intervals for contours. Clean-beam size is 0″.51×0″.76, with a PA of 76°.
In the text

	Fig. B.3 Same as Fig. B.1 for the C18O J=2-1 transition and 2.5 mJy beam−1 intervals for contours. Clean-beam size is 0″.51×0″.78, with a PA of 76°.
In the text

	Fig. B.4 Same as Fig. B.1 for the HCO+ J=2-1 transition and 10 mJy beam−1 intervals for contours. Clean-beam size is 0″.53×0″.76, with a PA of 86°.
In the text

	Fig. B.5 Same as Fig. B.1 for the HCN J=2-1 transition and 5 mJy beam−1 intervals for contours. Clean-beam size is 0″.53×0″.77, with a PA of 86°.
In the text

Fig. B.6 Moment zero maps (top panels) and integrated emission profiles (bottom panels) for the lines observed with the NOEMA interferometer: from left to right, 13CO J=2-1, C17O J=2-1, C18O J=2-1, HCO+ J=2-1, and HCN J=2-1. For the moment zero maps, the intensity is given in Jy beam−1. Contours are drawn at 10% to 90% by intervals of 10% of the maximum in the respective map. The contour step, 10% of maximum, in Jy beam−1, is given at the bottom of the panel. The colour scale is fixed between −10% and 100% of the maximum emission. Offsets are given along the major symmetry axis of the nebula (the maps have been rotated by 49° so the revolution axis of the nebula, at PA 311° is in the vertical direction). Integrated emission profile fluxes are given in Jy vs. LSR velocities in km s−1. To keep the same flux scale in all five panels, spectra other than that of 13CO J=2-1 have been multiplied by the factor indicated in the top right corner.

In the text