© ESO, 2014

1. Introduction

Radioastronomical instrumentation has undergone significant improvement in recent years, leading to notable advances in our understanding of the Universe at molecular wavelengths. On the one hand, the spatial resolution achievable with the Atacama Large Millimeter/submillimeter Array (ALMA) rivals that of the Hubble Space Telescope in the optical regime, with unprecedented sensitivities. On the other hand, the Heterodyne Instrument for the Far Infrared aboard Herschel (HIFI) is an invaluable tool as it has effectively opened a new window from which to probe warm molecular gas (~50–1000 K). While incapable of generating images of the source because of its low spatial resolution, and thus lacking information on the spatial distribution of its molecular gas, HIFI produced 1D high-resolution spectra in the range of 480−1250 GHz and 1410−1910 GHz, which includes transitions as high as ¹²CO and ¹³CO J = 16−15, unobservable from the ground – as is most of the sub-mm and far-infrared range. In operation until April 2013, Herschel/HIFI has contributed to the study of the excitation conditions of the warm molecular gas of many astrophysical nebulae (e.g. Bujarrabal et al. 2012).

The study of nebulae around evolved stars, where molecular gas is very common, can greatly benefit from the power of modern radioastronomical instrumentation, such as ALMA’s spatial resolution and sensitivity and the ability of Herschel/HIFI to probe warm molecular gas. Being aware of the growing importance of the development of tools for simplifying the analyses of molecular data from new-era observatories, we introduce the computer code shapemol¹. This complementary software enhances the ability of SHAPE, the often-used tool for performing spatio-kinematical modelling in the field of protoplanetary (PPNe) and planetary nebulae (PNe), by allowing the study of the excitation conditions of nebulae from the analysis of CO lines, which are known to be the best tracers of molecule-rich gas in these objects. Together, SHAPE+shapemol can generate high-resolution synthetic spectral profiles and maps to be compared with observational data. Although SHAPE+shapemol are primarily intended for the study of PPNe and PNe, any other astrophysical nebulae that fulfill some general conditions (see Sect. 3.1) can be studied with this software. shapemol has been succesfully tested in the study of the excitation conditions of the molecular envelope of the planetary nebula NGC 7027 using data from Herschel/HIFI and IRAM 30 m (Santander-García et al. 2012), and is applied here to the molecular envelope of the young PNe NGC 6302, showing the power of shapemol and yielding an accurate description of the complex structure and physical conditions of its molecular envelope.

The present work is structured as follows: We provide a brief description of SHAPE and the motivation for shapemol in Sect. 2. Section 3 comprehensively describes this complementary software, including both its treatment of the radiative transfer, based on the well-known large velocity gradient (LVG) approximation, and the implementation of the software itself. Section 4 outlines the uncertainties to be expected in the synthetic results, and Sect. 5 covers the limitations of shapemol. Finally, the science case is presented in Sect. 6 with the application of SHAPE+shapemol to the molecular envelope of NGC 6302.

2. SHAPE

SHAPE is an interactive 3D software tool for modelling complex gaseous nebulae (mainly planetary nebulae, but also supernova remnants, light echoes, emission nebulae from massive stars, high-energy phenomena, etc.). The distribution of density, velocity, and other physical properties is generated interactively using 3D mesh structures and other graphical and mathematical tools. From these data the program generates synthetic images, position-velocity diagrams, 1D spectral profiles, and channel maps for direct comparison with observations. Its versatility has made it a standard tool for the 3D reconstruction of planetary nebula (e.g. Steffen et al. 2011; Jones et al. 2010) and the analysis of hydrodynamical simulations (e.g. Steffen et al. 2009; Velázquez et al. 2011). SHAPE implements radiative transfer solving for atomic species using coefficients from the CHIANTI (Landi et al. 2012), Kurucz (Smith et al. 1996), and NIST (Reader et al. 2012) databases. However, molecular physics in thermalised and non-thermalised cases was not implemented in SHAPE until now. We designed shapemol to fill this gap.

While an earlier version of shapemol worked as a complement to SHAPE v4.5, it has been fully integrated into SHAPE v5. In its present state, shapemol enables radiative transfer in ¹²CO and ¹³CO lines. This is done by interpolating the absorption and emission coefficients from a set of pre-generated tables computed under the assumption of the LVG approximation.

This choice of approximate solutions over exact ones reflects the philosophy behind shapemol: while a typical 3D code with 10⁴ grid cells would take several minutes to compute the exact solution in current regular computers, SHAPE+shapemol will handle models consisting of 10⁶ grid-cells in a matter of 30−40 s. This, together with its integration within the user-friendly SHAPE software, already a standard in the field, makes shapemol a quick and versatile tool for analysing molecular data, as it allows the generation of almost any 3D structure and velocity fields with a few mouse clicks.

In the next section, we describe the specific treatment of the radiative transfer behind this complementary software as well as its implementation and interaction with SHAPE.

3. SHAPEMOL

shapemol produces synthetic spectral profiles and maps of ¹²CO and ¹³CO lines to be compared with sub-mm and mm-range observations. In particular, shapemol computes the absorption and emission coefficients k_ν and j_ν for a large number of cells defined within the nebula, as a function of a set of parameters given by the SHAPE model (i.e. positions, velocities, densities, and temperature of each point in the nebula) together with additional input (i.e. abundances, desired transition, and species) from the user. SHAPE then uses the computed coefficients for solving the radiative transfer equation, generating synthetic spectral profiles, and/or maps of the nebula at the given transition.

shapemol is included in SHAPE v5.0, although the tables on which it is based (see Sect. 3.2) have to be downloaded separately². The tables are freely available, the only requisite being that in publications that include models prepared with shapemol, credit should be given citing this work (as well as citing the paper describing SHAPE by Steffen et al. 2011).

3.1. Treatment of the radiative transfer and level population

shapemol solves the standard radiative transfer equation to calculate the line intensity along a large number of lines of sight and many frequencies around the resonant line: $\mathrm{d}{\mathit{I}}_{\mathit{\nu }}\mathit{/}\mathrm{d}\mathit{l}\mathrm{=}{\mathit{k}}_{\mathit{\nu }}{\mathit{I}}_{\mathit{\nu }}\mathrm{+}{\mathit{j}}_{\mathit{\nu }}\mathit{,}$(1)where k_ν and j_ν are the absorption and emission coefficients of a given point in the nebula, which depend on the populations of the levels joined by the transition. We assume that the frequency dependence of the normalized line profile φ_ν is the same for both coefficients, k_ν = kφ_ν, j_ν = jφ_ν. This total redistribution assumption is justified in our case and particularly for CO lines because of the high probability of the elastic collisions and the relatively low radiative probabilities. In the case of IR and mm-wave transitions from diffuse media (interstellar or circumstellar clouds), φ_ν is always given by the local velocity dispersion, that is, the thermal or microturbulent velocities. The (non-LTE) level populations used to solve the transfer equation are calculated using the well-known LVG approximation.

