Issue 
A&A
Volume 573, January 2015



Article Number  A56  
Number of page(s)  18  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201322348  
Published online  17 December 2014 
SHAPEMOL: a 3D code for calculating CO line emission in planetary and protoplanetary nebulae
Detailed modelfitting of the complex nebula NGC 6302^{⋆}
^{1} Observatorio Astronómico Nacional, Ap. de Correos 112, 28803 Alcalá de Henares, Madrid, Spain
email: m.santander@oan.es
^{2} Instituto de Ciencia de Materiales de Madrid (CSIC), 28049 Madrid, Spain
^{3} Department of Physics & Astronomy, University of Calgary, Calgary, Alberta, Canada
^{4} Instituto de Astronomía Universidad Nacional Autónoma de México, C.P. 22860, Ensenada, Mexico
Received: 23 July 2013
Accepted: 21 October 2014
Context. Modern instrumentation in radioastronomy constitutes a valuable tool for studying the Universe: ALMA has reached unprecedented sensitivities and spatial resolution, while Herschel/HIFI has opened a new window (most of the submm and farinfrared ranges are only accessible from space) for probing molecular warm gas (~50−1000 K). On the other hand, the software SHAPE has emerged in the past few years as a standard tool for determining the morphology and velocity field of different kinds of gaseous emission nebulae via spatiokinematical modelling. Standard SHAPE implements radiative transfer solving, but it is only available for atomic species and not for molecules.
Aims. Being aware of the growing importance of the development of tools for simplifying the analyses of molecular data from newera observatories, we introduce the computer code shapemol, a complement to SHAPE, with which we intend to fill the sofar underdeveloped molecular niche.
Methods. shapemol enables userfriendly, spatiokinematic modelling with accurate nonLTE calculations of excitation and radiative transfer in CO lines. Currently, it allows radiative transfer solving in the ^{12}CO and ^{13}CO J = 1−0 to J = 17−16 lines, but its implementation permits easily extending the code to different transitions and other molecular species, either by the code developers or by the user. Used along SHAPE, shapemol allows easily generating synthetic maps to test against interferometric observations, as well as synthetic line profiles to match singledish observations.
Results. We give a full description of how shapemol works, and we discuss its limitations and the sources of uncertainty to be expected in the final synthetic profiles or maps. As an example of the power and versatility of shapemol, we build a model of the molecular envelope of the planetary nebula NGC 6302 and compare it with ^{12}CO and ^{13}CO J = 2−1 interferometric maps from SMA and highJ transitions from Herschel/HIFI. We find the molecular envelope to have a complex, broken ringlike structure with an inner, hotter region and several “fingers” and highvelocity blobs, emerging outwards from the plane of the ring. We derive a mass of 0.11 M_{⊙} for the molecular envelope.
Key words: molecular data / radiative transfer / ISM: kinematics and dynamics / planetary nebulae: general
A copy of the code is only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/573/A56
© 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 highresolution spectra in the range of 480−1250 GHz and 1410−1910 GHz, which includes transitions as high as ^{12}CO and ^{13}CO J = 16−15, unobservable from the ground – as is most of the submm and farinfrared 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 newera observatories, we introduce the computer code shapemol^{1}. This complementary software enhances the ability of SHAPE, the oftenused tool for performing spatiokinematical 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 moleculerich gas in these objects. Together, SHAPE+shapemol can generate highresolution 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 (SantanderGarcí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 wellknown 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, highenergy 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, positionvelocity 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 nonthermalised 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 ^{12}CO and ^{13}CO lines. This is done by interpolating the absorption and emission coefficients from a set of pregenerated 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^{4} grid cells would take several minutes to compute the exact solution in current regular computers, SHAPE+shapemol will handle models consisting of 10^{6} gridcells in a matter of 30−40 s. This, together with its integration within the userfriendly 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 ^{12}CO and ^{13}CO lines to be compared with submm and mmrange 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^{2}. 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: (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 mmwave transitions from diffuse media (interstellar or circumstellar clouds), φ_{ν} is always given by the local velocity dispersion, that is, the thermal or microturbulent velocities. The (nonLTE) level populations used to solve the transfer equation are calculated using the wellknown LVG approximation.
The level populations depend on the collisional transition rates and the radiative excitation and deexcitation 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 wellknown 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 postAGB shells; see detailed comparison with nonlocal 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 nonlinear way on the density n and the temperature T of each point, and almost linearly on the product , 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, , and ϵ.
In our calculations, we used the most recent collisional rates from the LAMDA database^{3}, Schöier et al. (2005), and rates calculated by Yang et al. (2010) for collisions with ortho and paraH_{2}, 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_{2}.
Absorption and emission coefficients were obtained for the 17 lowest rotational transitions of the groundvibrational state, practically the only ones observed. In calculating them, a total of 40 rotational levels were considered, although in lowexcitation 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 groundvibrational 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 generalpurpose design of shapemol, as well as its userfriendliness 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 rotationalline 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 highJ 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 ^{12}CO 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 . Some of the transitions have been scaled up or down by the specified factors for displaying purposes. The nonlinear dependance on the density n can easily be seen here, as well as the almost linear dependance on the product . The value of ϵ, on the other hand, affects the j_{ν} coefficient only slightly. 
Fig. 2
Examples of the k_{ν} coefficient corresponding to the ^{12}CO 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 . Some of the transitions have been scaled up or down by the specified factors for displaying purposes. The nonlinear dependance on the density n can easily be seen here, as well as the almost linear dependance on the product . 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 timeconsuming 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 pregenerated tables of k_{ν} and j_{ν} as functions of n and T, each table corresponding to a species (^{12}CO and ^{13}CO in its present state), a value of ϵ, and a value of the 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 for each sampled position. Then, it selects the table with the closest value of ϵ and 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 , the software scales the computed absorption and emission coefficients according to the ratio of the sampled position value of 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 to that of the selected table. If, however, the option “Linear Interpolation” is set, shapemol selects the two tables whose 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 , and uses those to linearly interpolate the final k_{ν} and j_{ν} at the cell’s . This procedure is slightly more accurate, since no assumptions are made on the dependence of k_{ν} and j_{ν} on , 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 ^{12}CO and ^{13}CO. There are 199 values of T and 21 values of n in the available pregenerated tables, ranging from 10 to 1000 K in steps of 5 K for the temperature, and from 10^{2} cm^{3} to 10^{7} cm^{3} in multiplicative factors of 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 (SantanderGarcía et al. 2012).
3.3. Other molecular species
Future plans include the addition of more molecules, such as C^{18}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 usersupplied 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 – ^{12}CO or ^{13}CO – and transition for each structure of the nebula, and assign values for the molecular abundance, ϵ, and microturbulence, δ_{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 approximationbased 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 datacube consisting of a stack of images of the nebula, one per frequency resolution element (i.e. spectral channel). This datacube 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 HalfPower 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: (2)where s_{0} is the position of the edge of the target located farthest away from the observer, s is the position measured from s_{0} 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: . 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^{16} 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^{5} cm^{3}, and microturbulence δ_{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: (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 (4)For the ^{12}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^{3} image pixels that covered a cube with a side length 3 × 10^{16} 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 thermalisedcase 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^{16} 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^{13} 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 (%) = , 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 we are well within the assumption of the LVG approximation and can express the optical depth τ at a given position in the nebula as (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^{15} 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^{15} 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^{15} cm away from the centre of the nebula (4 arcsec on the plane of the sky) and 7.5 × 10^{15} 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^{15} 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 dotteddashed 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 pregenerated tables of the absorption and emission coefficients k_{ν} and j_{ν} for a broad range of abundances, ratios, temperatures, and densities. Once a table has been selected based on the product , the given coefficient is linearly interpolated from the adjacent tabulated values in temperature and density. Given the smooth but nonlinear 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^{2} and 10^{7} 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 secondorder term of the Taylorseries 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 lowJ 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 subexcited, highJ transitions. For instance, the theoretical emission coefficient j_{ν} for the ^{12}CO J = 16 − 15 transition at a temperature of 15 K (at which these levels have a very low population) is ~10^{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 highJ 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 subexcited transition, and no such emission is detected in actual cases. In most cases, moreover, the contribution from other nebular components that are not this subexcited 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 nonzero 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 photontrapping). 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, photontrapping will be overestimated, and so k_{ν} and j_{ν} will be, specially for highJ 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 photontrapping 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 , and mostly in very lowJ 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 arcsec^{2}, 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, dotteddashed line) and observations (black histogram) for the ^{12}CO and ^{13}CO transitions detected in NGC 6302 with HIFI. I_{factor} 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 lowexcitation 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_{⊙}. DinhVTrung et al. (2008) computed a mass of ~0.1 M_{⊙} through the means of ^{12}CO and ^{13}CO J = 2−1 interferometric mapping. More recently, Bujarrabal et al. (2012) presented highJ CO transitions observed with Herschel/HIFI and provided some constraints in density and temperature that were compatible with the findings by DinhVTrung 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 lowJ interferometric map by DinhVTrung et al. and the highJ 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 DinhVTrung 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 outwardprojected 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 DinhVTrung et al. 2008). The positions of several of these condensations show a dependence on their velocity. Careful consideration of the highJ 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 microturbulence velocity value, δ_{V}. On the other hand, we have set the ^{12}CO and ^{13}CO abundances to two unique values for the whole nebula.
Figure 5 shows a threedimensional view of the model with all structures visible. We fitted the model in the same iterative fashion as in the modelling of NGC 7027 (SantanderGarcía et al. 2012). The resulting bestfit of the Herschel/HIFI data is shown in Fig. 6, while the interferometric ^{12}CO and ^{13}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 SantanderGarcí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 bestfit 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 DinhVTrung 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.
Bestfit model parameters for the molecular ring components of NGC 6302.
Bestfit model parameters for the blob components of NGC 6302.
Bestfit model parameters for the finger components of NGC 6302.
The systemic velocity of the nebula is V_{LSR} = −31 km s^{1}. The relative abundance of the bestfit is 1.5 × 10^{4} for ^{12}CO and 2.5 × 10^{5} for ^{13}CO. The derived ^{12}CO/^{13}CO abundance ratio may appear to be low, but we note that it is not unusual among postAGB nebulae, particularly around Orich 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
Channeltochannel 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 DinhVTrung et al. and Bujarrabal et al. In particular, the temperature of the fastest blueshifted 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^{4} cm^{3} compared to ~10^{5} cm^{3} ) than found by these authors. The density of the inner regions of the torus is of the order of a few 10^{5} cm^{3}, however, in good agreement with the findings by Bujarrabal et al. (2012).
No signs of shockaccelerated 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 lowJ and highJ 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 COemitting 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 DinhVTrung et al. (2008). Note that these authors used a simplified method to estimate the total mass from the total emission of a lowJ^{13}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 ^{12}CO J = 2 − 1 transition. The size of each frame is 20 × 20 arcsec^{2}. 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 DinhVTrung 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 ^{12}CO J = 2 − 1 transition. The size of each frame is 20 × 20 arcsec^{2}. 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 DinhVTrung 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 largescale 3D structure of our model molecular envelope was determined mostly from the direct interpretation of the interferometric maps and highJ 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^{4} cm^{3} and 10^{5} 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, smallscale structure of the nebula must be full of clumps and voids, as sudden channeltochannel 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 “UNAMPAPIIT 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.: CSD200900038).
References
 Bujarrabal, V., & Alcolea, J. 2013, A&A, 552, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bujarrabal, V., Alcolea, J., SoriaRuiz, R., et al. 2012, A&A, 537, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bujarrabal, V., Alcolea, J., Van Winckel, H., SantanderGarcía, M., & CastroCarrizo, A. 2013, A&A, 557, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Castor, J. I. 1970, MNRAS, 149, 111 [NASA ADS] [Google Scholar]
 DinhVTrung, Bujarrabal, V., CastroCarrizo, A., Lim, J., & Kwok, S. 2008, ApJ, 673, 934 [NASA ADS] [CrossRef] [Google Scholar]
 Gomez, Y., Rodriguez, L. F., Moran, J. M., & Garay, G. 1989, ApJ, 345, 862 [NASA ADS] [CrossRef] [Google Scholar]
 Huggins, P. J., & Healy, A. P. 1989, ApJ, 346, 201 [Google Scholar]
 Jones, D., Lloyd, M., SantanderGarcía, M., et al. 2010, MNRAS, 408, 2312 [NASA ADS] [CrossRef] [Google Scholar]
 Lago, P., & Costa, R. D. D. 2014, in Asymmetrical Planetary Nebulae VI conference, Proceedings of the conference held 4−8 November, 2013, eds. C. Morisset, G. DelgadoInglada, & S. TorresPeimbert, 52 [Google Scholar]
 Landi, E., Del Zanna, G., Young, P. R., Dere, K. P., & Mason, H. E. 2012, ApJ, 744, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Matsuura, M., Zijlstra, A. A., & Fuller, G. A. 2005, in Planetary Nebulae as Astronomical Tools, eds. R. Szczerba, G. Stasińska, & S. K. Gorny, AIP Conf. Ser., 804, 211 [Google Scholar]
 Meaburn, J., Lloyd, M., Vaytet, N. M. H., & López, J. A. 2008, MNRAS, 385, 269 [NASA ADS] [CrossRef] [Google Scholar]
 O’dell, C. R., & Handron, K. D. 1996, AJ, 111, 1630 [NASA ADS] [CrossRef] [Google Scholar]
 Peretto, N., Fuller, G., Zijlstra, A., & Patel, N. 2007, A&A, 473, 207 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ramstedt, S., Schöier, F. L., Olofsson, H., & Lundgren, A. A. 2008, A&A, 487, 645 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rauber, A. B., Copetti, M. V. F., & Krabbe, A. C. 2014, A&A, 563, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reader, J., Kramida, A., & Ralchenko, Y. 2012, in Am. Astron. Soc. Meet. Abstr., 219, 443.01 [Google Scholar]
 Sánchez Contreras, C., Bujarrabal, V., CastroCarrizo, A., Alcolea, J., & Sargent, A. 2004, ApJ, 617, 1142 [NASA ADS] [CrossRef] [Google Scholar]
 SantanderGarcía, M., Bujarrabal, V., & Alcolea, J. 2012, A&A, 545, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Smith, P. L., Esmond, J. R., Heise, C., & Kurucz, R. L. 1996, in UV and Xray Spectroscopy of Astrophysical and Laboratory Plasmas, eds. K. Yamashita, & T. Watanabe, 513 [Google Scholar]
 SoriaRuiz, R., Bujarrabal, V., & Alcolea, J. 2013, A&A, 559, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Steffen, W., GarcíaSegura, G., & Koning, N. 2009, ApJ, 691, 696 [NASA ADS] [CrossRef] [Google Scholar]
 Steffen, W., Koning, N., Wenger, S., Morisset, C., & Magnor, M. 2011, IEEE Transactions on Visualization and Computer Graphics, 17, 454 [Google Scholar]
 Teyssier, D., Hernandez, R., Bujarrabal, V., Yoshida, H., & Phillips, T. G. 2006, A&A, 450, 167 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Velázquez, P. F., Steffen, W., Raga, A. C., et al. 2011, ApJ, 734, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, N. J., Barlow, M. J., Ercolano, B., & Rauch, T. 2011, MNRAS, 418, 370 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, B., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2010, ApJ, 718, 1062 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Effects of vibrational transitions on the emission of considered rotational transitions
In our general calculations, we only considered rotational levels within the groundvibrational 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 hotdust region or a stellar surface in cool stars. We verified that the assumption of a blackbody 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 moleculerich 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 (apostAGB source in which most moleculerich 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^{16} cm, relevant densities 10^{5} − 10^{6} 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 (SantanderGarcía et al. 2012): D = 1000 pc, F(4.7μ) = 50 Jy, V = 15 km s^{1}, r = 10^{17} cm, relevant densities ~10^{4}–10^{5} cm^{3}. Since in this source the
Fig. A.1 Model calculations of T_{R} 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 T_{R} as defined in Appendix A for the conditions corresponding to case B) in the text, for an ϵ value of 0.1.
Fig. A.3 Closeup of case B) for an ϵ value of 0.1.
Fig. A.4 Model calculations of T_{R} 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 T_{R} 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 T_{R} as defined in Appendix A for the conditions corresponding to case C) in the text, for an ϵ value of 1.
Fig. A.7 Closeup 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 (SoriaRuiz et al. 2013, Sánchez Contreras et al. 2004): D = 1000 pc, F(4.7 μ) = 15 Jy, V = 15 km s^{1}, r = 10^{16} cm, relevant densities ~ 5 × 10^{6} 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 lowJ 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 RaighleighJeans equivalent brightness temperature of ^{12}CO ^{13}CO J = 2 − 1 transition:
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 ^{12}CO 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 . Some of the transitions have been scaled up or down by the specified factors for displaying purposes. The nonlinear dependance on the density n can easily be seen here, as well as the almost linear dependance on the product . 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 ^{12}CO 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 . Some of the transitions have been scaled up or down by the specified factors for displaying purposes. The nonlinear dependance on the density n can easily be seen here, as well as the almost linear dependance on the product . 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 arcsec^{2}, 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, dotteddashed line) and observations (black histogram) for the ^{12}CO and ^{13}CO transitions detected in NGC 6302 with HIFI. I_{factor} 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 ^{12}CO J = 2 − 1 transition. The size of each frame is 20 × 20 arcsec^{2}. 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 DinhVTrung 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 ^{12}CO J = 2 − 1 transition. The size of each frame is 20 × 20 arcsec^{2}. 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 DinhVTrung 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 T_{R} 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 T_{R} 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
Closeup of case B) for an ϵ value of 0.1. 

In the text 
Fig. A.4
Model calculations of T_{R} 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 T_{R} 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 T_{R} 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
Closeup 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.