The level populations depend on the collisional transition rates and the radiative excitation and de-excitation rates, which in turn depend on the amount of radiation reaching the nebula point we are considering at the frequency of the line (averaged over angle and frequency within the local line profile). The calculation of this averaged radiation intensity requires a previous knowledge of the absorption and emission coefficients in the whole cloud to solve the radiative transfer equation in all directions and frequencies. Indeed, the populations of very many levels must be calculated simultaneously, since the population of each one depends on those of the others, and in all the points of the cloud. This renders the solution of the system extremely complex in the general case.

The problem is greatly simplified when there is an LVG in the cloud, introducing significant Doppler shifts between points that are sufficiently far away. When this shift is larger than the local velocity dispersions, the points cannot radiatively interact at large scales, so the radiative transfer is basically a local phenomenon. The interaction of molecules and radiation at a given point can be described by the well-known escape probability. The excitation equations can then be solved and k_ν and j_ν calculated at each point, independently of the rest of the cloud, in the frame of the LVG approximation. In any case, the level populations depend in a complex way on the (local) physical conditions, and the solution requires an iterative process. See Castor (1970) for a comprehensive version of the general LVG formalism. The LVG approximation includes the main ingredients of the problem (collisional and radiative rates, trapping when opacities are high, population transfer between different levels, etc.) and gives fast, sensible solutions. These excitation calculations are quite accurate, even when the LVG conditions are barely satisfied, particularly in the case of AGB and post-AGB shells; see detailed comparison with non-local exact calculations in Bujarrabal & Alcolea (2013), and further use of the LVG formalism by Teyssier et al. (2006) and Ramstedt et al. (2008). The approximation itself is not necessary to derive the resulting line profiles, which can be calculated solving the standard, full transfer equation using the level populations derived from the LVG approximation and a local velocity dispersion, as indeed we have done using SHAPE. The use of the LVG approximation to calculate the absorption and emission coefficients is basic to keep the simple and fast use of shapemol in very general 3D cases, since in this approximation the coefficients only depend on a few local parameters, and not on the general structure of the nebula.

According to the LVG approximation, k_ν and j_ν depend in a heavily non-linear way on the density n and the temperature T of each point, and almost linearly on the product $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$, where r is the distance of a given point (or grid cell) to the central star, V its velocity, and X the abundance of the species. In addition to these three parameters, the results of the LVG depend on the logarithmic velocity gradient, ϵ = dV/dr r/V, but only slightly (see Figs. 1 and 2). In its current state, shapemol lets the user select values of ϵ = 0.2, 1 and 3. The calculation of the table values are simplified in the most common case, when ϵ = 1, that is, a linear dependence of V on r, since the escape probability of a photon in that case is given by a simple analytical function of the opacity. Such velocity fields are found in many young planetary nebulae, which are basically expanding at high velocity following a “ballistic” velocity law (i.e. with a linear dependence on the distance from the central star). We stress that the LVG results (in the standard case) only depend on four parameters: n, T, $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$, and ϵ.

In our calculations, we used the most recent collisional rates from the LAMDA database³, Schöier et al. (2005), and rates calculated by Yang et al. (2010) for collisions with ortho- and para-H₂, for which we assumed an abundance ratio of 3. Collisions with other gas components, notably He, are not considered, so the derived density represents the total density and not only the density of H₂.

Absorption and emission coefficients were obtained for the 17 lowest rotational transitions of the ground-vibrational state, practically the only ones observed. In calculating them, a total of 40 rotational levels were considered, although in low-excitation cases only the lowest ones are significantly populated. We verified the convergence of the calculated coefficients on the total number of levels, which was expected because we only considered temperatures well below the energies of the highest considered levels. Except for tests of the calculation accuracy, we only considered the rotational transitions in the ground-vibrational state, since vibrational excitation has been found to have negligible effects in practical cases, see Appendix A. In fact, calculations including vibrational excitation are easy to perform (we just have to include more levels and transitions in the code), but this limitation is important for our calculations to be reasonably general. Otherwise, IR continuum sources in the studied nebulae, which cause most vibrational transitions, should be taken into account, resulting in line excitations that would depend in a complex way on the structure and intensity of IR sources, on the distance between them and the considered points, and on the absorption of IR photons in their trajectory between both. The general-purpose design of shapemol, as well as its user-friendliness and fast calculation, would be strongly limited in that case.

As we see in Appendix A, the negligible contribution of the IR excitation in fact holds in regions with significant emission. When the lines are very weak, for instance because of a very low density, the contribution of IR excitation can be important compared with the emission if IR is not included, but remaining very low in any case. In practical cases, that is, whenever the lines are detectable, the contribution of IR excitation to the total emission is very low.

These continuum sources also emit at rotational-line frequencies. In the current version of shapemol, the absorption by rotational transitions of photons from the cosmic background and the rotational transitions themselves are taken into account, but not those from the local continuum sources. In our calculations in Appendix A, we take into account the effects of this absorption on the molecular level population and line formation. We have checked that these effects are even lower than those due to vibrational excitation by IR radiation.

We recall that vibrational excitation can be important in other sources, for instance for high-J emission from circumstellar envelopes around AGB stars under certain conditions, because in these objects the densities (and therefore the collisional excitation) are lower and the continuum IR emission is higher. We finally note that for molecules other than CO, the vibrational excitation can be more intense, because carbon monoxide is easily excited by collisions (and its levels thermalised).

Fig. 1 Examples of the jν coefficient corresponding to the 12CO transitions with J = 2–1, 6–5, 10–9 and 16–15 as a function of temperature in six different physical conditions, with varying ϵ, n and $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$. Some of the transitions have been scaled up or down by the specified factors for displaying purposes. The non-linear dependance on the density n can easily be seen here, as well as the almost linear dependance on the product $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$. The value of ϵ, on the other hand, affects the jν coefficient only slightly.

Fig. 2 Examples of the kν coefficient corresponding to the 12CO transitions with J = 2–1, 6–5, 10–9 and 16–15 as a function of temperature in six different physical conditions, with varying ϵ, n and $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$. Some of the transitions have been scaled up or down by the specified factors for displaying purposes. The non-linear dependance on the density n can easily be seen here, as well as the almost linear dependance on the product $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$. The value of ϵ, on the other hand, affects the kν coefficient only slightly.

3.2. Implementation of shapemol

In principle, computing the k_ν and j_ν coefficients for each point of the model would require solving the population level equations once for the physical conditions of each single point via the convergence algorithm typical of LVG codes. In practice, however, this is a time-consuming task in a model typically consisting of a few millions of sampled points. Instead, the approach of shapemol consists of relying on a set of pre-generated tables of k_ν and j_ν as functions of n and T, each table corresponding to a species (¹²CO and ¹³CO in its present state), a value of ϵ, and a value of the $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$ product. The values of k_ν and j_ν in each of these tables have been previously computed for each n and T by iteratively solving the population level equations.

The flow diagram of shapemol is shown in Fig. 3. The user selects the desired species and transition in the SHAPE interface, together with abundance X, ϵ, and δ_V (the characteristic microturbulence) values of each structure of the model (lobes, polar jets, etc.). The code then samples the positions, velocities, densities, and temperatures of the whole model – with a spatial resolution chosen by the user –, and computes the product $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$ for each sampled position. Then, it selects the table with the closest value of ϵ and $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$ for the given species and transition. Once a table has been selected, shapemol computes the k_ν and j_ν for each sampled position by linear interpolation between the values for the two adjacent tabulated values of n and T. The steps in n and T were specifically chosen to be small enough so as to guarantee that a first order interpolation is a good approximation between two consecutive values (see Sect. 4). Finally, given the roughly linear dependence of k_ν and j_ν on the product $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$, the software scales the computed absorption and emission coefficients according to the ratio of the sampled position value of $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$ to that of the selected table; to avoid significant errors, calculations were performed for many values of this parameter.

By default, shapemol computes k_ν and j_ν by interpolation and applies a linear correction on the coefficients based on the ratio of the sampled position value of $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$ to that of the selected table. If, however, the option “Linear Interpolation” is set, shapemol selects the two tables whose $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$ is closest to that of the cell, linearly interpolates the values of k_ν and j_ν in both tables in order to have the coefficients at those two $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$, and uses those to linearly interpolate the final k_ν and j_ν at the cell’s $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$. This procedure is slightly more accurate, since no assumptions are made on the dependence of k_ν and j_ν on $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$, although it is slightly slower.

In its present state, shapemol is able to compute the J = 1 − 0, 2−1, 3−2,... up to the 17−16 transitions of ¹²CO and ¹³CO. There are 199 values of T and 21 values of n in the available pre-generated tables, ranging from 10 to 1000 K in steps of 5 K for the temperature, and from 10² cm^-3 to 10⁷ cm^-3 in multiplicative factors of ${\mathrm{4}}^{}\sqrt{\mathrm{10}}$ for the density. shapemol has been successfully used to study the excitation conditions of the molecular gas of the young planetary nebula NGC 7027, which was found to show a relatively complex physical conditions and structure (Santander-García et al. 2012).

3.3. Other molecular species

Future plans include the addition of more molecules, such as C¹⁸O, CS, and SiO, whenever the effect of vibrational excitation is low enough to warrant results as accurate as those presented here for CO (see Appendix A). In any case, the implementation of shapemol allows the inclusion of user-supplied sets of tables containing the absorption and emission coefficients of transitions of other molecular species. When a set of tables in a specific format has been placed in the user’s computer and the index text file has been properly modified, shapemol will be able to work with the supplied data without the need for compiling a new version of SHAPE. For questions or details on the format, contact the first author of this work.

For simple molecules and when the IR source can be assumed to be small, located in the centre of the nebula, and emitting isotropically, the effect of the IR vibrational excitation could be implemented, at a cost of a considerable increase in the volume of the tables (of the order of at least a factor 20) and a detailed study of the applicability of the treatment itself.

Fig. 3 Workflow of shapemol.

3.4. Interaction with SHAPE

This is a brief, conceptual description of the workflow when using SHAPE together with shapemol.

• 1.

  We model a nebula with several distinct outflows or structures using the standard tools in SHAPE. The nebula must have defined values of the density, temperature, and velocity modifiers. Next, we create the desired molecular species – ¹²CO or ¹³CO – and transition for each structure of the nebula, and assign values for the molecular abundance, ϵ, and micro-turbulence, δ_V.

• 2.

  When the model is run, the nebula is then sampled, point by point, with a resolution chosen by the user. Each of these points is characterised by an identification code (unique for each structure), a position, a velocity, a density, and a temperature. Only points that are actually populated are processed explicitly, this saves memory and processing time. shapemol then computes the absorption and emission coefficients using the LVG approximation-based tables.

• 3.

  SHAPE reads the output by shapemol and performs radiative transfer solving for the whole nebula, taking into account the emission from the cosmic microwave background, which is subtracted afterwards in order to compute the final line profile. The result is a data-cube consisting of a stack of images of the nebula, one per frequency resolution element (i.e. spectral channel). This data-cube is then convolved, channel by channel, with the telescope beam, which is simulated by a Gaussian with a Full Width at Half Maximum (FWHM) equal to the telescope Half-Power Beam Width (HPBW). The profile from the grid pixel corresponding to the telescope pointing is taken, and its intensity is divided by the projected area of the grid pixel on the plane of the sky. A final profile is then shown, in units of W s^-1 Hz^-1 m^-2 sr^-1 or main beam temperature (T_mb), allowing for direct comparison with observations.

4. Sources of uncertainty

4.1. Radiative transfer solving in SHAPE

Radiative transfer solving requires finding the intensity, for each frequency, that matches the following equation at each location of the target (i.e. a nebula) in the plane of the sky: ${\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{=}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{s}}_{\mathrm{0}}\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }\mathrm{\left(}\mathit{s}\mathrm{\right)}}\mathrm{+}{\mathrm{\int }}_{{\mathit{s}}_{\mathrm{0}}}^{\mathit{s}}{\mathit{S}}_{\mathit{\nu }}{\mathrm{e}}^{\mathrm{\left[}\mathit{\tau }\mathrm{\right(}\mathit{s}\mathrm{)}\mathrm{-}\mathit{\tau }\mathrm{(}\mathit{x}\mathrm{\left)}\mathrm{\right]}}\mathrm{d}\mathit{\tau }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{,}$(2)where s₀ is the position of the edge of the target located farthest away from the observer, s is the position measured from s₀ along the line of sight towards the observer, τ is the optical depth, with its differential dτ(x) and x the integration variable, also running along the line of sight towards the observer, I is the intensity emitted per unit of time, frequency, and solid angle, and S the source function, which is the ratio of the emission to the absorption coefficient: $\mathit{S}\mathrm{=}\frac{{\mathit{j}}_{\mathit{\nu }}}{{\mathit{k}}_{\mathit{\nu }}}$. In general, S varies along dτ(x) and thus the solution to the equation is often complex, and can seldom be expressed in analytical form.

The approach that SHAPE uses consists of sampling the nebula point by point according to the resolution chosen by the user. Provided that we know the absorption and emission coefficients (and therefore the source function) at each of these points, solving the radiative transfer equation implies accumulating the effects of the emission and absorption coefficients at each point, from i = 0 (i.e. the points located farthest away from the observer, along the line of sight) to the last point in the system (i.e. the closest one to the observer). At each step, the emission coefficient contributes to increasing the accumulated intensity, while the absorption coefficient does the opposite (unless it is negative, in which case we have a maser).

This method provides a good approximation to the real case, provided that the number of points (i.e. the spatial resolution of the model nebula) is large enough. However, the accumulated sum of the effects of emission and subtraction can result in computational rounding errors. To test to which extent this occurred within SHAPE, we modelled a simple nebula for which an analytical solution exists under several assumptions.

The common ground for each of these models was a spherical nebula with radius r = 1.2 × 10¹⁶ cm located 100 pc away (projecting to a radius of 8 arcsec on the plane of the sky) and filled with gas with homogeneous temperature T = 100 K and density n = 1.5 × 10⁵ cm^-3, and micro-turbulence δ_V = 2.5 km s^-1. Under the assumption of homogeneity, it can be shown that in the optically thick case, the emitted intensity is given by the source function: ${\mathit{I}}_{\mathrm{thick}}\mathrm{(}\mathrm{W}\mathrm{H}{\mathrm{z}}^{-1}{\mathrm{m}}^{-2}\mathrm{s}{\mathrm{r}}^{-1}\mathrm{)}\mathrm{=}\frac{\mathrm{2}\mathit{h}{\mathit{\nu }}^{\mathrm{3}}}{{\mathit{c}}^{\mathrm{2}}}\frac{\mathrm{1}}{{\mathrm{e}}^{\frac{\mathit{h\nu }}{{\mathit{K}}_{\mathit{B}}{\mathit{T}}_{\mathrm{ex}}}}\mathrm{-}\mathrm{1}}\mathit{,}$(3)where I is the emitted intensity per unit of time, frequency, surface, and solid angle, T_ex the excitation temperature, h is the Planck constant and k_B is the Boltzmann constant. In the general case, the intensity is decreased from that in the thick case by a factor that depends on the value of τ as $\mathit{I}\mathrm{=}{\mathit{I}}_{\mathrm{thick}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{\tau }}\mathrm{)}\mathit{.}$(4)For the ¹²CO J = 6 − 5 transition and the given parameters, I_thick = 1.24 × 10^-14 W Hz^-1 m^-2 sr^-1.

In all our tests, we used a grid consisting of 64³ image pixels that covered a cube with a side length 3 × 10¹⁶ cm inside of which the spherical nebula was located. We always considered a microturbulence value of δ_V = 2.5 km s^-1 and performed the calculations in SHAPE using 200 spectral bands for adequately sampling the given spectral profile.

4.1.1. Static case

We first considered a nebula without expansion and in thermal equilibrium, so the level population is adequately described within the thermalised-case hypothesis and the absorption and emission coefficients k_ν and j_ν are trivial to compute. It must be stressed that we did not make use of the LVG approximation in these tests (i.e. the coefficients were not computed by shapemol), since we are hereby checking the accuracy of radiative transfer itself within SHAPE. In this case, the optical depth is simply τ = k_νL, where L is the total length of the nebula along the line considered (i.e. in a spherical nebula, L = 2 × r if we look at the centre, and progressively smaller as we consider points nearer the edge).

To test the model, we integrated the emission over a square box with side 3 × 10¹⁶ cm (6 grid pixels, or 1.875 arcsec at the adopted distance) centred on the nebula and compared the values with the theoretical prediction mentioned above for several values of the optical depth (achieved by varying the abundance X over a wide range). The intensity of our SHAPE+shapemol model almost perfectly matched the analytical value at large τ (i.e. the optically thick case), 9.8 × 10¹³ W Hz^-1 sr^-1 according to Eq. (3). Errors slowly increase (in absolute value) with decreasing optical depth, as shown in Fig. 4 (black thick line), where the error is defined as Error (%) = $\mathrm{100}\mathrm{\times }\frac{{\mathit{I}}_{\mathrm{S}}\mathrm{-}{\mathit{I}}_{\mathrm{P}}}{\frac{{\mathit{I}}_{\mathrm{S}}\mathrm{+}{\mathit{I}}_{\mathrm{P}}}{\mathrm{2}}}$, where I_S is the intensity as computed via radiative transfer by SHAPE, and I_P the predicted, analytically obtained intensity.

At large τ, the radiative transfer very well matches the analytical prediction, with the error limited to −0.5% (negative meaning that the analytically calculated intensity is higher than that obtained via radiative transfer). This systematic discrepancy is to be expected, given the difference in the method used for computing the intensities: while the analytical prediction is based on the source function, the absorption and emission coefficients we have used in the radiative transfer come from population level computations based on the theoretical coefficients under the assumption of the thermalised case.

The behaviour of the error changes with decreasing τ, growing in absolute value until it stabilises at around −2% with small τ. Part of this uncertainty can be explained by the fact that the τ value used in the analytical prediction was the one corresponding to the centre of the box, while the curvature of the nebula makes the product kL slightly smaller as one considers points offset from the centre. This slight overestimate of the analytical τ results in a higher opacity and therefore a higher intensity. In this case, with the square box centred on the nebula, this phenomenon can explain approximately one third of the total uncertainty. We can consider the rest as a combination of software rounding errors and the limited spatial resolution (i.e. discretization).

Fig. 4 Error behaviour in SHAPE radiative transfer as a function of τ for several cases discussed in the text. The nebula is always a filled sphere with radius 8 arcsec. In the cases with offset, the first number corresponds to a pointing offset from the centre of the spherical nebula on the plane of the sky, while the second number refers to the offset towards the observer (i.e. the intensity is not measured at the profile peak at zero velocity, but at a certain point in the line wings, at a velocity corresponding to the emission in the examined region).

4.1.2. Constant velocity case

We next tested the case of a filled, spherical nebula expanding at a constant velocity of V = 25 km s^-1. With a ratio $\frac{\mathit{V}}{{\mathit{\delta }}_{\mathit{V}}}\mathrm{=}\mathrm{10}$ we are well within the assumption of the LVG approximation and can express the optical depth τ at a given position in the nebula as ${\mathit{\tau }}_{\mathrm{LVG}}\mathrm{=}\mathit{k}\sqrt{\mathit{\pi }}{\mathit{\delta }}_{\mathit{\nu }}\frac{\mathit{r}\mathit{/}\mathrm{\left(}\mathit{\nu }\frac{{\mathit{V}}_{\mathrm{r}}}{\mathit{c}}\mathrm{\right)}}{\mathrm{1}\mathrm{+}{\mathit{\mu }}^{\mathrm{2}}\mathrm{(}\mathit{ϵ}\mathrm{-}\mathrm{1}\mathrm{)}}\mathit{,}$(5)where k is the normalized absorption coefficient, V_r is the radial component of the velocity, r is the linear distance from the centre of the nebula to the point considered, and μ = cosα, with α the angle between the velocity vector and the line of sight.

In this case, we integrated the emission over a box the same size as in the static case, but centred 3 × 10¹⁵ cm away from the centre of the spherical nebula (2 arcsec on the plane of the sky), and considered the emitted intensity not at the centre of the spectral profile (i.e. corresponding to V_r = 0 km s^-1) like in the previous case, but at a certain point in the profile wings. In particular, we examined the region located 6 × 10¹⁵ cm towards the observer, which, together with the offset of the box on the plane of the sky, corresponds to the emission at −22.4 km s^-1 of the integrated spectral profile.

The error (see Fig. 4, red dotted line) is the same as in the static case at large τ (i.e. optically thick nebula), and is due to the same reason. The absolute value of the error slowly grows as τ decreases, reaching values between −1% and −2.3%. The reasons for this uncertainty are the same as in the previous, static case.

4.1.3. Linear velocity case

Most PNe and PPNe have an approximately ballistic expansion (i.e. under the effect of no forces) after a brief period of shaping; they are said to follow a linear velocity law, with the velocity of each gas particle proportional to the distance of the particle to the star. To reflect this behaviour, we tested the case of a filled, spherical nebula similar to the previous case except for its velocity, with V ∝ r. The velocity was 0 at the location of the central star, and 25 km s^-1 at the edge of the nebula. The LVG approximation can be used here, with the advantage that the equation for τ_LVG becomes simpler since ϵ = 1 and the dependence on μ is removed.

As in the previous case, we considered the emitted intensity at a location away from the centre of the nebula both on the plane of the sky and towards the observer. We examined two such locations, the former centred 6 × 10¹⁵ cm away from the centre of the nebula (4 arcsec on the plane of the sky) and 7.5 × 10¹⁵ cm towards the observer (corresponding to the emission at −15.6 km s^-1); the latter aligned with the centre of the nebula but offset 9.6 × 10¹⁵ cm towards the observer (corresponding emission at −20 km s^-1).

The errors behave, in both cases, in a similar way to that of the static case: The error is fixed at −0.5% at large τ, and moderately grows with decreasing τ until reaching stable values around −2% in the optically thin regime (see Fig. 4, blue dotted-dashed and green dashed lines). These errors are attributable to the same factors as in the previously discussed static case.

4.2. Linear interpolation

shapemol relies on pre-generated tables of the absorption and emission coefficients k_ν and j_ν for a broad range of abundances, $\frac{\mathit{r}}{\mathit{V}}$ ratios, temperatures, and densities. Once a table has been selected based on the product $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$, the given coefficient is linearly interpolated from the adjacent tabulated values in temperature and density. Given the smooth but non-linear distribution of the values and the discretization applied (199 values of the temperature between 10 and 1000 K, 21 values of the density between 10² and 10⁷ cm^-3), this introduces an additional source of uncertainty.

To test the relative importance of this source of uncertainty, we considered the difference between the linearly interpolated value and the “real” value at a given point as the second-order term of the Taylor-series development of the desired coefficient k_ν or j_ν as a function of either the density or the temperature, with the other one remaining fixed. We estimated this term by computing the second derivative of the function in the neighbouring points of the table for every point. The associated error was given as a percentage of the value of the k or j coefficient at that point. However, since such a procedure would lead to divisions by almost zero wherever the coefficient is close to zero, and thus to unrealistic errors, we only took into account values within the top 80% of the highest value of the coefficient k_ν or j_ν for the given function.

The resulting errors for the k_ν and j_ν coefficients are below 1% for the vast majority of the tables, with some peculiarities worth noting: the error is slightly larger at low-J transitions (i.e. J = 1 − 0 to 4−3) and very low temperatures, in particular 15 K. In that case, the error reaches 2.8% for both k_ν and j_ν, except in the 1−0 transition at 15 K, where the k coefficient shows a 7% error.

5. Limitations

This section briefly describes the limitations of shapemol, which the user should be aware of.

5.1. Very faint emission

shapemol is not meant for studying very faint emission from extremely sub-excited, high-J transitions. For instance, the theoretical emission coefficient j_ν for the ¹²CO J = 16 − 15 transition at a temperature of 15 K (at which these levels have a very low population) is ~10¹⁰ larger than at 10 K. In shapemol, we have set every value of j_ν and k_ν below a certain very restrictive threshold to 0 to avoid computational noise. In particular, the coefficients are set to 0 wherever the energy (in K units) of the levels involved is ≳28 times higher than the corresponding equivalent rotational temperature describing the level population. In practice, however, cutting off the emission and absorption from high-J transitions at very low temperatures is unlikely to become a problem, since observational noise is usually much higher than the emission one would get from such a sub-excited transition, and no such emission is detected in actual cases. In most cases, moreover, the contribution from other nebular components that are not this sub-excited will completely dominate the total emission. The user should nevertheless be aware that for very faint emission (i.e. at a level below ~5% of the highest value of the coefficient j_ν or k_ν in a given table and transition), errors in the linear interpolation discussed in Sect. 4.2 will be significantly larger than stated there. In particular, errors in the interpolation will be large when crossing the jump from a zero value to the first non-zero value.

At these faint levels of emission, the contribution of the IR continuum from the central star or surrounding dust (i.e. the effect of vibrational transitions, see Appendix A), which is not accounted for by shapemol, could be significant.

5.2. Abuse of LVG conditions

The LVG approximation adequately describes the physical excitation conditions in most cases, even when the velocity is similar to the microturbulence value, δ_V (see details in Bujarrabal & Alcolea 2013). It can be applied to almost all PNe and PPNe, except in the most peculiar objects. However, in some cases, this approximation can give no more than a rough, qualitative idea of the excitation conditions of the target. The reason is that the emission and absorption coefficients, j_ν and k_ν, might be overestimated or underestimated when the LVG coherence zone (the region where particles might have similar velocities and thus are prone to photon absorption and interaction) is much larger or smaller than in the actual nebula.

Consider, for example, a shell expanding at a very low velocity (i.e. similar to δ_V). For any considered location in the shell, the corresponding LVG coherence zone will be a large, imaginary sphere centred on that location. The LVG approximation will consider photons coming from every point inside that sphere in the absorption and emission coefficient calculation (i.e. the effect usually called photon-trapping). However, if the shell is thin and the velocity is low enough, the LVG coherence zone sphere will be larger than the width of the shell, thus reaching farther towards empty space. shapemol, using the LVG approximation in the level population calculations, will then falsely consider absorption and emission from photons coming from an empty region. In other words, photon-trapping will be overestimated, and so k_ν and j_ν will be, specially for high-J transitions.

Another example worth considering is that of a disk in Keplerian rotation. Here, when V is very low or high (at large or small distances from the star respectively), the LVG coherence zone is too large or small compared to the actual photon-trapping region, and k_ν and j_ν will be incorrectly computed. This leads to a result that, although it can give an idea of the excitation conditions in the disk, cannot be considered quantitatively accurate.

Note, however, that this problem only affects the LVG calculations of j_ν and k_ν; instead, radiative transfer solving in SHAPE remains unaffected and is as precise as in any other case.

5.3. Masers

Maser emission will occur whenever the absorption k_ν coefficient is <0. This occurs for rare, specific values of n, T, and $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$, and mostly in very low-J CO transitions. Each point along the line of sight with such an absorption coefficient will amplify the arriving intensity. Real masers in CO are never very intense (a maser is considered intense when τ< − 1) under conditions of collisional excitation, since the pumping of CO masers by collisions is always inefficient. However, shapemol could spuriously yield high amplification if many of the model points share those specific values of the parameters and thus many points with negative k_ν are accumulated along the line of sight. In such a situation, which is in fact extremely improbable, a false maser several orders of magnitude more intense than the true emission will be conspicuous in the synthetic spectral profile. Note, however, that this problem can only occur when the LVG calculations of j_ν and k_ν are used for real conditions that would in fact not allow the use of the LVG approximation.

In the current version of SHAPE, a warning is displayed at the radiative transfer solving stage whenever τ< − 1 in a given line of sight and velocity.

6. Applying shapemol: the case of NGC 6302

Fig. 5 3D mesh view of the model of the molecular envelope of NGC 6302. The size of each frame is 40 × 40 arcsec2, and its label corresponds to the usual convention within SHAPE: front view, top view, and side view correspond to the view from Earth (north up, east to the left), from the direction defined by north in the plane of the sky, and from the direction defined by west in the plane of the sky (so the observer is on the left).

Fig. 6 Resulting synthetic spectra (red, dotted-dashed line) and observations (black histogram) for the 12CO and 13CO transitions detected in NGC 6302 with HIFI. Ifactor refers to the intensity factor applied to the model to account for the uncertainties in the radiative transfer solving and in calibrating the observations.

6.1. Motivation

NGC 6302 is a spectacular young PN. Located 1.17 ± 0.19 kpc away (Meaburn et al. 2008), it hosts a very hot central star with a T_eff of at least 150 000 K (Wright et al. 2011). The nebula shows an extreme bipolar morphology in optical images, along with an equatorial dust lane that most likely is associated with a molecular torus or disk (Matsuura et al. 2005). Authors have derived different values of the mass of this molecular structure, which, when scaled to the distance adopted in this work (with mass depending on the distance squared), can be used to establish valid comparisons. In particular, Peretto et al. (2007) analysed low-excitation CO interferometric observations and found an expanding torus with a mass of ~2.7 M_⊙. This is completely at odds with the masses derived by other authors: Gomez et al. (1989) found ~0.14 M_⊙, while Huggins & Healy (1989) found a slightly higher mass of 0.3 M_⊙. Dinh-V-Trung et al. (2008) computed a mass of ~0.1 M_⊙ through the means of ¹²CO and ¹³CO J = 2−1 interferometric mapping. More recently, Bujarrabal et al. (2012) presented high-J CO transitions observed with Herschel/HIFI and provided some constraints in density and temperature that were compatible with the findings by Dinh-V-Trung et al.

Careful consideration of the data leads to the conclusion that the molecular component of the nebula must be broken into clumps of dense material, with other regions completely devoid of molecules. In other words, a simple toroidal model with constant values of the density and temperature provides a rudimentary fit to the data, at most. In an effort to better constrain the mass and excitation conditions of the nebula, we have built a model of the molecular component of NGC 6302 and compared it with both the low-J interferometric map by Dinh-V-Trung et al. and the high-J observations presented by Bujarrabal et al. (2012).

The main motivation behind this model is to show the potential of shapemol. Because it is integrated within SHAPE and because of the 3D engine the model can be modified and adapted with a few mouse clicks, including its geometry, which can be as complex and detailed as needed. On the other hand, radiative transfer modelling with shapemol runs in a few tens of seconds on current standard computers, allowing quick iteration while still providing a reasonable approximation to the actual physical conditions.

Additionally, this model shows the complexity of the molecular envelope of NGC 6302 and allows us to determine the physical and chemical properties of its molecular envelope, as well as to provide a full description of the 3D morphology of the nebula assuming a homologous expansion for the velocity field. We investigate the characteristic densities and temperatures of the nebula and provide an independent, alternate determination of the molecular mass of the nebula.

6.2. The model

The model consists of a broken torus with the same orientation as the simple model by Dinh-V-Trung et al. plus a series of fingers and blobs expanding outwards from the central star (see the sketch in Fig. 5). This torus is broken into two main regions, named blue and red according to their Doppler shift with respect to the central star. Each of these regions is a piece of a cylindrical ring, each different in size and span, and is itself split into two components: an inner, thinner, hotter, and much denser component, and an outer, larger, cooler, and more tenuous one. The model also contains three spherical blobs with different physical conditions, and two outward-projected cylindrical filaments we have called fingers, each one split into an inner and an outer section.

All these components are in fact directly identifiable in the interferometric maps, which show emission from a structure resembling a broken torus with smaller condensations at different positions and velocities (see Fig. 1 in Dinh-V-Trung et al. 2008). The positions of several of these condensations show a dependence on their velocity. Careful consideration of the high-J profiles from HIFI allows us to associate some of these structures with their profile counterparts, revealing that the torus must have an inner region with a higher temperature and density.

To keep the model as simple as possible while providing a good fit to the data, we have fixed the temperature and density of each structure to a single, representative value. All the structures are expanding away from the star with velocity laws that are either constant values (in structures that are too thin, in the radial direction, to produce an observable difference in the profile) or ballistic velocity patterns (in radially thicker structures). Each structure has its own micro-turbulence velocity value, δ_V. On the other hand, we have set the ¹²CO and ¹³CO abundances to two unique values for the whole nebula.

Figure 5 shows a three-dimensional view of the model with all structures visible. We fitted the model in the same iterative fashion as in the modelling of NGC 7027 (Santander-García et al. 2012). The resulting best-fit of the Herschel/HIFI data is shown in Fig. 6, while the interferometric ¹²CO and ¹³CO J = 2 − 1 synthetic maps are shown in Figs. 7 and 8. Because of the uncertainty in the HIFI data fluxes, which are of the order of 20−25% (see Santander-García et al. 2012), we applied different free intensity scale factors to each of the model spectrum transitions (i.e. the factor to be applied to the intensity of the modelled profile to account for computational errors, and to be considered as the conservative error of the intensity of our model). In the best-fit model, the deviation of these factors from unity is within 30% (i.e. the factors range from 0.77 to 1.3). Similar uncertainties are expected in the interferometric observations by Dinh-V-Trung et al. (2008), with the additional problem that about 20% of the flux is lost in those observations in a hardly predictable manner. Our model fits these maps fairly well, with differences between model calculations and data within observational uncertainties, although some strong observational variations from channel to channel could not be reproduced by our model.

The systemic velocity of the nebula is V_LSR = −31 km s^-1. The relative abundance of the best-fit is 1.5 × 10^-4 for ¹²CO and 2.5 × 10^-5 for ¹³CO. The derived ¹²CO/¹³CO abundance ratio may appear to be low, but we note that it is not unusual among post-AGB nebulae, particularly around O-rich stars (see the thorough discussion in Bujarrabal et al. 2013).

The results are summarised in Tables 1–3, along with a conservative estimate of the uncertainties. These were computed by varying a given parameter while leaving the others unchanged, so that the deviation from the data remained within the aforementioned 30%. For the velocities, the uncertainty was derived by eye, by varying the velocity until the model did no longer provide a good fit. Note that, given the lack of information on the geometry of the nebula (because of the poor spatial resolution of the observations) and the main motivation of the model, we did not consider uncertainties on the geometry parameters.

6.3. Discussion

Channel-to-channel variations in the interferometric data of NGC 6302 show that the real nebula must be full of clumps, matter condensations and voids. Our model, however, provides an approximate description of the general properties and dynamics of the molecular envelope of this nebula and its characteristic physical conditions.

The densities and temperatures of the different structures of our model are in general compatible with previous results by Dinh-V-Trung et al. and Bujarrabal et al. In particular, the temperature of the fastest blue-shifted structures (referred to as fast winds by Bujarrabal et al.) is slightly higher (50−80 K compared to 40 K), and the densities somewhat lower (~10⁴ cm^-3 compared to ~10⁵ cm^-3 ) than found by these authors. The density of the inner regions of the torus is of the order of a few 10⁵ cm^-3, however, in good agreement with the findings by Bujarrabal et al. (2012).

No signs of shock-accelerated structures are apparent in our modelling, other than the fingers, which lie outside the plane of the torus and move away from the star, reaching velocities of up to 28 and 38 km s^-1. Note that although these structures might be at a different position along the same line of sight, the velocities are derived from the interferometric maps and thus are well constrained. The fingers seem to show slightly higher velocities (1−2 km s^-1) in their inner, hotter regions than in the neighbouring boundaries with the cooler, outer part of the same fingers. This might be indicative of a slow shock developing through the fingers, although there is not enough evidence for this assertion. The fingers share the main properties of the rest of the nebula at optical wavelengths: their projected expansion velocity is similar to the observed expansion in the same region at optical wavelengths (Lago & Costa 2014), while their densities are similar to those found in the ionized region by Rauber et al. (2014). To some extent, the fingers resemble the cometary tails found in the Helix (O’dell & Handron 1996), with a filamentary structure and a shielded region undergoing erosion.

Modelling low-J and high-J transitions together has allowed us to better characterise the physical conditions of the molecular envelope of NGC 6302. We derive a mass of 0.08M_⊙ for the torus and 0.03 M_⊙ for the blobs and fingers, giving a total mass for the CO-emitting region of 0.11 M_⊙. This result is completely at odds with the mass computed by Peretto et al. (2007), but in full agreement with the ~0.1 M_⊙ found by Dinh-V-Trung et al. (2008). Note that these authors used a simplified method to estimate the total mass from the total emission of a low-J¹³CO line, which is expected to be optically thin and only slightly dependent on the excitation state.

Fig. 7 Resulting synthetic interferometric map of NGC 6302 in the 12CO J = 2 − 1 transition. The size of each frame is 20 × 20 arcsec2. The LSR velocity is indicated in the upper left corner of each frame. The size and orientation of the synthetic beam, and the contour levels have been adapted for approximate visual comparison with the observations displayed in Fig. 1 in Dinh-V-Trung et al. (2008); emission is represented in a logarithmic scale, with the following maxima (10th contour) depending on the label: 2.4 Jy beam-1 for “c1” channels, 6.4 Jy beam-1 for “c2” channels, and 11 Jy beam-1 for “c3” channels.

Fig. 8 Resulting synthetic interferometric map of NGC 6302 in the 12CO J = 2 − 1 transition. The size of each frame is 20 × 20 arcsec2. The LSR velocity is indicated in the upper left corner of each frame. The size and orientation of the synthetic beam, and the contour levels have been adapted for approximate visual comparison with the observations displayed in Fig. 2 in Dinh-V-Trung et al. (2008); emission is represented in a logarithmic scale, with the following maxima (10th contour) depending on the label: 2.7 Jy beam-1 for “c1” channels, 5.4 Jy beam-1 for “c2” channels, and 7.5 Jy beam-1 for “c3” channels.

6.4. Summary

We have built a model of the molecular envelope of NGC 6302 to investigate its physical and excitation conditions. The large-scale 3D structure of our model molecular envelope was determined mostly from the direct interpretation of the interferometric maps and high-J line profiles. Provided a simple velocity pattern, as typically observed in PNe, this structure is the simplest one that can account for all the observations at hand.

The nebula was modelled as a broken ring with an inner, hotter region and an extended, colder, outer one, together with a series of blobs and fingers expanding outwards from the plane of the ring. The velocity follows a ballistic expansion pattern except in the thinner structures, which have single, constant values of the velocity. The densities are in the range of 10⁴ cm^-3 and 10⁵ cm^-3 for the fingers and ring, while the temperatures range from a few tens of K in the outer regions of the ring to 300 K in the inner, exposed region. The mass of the molecular envelope is 0.11M_⊙.

Our model provides a fair fit to the data. It is clear, however, that the intrinsic, small-scale structure of the nebula must be full of clumps and voids, as sudden channel-to-channel variations, which our model fails to reproduce, seem to imply in the interferometric maps.

Finally, we stress that complex 3D structures and dynamics can be modelled with a reasonable effort because shapemol and its implementation within SHAPE is fast and flexible.

Acknowledgments

This work was partially supported by grant “UNAM-PAPIIT 100410 and 101014” and by the Spanish MICINN within the program CONSOLIDER INGENIO 2010, under grant “Molecular Astrophysics: The Herschel and ALMA Era, ASTROMOL” (ref.: CSD2009-00038).

References

Appendix A: Effects of vibrational transitions on the emission of considered rotational transitions

In our general calculations, we only considered rotational levels within the ground-vibrational state, assuming that vibrational transitions have little efect on the emission of the relevant v = 0 rotational lines. To validate this assumption, we performed calculations in some relevant cases, also including the v = 1 rotational levels (40 additional levels) and an IR continuum source. The IR source was always assumed to be central and relatively small compared to the distance at which the studied point is placed in the nebula (both assumptions are obviously not relevant to our goals). The IR spectrum is assumed to be that of a black body for a given temperature of about 500−10 000 K,

whose total size is calculated to yield total IR intensities at 4.7 microns equal to the flux actually observed in the studied nebulae (taking into account the distance to them). This IR source can be a compact hot-dust region or a stellar surface in cool stars. We verified that the assumption of a black-body radiation spectrum and its temperature have negligible effects on the final results, the relevant parameters being (as expected) the IR source intensity and the distance between this source and the considered element of molecule-rich gas. We also verified that in all the studied cases the collisional vibrational excitation has negligible effects on the excitation of the v = 0 levels.

We considered three different model nebulae (see definitions of the parameters in Sect. 3):

• A)

  A case similar to conditions expected in the Red Rectangle (apost-AGB source in which most molecule-rich gas is placed in arotating disk, Bujarrabal & Alcolea 2013):distance to the object,D = 710 pc, F(4.7μ) = 100 Jy, V = 1.7 km s^-1, r = 10¹⁶ cm, relevant densities 10⁵ − 10⁶ cm^-3. Since in the Red Rectangle rotation dominates dynamics, we have assumed a low value of ϵ, 0.1 to simulate absorption by a long path between the centre and the considered point. We recall that the LVG approach is not accurate in the case of rotation, but the results of our tests give an idea of the effects of the IR continuum in this interesting example.

• B)

  A case similar to conditions expected in the young PN NGC 7027 (Santander-García et al. 2012): D = 1000 pc, F(4.7μ) = 50 Jy, V = 15 km s^-1, r = 10¹⁷ cm, relevant densities ~10⁴–10⁵ cm^-3. Since in this source the

  Fig. A.1 Model calculations of TR as defined in Appendix A for the conditions corresponding to case A) in the text and an ϵ value of 0.1.

  Fig. A.2 Model calculations of TR as defined in Appendix A for the conditions corresponding to case B) in the text, for an ϵ value of 0.1.

  Fig. A.3 Close-up of case B) for an ϵ value of 0.1.

  Fig. A.4 Model calculations of TR as defined in Appendix A for the conditions corresponding to case B) in the text, for an ϵ value of 1.

  Fig. A.5 Model calculations of TR as defined in Appendix A for the conditions corresponding to case C) in the text, for an ϵ value of 0.1.

  Fig. A.6 Model calculations of TR as defined in Appendix A for the conditions corresponding to case C) in the text, for an ϵ value of 1.

  Fig. A.7 Close-up of case C) for an ϵ value of 1. Deviations from the model are only significant for densities much lower than those found in the real nebula (see references in text).

  velocity gradients are often large, we took ϵ = 1, but we also present calculations for a very low velocity gradient ϵ = 0.1 to see possible effects due to this change.

• C)

  A case similar to conditions expected in the PPNe CRL 618 (Soria-Ruiz et al. 2013, Sánchez Contreras et al. 2004): D = 1000 pc, F(4.7 μ) = 15 Jy, V = 15 km s^-1, r = 10¹⁶ cm, relevant densities ~ 5 × 10⁶ cm^-3. Some components of this source (e.g. the very fast outflow) show a strong velocity gradient, but in others (the fossil AGB envelope) the velocity increases slowly, so in this case we also considered ϵ = 1, 0.1.

We mostly calculated the intensity for a J = 16−15 transition, which has been widely observed with Herschel/HIFI, since low-J transitions are more easily thermalised and their emissivities are more easy to describe. We present in Figs. A.1 to A.5 predictions for the A), B), C) cases of the Raighleigh-Jeans equivalent brightness temperature of ¹²CO ¹³CO J = 2 − 1 transition:

${\mathit{T}}_{\mathrm{R}}\mathrm{=}\frac{\mathit{k}}{\mathit{h\nu }}\mathrm{[}\mathit{I}\mathrm{-}\mathit{I}\mathrm{(}\mathrm{BG}\mathrm{\left)}\mathrm{\right]}\mathit{,}$where I is the intensity at the line centre emitted in the direction perpendicular to the radial one and I(BG) is the intensity of the cosmic background at the line frequency.

All Tables

All Figures

Fig. 1 Examples of the jν coefficient corresponding to the 12CO transitions with J = 2–1, 6–5, 10–9 and 16–15 as a function of temperature in six different physical conditions, with varying ϵ, n and $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$. Some of the transitions have been scaled up or down by the specified factors for displaying purposes. The non-linear dependance on the density n can easily be seen here, as well as the almost linear dependance on the product $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$. The value of ϵ, on the other hand, affects the jν coefficient only slightly.

In the text

Fig. 2 Examples of the kν coefficient corresponding to the 12CO transitions with J = 2–1, 6–5, 10–9 and 16–15 as a function of temperature in six different physical conditions, with varying ϵ, n and $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$. Some of the transitions have been scaled up or down by the specified factors for displaying purposes. The non-linear dependance on the density n can easily be seen here, as well as the almost linear dependance on the product $\frac{\mathit{r}}{\mathit{V}}\mathit{X}$. The value of ϵ, on the other hand, affects the kν coefficient only slightly.

In the text

	Fig. 3 Workflow of shapemol.
In the text

Fig. 4 Error behaviour in SHAPE radiative transfer as a function of τ for several cases discussed in the text. The nebula is always a filled sphere with radius 8 arcsec. In the cases with offset, the first number corresponds to a pointing offset from the centre of the spherical nebula on the plane of the sky, while the second number refers to the offset towards the observer (i.e. the intensity is not measured at the profile peak at zero velocity, but at a certain point in the line wings, at a velocity corresponding to the emission in the examined region).

In the text

Fig. 5 3D mesh view of the model of the molecular envelope of NGC 6302. The size of each frame is 40 × 40 arcsec2, and its label corresponds to the usual convention within SHAPE: front view, top view, and side view correspond to the view from Earth (north up, east to the left), from the direction defined by north in the plane of the sky, and from the direction defined by west in the plane of the sky (so the observer is on the left).

In the text

	Fig. 6 Resulting synthetic spectra (red, dotted-dashed line) and observations (black histogram) for the 12CO and 13CO transitions detected in NGC 6302 with HIFI. Ifactor refers to the intensity factor applied to the model to account for the uncertainties in the radiative transfer solving and in calibrating the observations.
In the text

Fig. 7 Resulting synthetic interferometric map of NGC 6302 in the 12CO J = 2 − 1 transition. The size of each frame is 20 × 20 arcsec2. The LSR velocity is indicated in the upper left corner of each frame. The size and orientation of the synthetic beam, and the contour levels have been adapted for approximate visual comparison with the observations displayed in Fig. 1 in Dinh-V-Trung et al. (2008); emission is represented in a logarithmic scale, with the following maxima (10th contour) depending on the label: 2.4 Jy beam-1 for “c1” channels, 6.4 Jy beam-1 for “c2” channels, and 11 Jy beam-1 for “c3” channels.

In the text

Fig. 8 Resulting synthetic interferometric map of NGC 6302 in the 12CO J = 2 − 1 transition. The size of each frame is 20 × 20 arcsec2. The LSR velocity is indicated in the upper left corner of each frame. The size and orientation of the synthetic beam, and the contour levels have been adapted for approximate visual comparison with the observations displayed in Fig. 2 in Dinh-V-Trung et al. (2008); emission is represented in a logarithmic scale, with the following maxima (10th contour) depending on the label: 2.7 Jy beam-1 for “c1” channels, 5.4 Jy beam-1 for “c2” channels, and 7.5 Jy beam-1 for “c3” channels.

In the text

	Fig. A.1 Model calculations of TR as defined in Appendix A for the conditions corresponding to case A) in the text and an ϵ value of 0.1.
In the text

	Fig. A.2 Model calculations of TR as defined in Appendix A for the conditions corresponding to case B) in the text, for an ϵ value of 0.1.
In the text

	Fig. A.3 Close-up of case B) for an ϵ value of 0.1.
In the text

	Fig. A.4 Model calculations of TR as defined in Appendix A for the conditions corresponding to case B) in the text, for an ϵ value of 1.
In the text

	Fig. A.5 Model calculations of TR as defined in Appendix A for the conditions corresponding to case C) in the text, for an ϵ value of 0.1.
In the text

	Fig. A.6 Model calculations of TR as defined in Appendix A for the conditions corresponding to case C) in the text, for an ϵ value of 1.
In the text

	Fig. A.7 Close-up of case C) for an ϵ value of 1. Deviations from the model are only significant for densities much lower than those found in the real nebula (see references in text).
In the text