Issue 
A&A
Volume 616, August 2018



Article Number  A131  
Number of page(s)  9  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201731943  
Published online  04 September 2018 
MOLPOPCEP: an exact, fast code for multilevel systems
^{1}
Instituto de Astrofísica de Canarias, 38205, La Laguna, Tenerife, Spain
email aasensio@iac.es
^{2}
Departamento de Astrofísica, Universidad de La Laguna, 38205 La Laguna, Tenerife, Spain
^{3}
Department of Physics & Astronomy, University of Kentucky, Lexington, KY 40506, USA
^{4}
Astronomy Dept., University of California, Berkeley, CA 94720, USA
Received:
13
September
2017
Accepted:
5
June
2018
We present MOLPOPCEP, a universal line transfer code that allows the exact calculation of multilevel line emission from a slab with variable physical conditions for any arbitrary atom or molecule for which atomic data exist. The code includes error control to achieve any desired level of accuracy, providing full confidence in its results. Publicly available, MOLPOPCEP employs our recently developed coupled escape probability (CEP) technique, whose performance exceeds other exact methods by orders of magnitude. The program also offers the option of an approximate solution with different variants of the familiar escape probability method. As an illustration of the MOLPOPCEP capabilities we present an exact calculation of the Spectral Line Energy Distribution (SLED) of the CO molecule and compare it with escape probability results. We find that the popular largevelocity gradient (LVG) approximation is unreliable at large CO column densities. Providing a solution of the multilevel line transfer problem at any prescribed level of accuracy, MOLPOPCEP is removing any doubts about the validity of its final results.
Key words: radiative transfer / line: formation / ISM: lines and bands / methods: numerical
© ESO 2018
1. Introduction
Much of the information about any astronomical source comes from its spectral line emission. Lines are the only probe of detailed kinematics, and provide the tightest constraints on density and temperature. Determining the spectral line emission from a multilevel system requires solution of the level population equations for all levels coupled with the radiative transfer equation for every line connecting them. Because of the complexity and computational demands of exact solution methods, many analysis codes bypass altogether solution of the radiative transfer equation, employing instead the escape probability approximation. This approach requires uniform physical conditions and is predicated on the conjecture that the effects of radiative transfer can be lumped into a multiplicative correction to the spontaneous decay rate. This single “escape probability” multiplicative factor is supposed to mimic the effects of radiative transfer in the entire source, and its functional form is posited from some plausibility arguments. Only the level populations are considered, calculated from rate equations augmented by these photon escape factors.
The escape probability approach amounts to an uncontrolled approximation without internal error estimate because it is not derived from first principles but instead is founded on a plausibility assumption right from the start. The only way to assess its error is to repeat the calculation with an exact method and compare the results^{1}. Nevertheless, this inherent shortcoming is often tolerated because of the simplicity and usefulness of the escape probability approach and the near impracticality of exact methods. Indeed, current calculations of line emission from active galactic nuclei, shock fronts and photodissociation regions (PDRs) are only performed in the escape probability approximation, which requires the underlying physical conditions to be uniform across the emitting region. Therefore these calculations are forced to replace in each case the variable conditions with a single value derived from either averaging or from plausibility arguments that attempt to pick out the most significant region. Even under these assumptions, the final results are inexact and there is no way of finding out how large, or how small, the errors involved are. Virtually all model results reported in the literature are afflicted by these problems. Furthermore, the escape probability is plainly useless in analysis of spectral line shapes—this method predicts flattop profiles for all optically thick quiescent regions (without largescale ordered motions), but such sources can in fact produce doublepeaked profiles (Elitzur et al. 2012, and references therein). Because of this inherent shortcoming, as long as the analysis is limited to escape probability calculations, much of the available information cannot be extracted from spectral data. This is an especially severe handicap when angular resolution is limited and line profiles provide the only handle on detailed structure.
We have recently developed the coupled escape probability (CEP), a new radiative transfer method that provides an exact solution for the multiline problem while retaining all the advantages of the naive escape probability approach (Elitzur & Asensio Ramos 2006; CEP06 hereafter). In this new technique the source is divided into zones and formal level population equations, including interactions with the transferred radiation, are derived rigorously from first principles. The final formulation does not contain the radiative transfer equation explicitly, only a set of coupled nonlinear level population equations. These equations are identical in form to those employed in standard escape probability calculations, but the naive photon escape factors are replaced by terms derived formally from the exact equations, including those for the radiative transfer of all lines. These terms introduce coupling between all zones, and it is this coupling which makes CEP an exact method. Solution of this set of algebraic equations determines in each zone level populations that are selfconsistent with the line radiation they generate. Once the correct level populations are derived, line profiles and fluxes are computed from straightforward summations over zones. The numerical error is controlled by reducing the zone sizes. Any level of accuracy can be achieved by increasing the number of zones.
Originally formulated for the slab geometry, the CEP method was subsequently extended to 3D, first to spherical sources both in hydrostatic equilibrium and subject to large scale motions (Yun et al. 2009; Yun & Park 2012). As in the slab case, the spherical implementation of CEP yields significant gains in accuracy and efficiency. Gersch & A’Hearn (2014) and Debout et al. (2016) have further extended CEP to asymmetrical situations and employed it to model optically thick line emission from cometary comae. The power of the CEP approach has been recognized also outside the astronomical community: Masnavi et al. (2013) employed CEP to model in detail extremeUV laser based lithography, an emerging technology in highvolume semiconductor chip production. And Rezaei et al. (2013) employed CEP in numerical calculations of the temporal behavior of plasma emission after laser irradiation.
Since CEP brings a significant speed improvement over other methods, exact calculations are becoming practical for most line analysis problems. We have implemented our original method in the code MOLPOPCEP and it is now publiclyavailable (see Sect. 3). MOLPOPCEP is a universal line transfer code that allows the exact calculation of multilevel line emission from a slab with variable physical conditions for any arbitrary atom or molecule for which atomic data exist. The code includes error control to achieve any desired level of accuracy, providing full confidence in its results.
MOLPOPCEP has already been employed in numerous studies, including analysis of the H_{2}D^{+} resonance line emission from protoplanetary disks (Asensio Ramos et al. 2007); doublepeaked line profiles in rotating disks (Elitzur et al. 2012); maser observations in supernova remnants (Pihlström et al. 2014; McEwen et al. 2014, 2016); farinfrared tracers of oxygen chemistry in diffuse clouds (Wiesemeyer et al. 2016); and ALMA multipletransition molecular line observations of an Ultraluminous Infrared Galaxy (Imanishi et al. 2017). The code is available from a dedicated web site^{2}.
This paper serves to introduce MOLPOPCEP and its capabilities. For completeness, in Sect. 2 we describe briefly the CEP radiative transfer formalism and introduce our notations for the key quantities. Section 3 describes the working of the program and discusses its implementation. A specific example is provided in Sect. 4, presenting detailed calculations of CO line emission and comparisons of the exact results with various variants of the escape probability. The paper concludes in Sect. 5 with suggested future directions of astrophysical radiative transfer modeling.
2. The line transfer problem
Consider a molecular or atomic multilevel system, with levels k = 1, 2, …, L ordered by energy. Determining the spectral line emission from this system requires the coupled solution of the level population equations and the radiative transfer equations for all the lines connecting them. In the CEP approach, all radiative quantities are expressed in terms of the level populations through the formal solution of the radiative transfer equation. As a result, the complete problem is formulated purely in terms of a set of nonlinear, selfconsistent level population equations. Here we summarize briefly the CEP solution formalism in the slab geometry; for full details, see CEP06.
2.1. CEP formalism in a slab
Consider the planeparallel geometry, so that physical properties vary only perpendicular to the surface. The slab is divided into z zones, sufficiently small that all properties can be considered constant within each zone. The population per substate of level k in zone i (=1, 2, …, z) is , where subscripts denote levels and superscripts zones. Then the system overall population in zone i is(1)
where g_{k} is the level degeneracy. Denote by ℓ^{i} the width of the ith zone and by Φ_{i}(ν) the absorption profile in the zone, normalized through ∫ Φ_{i}(ν)dν = 1 and assumed to have the same functional form for all transitions. For a transition between lower level l and upper level u with energy separation E_{ul} = hν_{ul} and Doppler width , introduce the dimensionless frequency shift from line center ; the linewidth Δν^{i} may reflect either thermal motions or the dispersion of a microturbulent velocity field. Then the zone optical thickness at frequency x along a ray slanted at θ = cos^{−1}μ from normal is where(2)
and B_{ul} is the coefficient of stimulated emission for the transition. Denote by A_{ul} the corresponding spontaneous transition rate and by the collision rate in zone i. Then the CEP level population equations are(3)
where d/dt = 0 in steady state. Here(4)
and where and , with the function β defined from(6)
This function was first introduced by Capriotti (1965), who also provided numerical approximations that we found useful for its accurate, efficient calculation. Plotted in Fig. 1 for the Doppler profile, Φ(x) = π^{−1/2} exp(−x^{2}), the function β contains the essence of radiative transfer in the CEP approach. It represents the probability for photon escape from a slab of optical thickness τ, averaged over the injected photon direction, frequency and position in the slab. The first term in the expression defining p^{i} (Eq. (4)) is thus the average probability for photon escape from zone i, reproducing one of the common variants of the escape probability method in which the whole slab is treated as a single zone (e.g., Krolik & McKee 1978). The subsequent sum in Eq. (4) describes the effect on the level populations in zone i of radiation produced in all other zones. Each term in the sum has a simple interpretation in terms of the probability that photons generated elsewhere in the slab traverse every other zone and get absorbed in zone i, where their effect on the level populations is similar to that of external radiation.
Fig. 1.
Left: optical depth variation of β (Eq. (6)), plotted in green for the Doppler profile; this is the fundamental CEP function describing the average probability for escape from a uniform slab with optical depth τ across its full width (see text). Also plotted are escape probability approximations for a uniform, static sphere (Eq. (13)) and LVG approximations for spherical (Eq. (14)) and planeparallel (Eq. (15)) geometries. In the LVG cases, the appropriate velocity gradient replaces the ratio Δυ/ℓ in the definition of optical depth. Right: the ratio of each escape probability approximation and the fundamental slab function. 
When external continuum radiation exists, each term in the sums on the righthandside of Eq. (3) is supplemented by , where is the profile and angleaveraged intensity of the external radiation in zone i. When the external radiation corresponds to the emission from dust permeating the source, is simply the angleaveraged intensity of the local dust emission in the ith zone. When the external radiation originates from outside the slab and has an isotropic distribution with intensity I_{e} (=J_{e}) in contact with the τ = 0 face,(7)
When the slab is illuminated by parallel rays with intensity I_{e} (=4πJ_{e}) entering at direction (μ_{0}, ϕ_{0}) to the τ = 0 face,(8)
Equation (3) provides a set of L − 1 independent equations for the L unknown populations in each zone, . Equation (1) for the overall density in the zone closes the system of equations. It is convenient to switch to the scaled quantities as the unknown variables and introduce the overall column density(10)
Neither densities nor physical dimensions need then be specified since only 𝒩 enters as an independent variable, with the zone partition done in terms of 𝒩 rather than ℓ. Apart from external radiation, the problem is fully specified by three input properties: temperature and density of collision partners, which together determine the collision terms, and 𝒩, which sets the scale for all optical depths.
Solution of the set of Eqs. (1) and (3) yields the full solution of the transfer problem for all lines by considering only level populations; the computed populations are selfconsistent with the radiation field in the slab, including the internally generated diffuse radiation, even though the radiative transfer equation is not handled at all. Once the populations are found, radiative quantities can be calculated in a straightforward manner from summations over the zones. For example, the contribution of the u → l transition to the slab cooling rate per unit area, accounting for the emission from both faces of the slab, is determined by its surface flux F_{ν,ul} and can be characterized by the line cooling coefficient(11)
introduced for convenience when Δν is constant in the slab. Once the level populations have been determined, j can be calculated from(12)
where is the line source function in the ith zone.
2.2. Solution
Solution of the overall system of z · L nonlinear algebraic equations (1) and (3) in all zones determines level populations that are consistent with the radiation field generated everywhere inside the slab. In actual numerical calculations, the solution of these equations provides the exact solution of the radiative transfer problem for all levels when for every i; the only approximation is the finite size of the discretization, i.e., the finite number of zones. Therefore to solve the problem at any desired precision, start with some initial number of zones and keep refining the divisions until the relative change in level populations decreases everywhere below the prescribed tolerance.
Because of the nonlocal nature of radiative transfer, the system of equations is highly nonlinear, making it necessary to apply suitable iterative methods. We find the multidimensional Newton method most suitable for solution of the nonlinear CEP equations. The efficiency of the Newton method is enhanced in the CEP approach because the Jacobian matrix can be calculated analytically, since the functional dependence on population is known explicitly for all terms. Because the Newton method requires inversion of the Jacobian matrix, the number of operations in this process increases as the third power of the matrix dimension and can degrade the performance in cases of very large numbers of levels and zones. Matrix inversion is avoided in the BiCGSTAB iterative scheme designed by Van der Vorst (1992) for solution of the linear system of the Newton method. In this scheme, geared toward sparse Jacobian matrices, only the nonzero matrix elements are stored and used. It is particularly suitable for the CEP technique because multilevel problems tend to produce sparse matrices, as each level generally couples to only a limited number of other levels. MOLPOPCEP switches automatically to this method when the matrix size exceeds 10^{3} × 10^{3}. Matrix inversion can also be avoided by adopting a Λiteration approach, instead of the Newton method, to solving the set of nonlinear level population equations. As described in CEP06 (see Sect. 5.2 in that paper), the CEP method is well suited for acceleration schemes and this approach can be selected in MOLPOPCEP when the Newton method slows down.
Whatever the iteration technique, it is always advantageous to start from a good initial guess. To this end we have implemented two variants of a simple strategy. The first is to start from the 𝒩 → 0 limit with p = 1 everywhere, then the level population equations are linear (see Eq. (3)). The solution of this set of linear equations serves as the initial guess for a low column density 𝒩 in which all lines are optically thin. Then N is increased in small steps, with the previous solution taken as the initial guess for the increased column. The steps are repeated until the desired column density is reached. The complementary approach is to start from the 𝒩 → ∞ limit, which yields the Boltzmann distribution, use that as the initial guess for a very large column in which all lines are optically thick and decrease 𝒩 toward the target value. Going either way, this incremental strategy aids convergence and provides the solutions for many intermediate cases as a byproduct.
3. The computer code MOLPOPCEP
MOLPOPCEP is based on a code originally developed by M. Elitzur for solving the nonLTE level populations of an arbitrary species in the escape probability approximation. The CEP method has been fully implemented into the original code, which has been further modularized and ported to Fortran 90. For a slab with variable physical conditions, MOLPOPCEP can calculate at a prescribed level of accuracy the line emission by any atom or molecule for which atomic data exist. We outline the program capabilities through a brief description of its input and output.
3.1. Basic input
The input file has a free format, text and empty lines can be entered arbitrarily. All lines that start with the ‘*’ sign are copied to the output, and can be used to print out notes and comments. This option can also be useful when the program fails for some mysterious reason, enabling the user to compare its output with an exact copy of the input line as it was read in before processing by MOLPOPCEP.
A single MOLPOPCEP run can process an unlimited number of models. To accomplish this, the program always calls a master input file (molpop.inp) containing a list of the individual cases that are launched sequentially. These individual models can reside in separate directories, with each model output produced in the corresponding directory. The information in every input file can be roughly divided into six groups, which we now describe.
3.1.1. Radiative transfer method
MOLPOPCEP provides a full CEP calculation in the slab geometry. Although the CEP formalism can handle any profile Φ, currently only the Doppler profile is implemented. This is sufficient for most cases.
Through the use of keywords, MOLPOPCEP also offers the choice of a standard escape probability approximation, included because it provides handling of three special situations that have not yet been implemented in CEP:
Line overlap. Under certain circumstances, linewidths become comparable to the separation between different lines so that photons emitted in one transition can be absorbed in another. One commonly affected molecule is OH (e.g., Guilloteau et al. 1981), and this effect has been shown to explain the differences between the emission patterns of OH megamasers and their Galactic counterparts (Lockett & Elitzur 2008). The effect is handled with the method of Lockett & Elitzur (1989) and can be turned on with a proper flag when MOLPOPCEP is run in the escape probability mode. Other molecules that have line overlap effects, due to the presence of hyperfine structure, are HCN (Cernicharo et al. 1984; GonzálezAlfonso & Cernicharo 1993) and N_{2}H^{+} (Daniel et al. 2006). Another instance of line overlap appears when lines from different species coincide in wavelength (e.g., Cernicharo et al. 1991; Cernicharo & Bujarrabal 1992; GonzálezAlfonso et al. 1996). Such an effect, not currently included in MOLPOPCEP, could be implemented with relatively minor modifications.
Dust contribution. Molecular gas is mixed with dust, therefore line photons can be absorbed by the dust when it is optically thick at the same wavelength. Such absorption and the corresponding emission affect the line transfer (see, e.g., Elitzur 1992, pp. 202–203). These modifications are incorporated in the code and can be turned on by setting the appropriate flag.
Maser saturation. Transitions with inverted population (n_{u} > n_{l}) produce maser radiation, which can strongly affect the population inversion when the maser becomes saturated. The saturation effect, strongly dependent on the geometry because of maser beaming, can be described with the escape probability approach (Elitzur 1990). In MOLPOPCEP, saturation can be either neglected (the unsaturated maser domain) or handled approximately with the escape probability in a rudimentary manner that ignores maser beaming; a proper treatment of maser saturation requires a more accurate handling of the beaming effect (e.g., Daniel & Cernicharo 2013).
When run in the escape probability mode, MOLPOPCEP offers four variants of this approximation. The first is a uniform, static slab, with the escape probability β given in Eq. (6). While this is exactly the same as performing a CEP calculation with a single zone, the special effects just described are only available when the escape probability option is selected. Another variant is the escape probability from a uniform, static sphere(13)
where τ is the optical depth across the diameter (van der Tak et al. 2007; note that this is the escape probability used in their online RADEX code). The two final variants involve approximate solutions for large velocity gradients (LVG). MOLPOPCEP implements the escape probability for radial flows in spherical geometry, with the local velocity variation ε = d ln υ/d ln r an input parameter in the range 0 < ϵ ≤ 1 (see, e.g., Elitzur 1992, pp. 39–41). In the literature, the term “LVG” invariably refers to a Hubbletype velocity field (υ ∝ r; ϵ = 1), which yields(14)
Here the optical depth τ is given by Eq. (2), with the ratio Δυ/ℓ replaced by the radial velocity gradient dυ/dr. MOLPOPCEP also handles velocity gradients in planeparallel geometry, implemented with the Scoville & Solomon (1974) expression(15)
where in this case the relevant velocity gradient is taken normal to the plane.
Figure 1 plots the various escape probabilities as functions of optical depth and compares them with the slab fundamental function. The differences from the slab case can be large. Even for the LVGPP case, that uses the same geometry, the ratio β_{LVG−PP}/β_{slab} is significantly different from unity at all τ > 0. Furthermore, the ratio keeps decreasing as τ is increasing, as do the ratios for other escape probabilities. The reason is that at large τ, the various escape probability approximations vary as 1/τ while the proper escape probability for a slab decreases more slowly as . The plausibility arguments employed to derive the various escape probability approximations are too crude to capture the logarithmic dependence of radiative transfer, rendering them increasingly unreliable as τ increases.
3.1.2. The molecule/atom
An atomic or molecular species is defined by the energy levels, statistical weights and Acoefficient tabulated in an ordinary text file. Each data file is identified by its name and directory, allowing for comparison of different sets of data for the same species. The MOLPOPCEP distribution package provides a set of entries from the Leiden Atomic and Molecular Database (LAMDA)^{3}, and it also includes a tool that downloads and updates the MOLPOPCEP input files with the latest entries from this database. As an example of multiple databases, we include also sample data files from the Basecol^{4} database. Installing additional species requires only expansion of the MOLPOPCEP atomic database; the code itself remains untouched.
3.1.3. Collisions
MOLPOPCEP allows for collisions with up to ten collisional partners. Tabulated collision rates are invoked through their file names. Since collisional data are generally available only for a small set of temperatures, MOLPOPCEP offers a number of different interpolation methods. Because rate coefficients are not always available, the code offers some simple analytic approximations (e.g., hard sphere collisions) that can be invoked with keywords.
3.1.4. Physical conditions
Uniform physical conditions (density, temperature and molecular abundance) can be specified in the input file. This is the only option when MOLPOPCEP is run in the escape probability mode. The full CEP mode allows also variable conditions in the slab, in which case the spatial profiles of the physical parameters have to be tabulated in a separate file. The number of CEP zones must then be at least as large as the number of entries in these profiles.
3.1.5. External radiation
The cosmic blackbody radiation at a temperature of 2.7 K is always included. Additional radiation fields can also be specified, and in the slab case it is possible to illuminate each side with different radiation. The external radiation can include an arbitrary number of diluted blackbody components, each parameterized in terms of its temperature and dilution factor; this option covers all cases of illumination by nearby stars. One can also specify an arbitrary number of radiation fields of the form (1 − e^{−τλ})B_{λ}(T_{d}), where B is the Planck function and τ_{λ} optical depth with the spectral variation of interstellar dust. This can be used as a crude approximation for emission by dust, characterized by its temperature T_{d} and optical depth at visual, that is mixed with the gas. Additionally, it is possible to specify a radiation field with an arbitrary spectral shape through a tabulation in a file. Such tabulations can be generated by some other radiative transfer code such as, e.g., DUSTY (Ivezic et al. 1999).
3.1.6. Numerics
MOLPOPCEP offers great control over all aspects of its operation, including various accuracy parameters. In the multizone case, the solution technique can be chosen between Newton method with analytical derivatives (both matrix inversion and noninversion) and an accelerated Λiteration. Convergence is tested with(16)
where n(itr) refers to a vector containing populations for all the levels in all zones at iteration number itr. Convergence is attained when R_{c}(itr) decreases below an inputspecified threshold. For safety, a maximum number of iterations is also specified to stop execution in case of a runaway calculation. An orderly completion of a run occurs when the column density reaches the maximum/minimum column specified in the input for the scheme that iteratively increases/decreases the column density to improve the convergence properties (Sect. 2.2). If the CEP method is used, it is also possible to turn on the grid convergence. The problem is solved in grids of increasing number of zones until the maximum relative change between two consecutive grids decreases below the userdefined threshold.
3.2. Output and its control
A single MOLPOPCEP run generates prodigious amounts of data, more than might be needed in any particular application of the code. Various on/off switches provide control over the output. The numerous options are described in detail both in the manual and in the sample input files provided with the MOLPOPCEP package. As noted above, every input line that starts with * is echoed in the output, making is easy to incorporate notes and comments.
4. An illustrative example
As a demonstration of the capabilities of MOLPOPCEP we present here a detailed analysis of CO emission, calculating what has become known as its Spectral Line Energy Distribution (SLED; Weiß et al. 2005; Papadopoulos et al. 2012). Several transitions of the CO rotational ladder are observable in highredshift galaxies (e.g., Solomon & Vanden Bout 2005). Analysis of diagrams that plot line emission vs. rotational quantum number, i.e., CO SLED, offers one of the most direct diagnostic techniques to study the excitation conditions of the molecular gas. In LTE, the level populations depend only on the local temperature T, while the line emission depends additionally on the CO column density per unit band width. In that case, the shape of the CO SLED normalized to the J = 1 → 0 emission uniquely determines T, and the CO column follows directly from the magnitude of the measured flux. Outside of LTE, the level population distribution is determined through the competition between collisional and radiative interactions, bringing in also the local H_{2} number density and necessitating a full radiative transfer calculation.
To our knowledge, the CO SLED has only been computed in the LVG escape probability approximation (Eq. (14)). Inherently, such calculations cannot be exact even when the physical conditions are uniform: The strength of radiative interactions varies with distance to the surface for optically thick lines, therefore the population distribution of transition levels does depend on position in the source. This effect cannot be captured by the escape probability approach, where the entire source is described by a single excitation temperature. Furthermore, every line on the CO rotational ladder has a different optical depth, therefore they are affected differently by this variation. As a result, the reliability of studies based on SLED analysis (e.g., Papadopoulos et al. 2012; da Cunha et al. 2013) cannot be assessed because the escape probability method does not contain an error estimate so it is impossible to know whether its deviations from the correct results fall within the observational uncertainties. An example where the conclusions based on escape probability analysis are in serious error is the ratio of the two ^{3}P oxygen cooling lines at 63 and 145 μm (see CEP06). As is always the case, it is impossible to know whether the escape probability results are meaningful without comparison with an exact calculation.
Here we perform such a comparison for the CO SLED. For a meaningful comparison, one must utilize the same computed quantity. While the CEP method provides rigorous expressions for all radiative quantities, the escape probability approximation can be expected to produce useful results only in some average sense. This makes the line cooling coefficient j(Eq. (11))especially useful. In the case of a slab, the CEP formalism yields Eq. (12). Since a singlezone CEP calculation is equivalent to a slab escape probability calculation, we can read off immediately the proper escape probability expression by inserting z = 1. This yields the familiar escapeprobability result j_{ul} = hν_{ul}A_{ul}β_{ul}N_{u}/4πΔν_{ul}, where N_{u} is the column density of the upper level (e.g., Elitzur 1992, p 28).
We have conducted an extensive set of MOLPOPCEP calculations for a uniform slab, varying the temperatures from 50 K to 150 K and the molecular hydrogen density from 10^{2} to 10^{6} cm^{−3}. The slab thickness, linewidth and CO abundance, [CO] = n(CO)/n(H_{2}), enter only through the combination γ = [CO] ℓ/Δν; the scale of all optical depths is controlled by the product γn (Eq. (2)). Beginning with Weiß et al. (2005), CO SLED calculations assumed a constant value for this parameter, and for compatibility we adopt this assumption, presenting results for γ = 10^{−4} and 10^{−5} pc (km s^{−1})^{−1}; these values bracket the likely range in local and highz galaxies (Weiß et al. 2007). Collision rates are from Yang et al. (2010) with orthotopara H_{2} ratio of 3. The only external radiation is the cosmic microwave background.
4.1. Exact solutions
Figures 2 and 3 show the resulting SLEDs in terms of 4π j, the line cooling coefficient in units of flux density, for the two studied values of γ. Density varies among panels in each column, temperature between the columns. The exact solutions are shown with solid blues lines; we have verified that in all cases the CEP calculations with 20 zones are fully converged—increasing further the number of zones does not change the results. It is interesting to note how different the actual SLEDs are from LTE predictions. If all models were in LTE, the SLEDs would be identical in each column and would differ only between the columns. The actual behavior is just the opposite: CO SLEDs vary only moderately with temperature in the displayed range but are strongly dependent on density. Also, some combinations of T and n(H_{2}) produce inversions in low J → J−1 transitions, as noted already by Goldreich & Kwan (1974). These cases were ignored in the figures.
Fig. 2.
MOLPOPCEP results for the CO SLEDs of uniform slabs, computed with different methods (the line cooling coefficient j is defined in Eq. (11)). The yaxis scale of each panel is listed at its topleft corner. Columns have the same temperature, listed at the top, rows the same H_{2} density, marked in each leftmost panel. Optical depths are set from [CO] ℓ/Δυ = 10^{−4} pc (km s^{−1})^{−1}, where [CO] is the CO abundance. Full CEP calculations using grids with 10 and 20 zones are plotted with solid lines of different colors, as marked. When the 20zones line (blue) is missing, it is because it coincides with the 10zones result (red), demonstrating CEP convergence in practically all cases. Dashed lines show results for different escape probabilities (see Fig. 1). Note that the slab escape probability calculation is the same as singlezone CEP. 
Each panel also shows CEP calculations with 10 zones (solid red lines). Doubling the number of zones from 10 to 20 hardly affects the results, in many cases the corresponding plots are barely distinguishable from each other. In Fig. 2, where γ = 10^{−4}, the largest differences from the exact solution are at n(H_{2}) = 10^{4} cm^{−3}, declining towards both smaller and larger densities. Even though CO is a multilevel system, these trends can be understood from the solution of the twolevels problem: As shown in CEP06, the largest deviations of the slab escape probability (i.e., single zone) calculation from the exact solution of the twolevels problem occur around τn/n_{crit} ~ 1, where n_{crit} is the transition critical density. When γ = 10^{−4}, this condition is fulfilled around the SLED peak at about n(H_{2}) = 10^{4} cm^{−3}, and that is where the largest deviations occur from the exact solutions. Since τ ~ γn (Eq. (2)) and γ is held fixed, τn ∝ n^{2} along each column, and the deviations from the exact solution are rapidly diminished moving away from n(H_{2}) = 10^{4} cm^{−3} in either direction. In Fig. 3, where γ = 10^{−5}, optical depths are smaller and the 10zone results are practically identical to the exact ones in all cases. Thus the CEP calculations with only 10 zones are practically convergent in all the presented solutions.
4.2. Escape probability calculations
Solutions obtained with various escape probability approximations are shown with dashed lines. In Fig. 2, virtually all escape probability approximations reproduce rather well the exact solutions in the top two rows, where optical depths are the smallest. The sole outlier is the slab LVG approximation (LVGPP; Eq. (15)). Even though this approximation purports to describe the very same slab geometry, it fails to reproduce its proper solution even at the smallest density, and its deviations from the exact results keep increasing continuously with density (and optical depths). Evidently, the plausibility arguments that yielded this approximation (Scoville & Solomon 1974) were off; the planeparallel LVG escape probability is inadequate for CO radiative transfer calculations.
The slab escape probability approximation (Eq. (6)), equivalent to a singlezone CEP calculation, is doing much better, producing reasonably good agreement with the exact solution. Similar to the 10zone CEP calculation, and for the same reasons (Sect. 4.1), the largest deviations from the exact solution occur around n(H_{2}) = 10^{4} cm^{−3} when γ = 10^{−4} (Fig. 2). Yet even in this case the 1zone calculation is within ~27% of the exact solution, and its deviations decrease towards both lower and higher densities. When γ = 10^{−5} (Fig. 3), the slab escape probability similarly provides an adequate approximation in all cases except for the lowest density of 10^{2} cm^{−3}, where the discrepancy is more than a factor of 2. Since λ is a factor of 10 lower than in Fig. 2, these models with n(H_{2}) = 10^{2} cm^{−3} stand out as having the lowest optical depths of all the cases presented here, implying that large discrepancies occur only at low optical depths. All in all, the results show that the slab escape probability provides an excellent approximation for CO SLEDs in all uniform slabs with n(CO)ℓ/Δυ ≳ 10^{−3} cm^{−3} pc (km s^{−1})^{−1}.
The figures also show the results of calculations with the two escape probabilities for spherical geometry offered by MOLPOPCEP: a uniform, static sphere (Eq. (13)) and LVG spherical expansion (Eq. (14)). Both versions provide excellent approximations for the exact results at densities n(H_{2}) ≤ 10^{4} cm^{−3} when γ = 10^{−4}. Beyond that, the deviations from the exact solution increase with density (and optical depths), with the SLED peaks occurring at the wrong J. These large deviations cannot be attributed to the effect of a different geometry on radiative transfer—the two escape probabilities are successful at n(H_{2}) = 10^{4} cm^{−3}, where optical depths are already significant. Rather, these versions of β are afflicted by the same weakness as the LVGPP escape probability, which produces the same SLED peaks: all three yield β ~ 1/τ at large τ, failing to reproduce the additional logarithmic variation, which becomes important when τ is sufficiently larger. The various escape probability approximations are too crude to capture the logarithmic dependence of radiative transfer, losing their reliability at large optical depths. As seen in Fig. 3, all escape probabilities are off at n(H_{2}) = 10^{2} cm^{−3} when γ = 10^{−5}. Therefore, the common LVG escape probability, as well as the staticsphere escape probability, provides an adequate approximation for CO radiative transfer only in the range 10^{−3} cm^{−3} pc (km s^{−1})^{−1} ≳ n(CO)ℓ/Δυ ≲ 10^{−1} cm^{−3} pc (km s^{−1})^{−1}; the utility range of the two spherical escape probabilities is limited both at small and large CO column densities.
5. Discussion
Highresolution facilities such as ALMA enable detailed observations that can only be understood through realistic 3D modeling. However, on top of complexity and heavy computational demands, proper definitions of 3D models involve a large number of uncertainties. While such models are clearly necessary in many cases, the widespread success of the RADEX code (as of this writing, van der Tak et al. 2007 has received more than 700 citations) shows there remains a large demand for efficient, if simple, radiative transfer tools. MOLPOPCEP offers such a tool, providing an efficient, exact and verifiable solution of radiative transfer.
While species such as N_{2}H^{+}, HCN, HCO^{+}, CS or NH_{3} are required for probing densities higher than about 10^{4} cm^{−3}, the two most common indicators of physical conditions in molecular clouds and PDRs are CO, whose SLED analysis is presented here, and OI cooling lines, studied earlier in CEP06. In both cases LVG escape probability calculations can produce erroneous results, therefore this widespread technique is an unreliable analysis tool. Although the escape probability is occasionally adequate in some limited regions of parameter space, it is impossible to discern any trends or establish guidelines that would enable reliance on it alone for calculating line intensities for multilevel systems. The only way to get the correct result is to perform a correct calculation; there is no viable shortcut through the escape probability approximation. MOLPOPCEP eliminates the guesswork, offering a public tool for exact solution of the line transfer problem. While the current version handles only slabs, this geometry suffices for accurate modeling of both shock and photodissociation fronts, covering the needs of many interstellar studies.
Dumont et al. (2003) provide a detailed discussion of the escape probability method and comparison with ALI calculations.
Acknowledgments
We thank Paul Goldsmith, Massa Imanishi and José Cernicharo, the referee. for useful comments and suggestions. AAR acknowledges financial support by the Spanish Ministry of Economy and Competitiveness (MINECO) through project AYA201460476P.
References
 Asensio Ramos, A., Ceccarelli, C., & Elitzur, M. 2007, A&A, 471, 187 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Capriotti, E. R. 1965, ApJ, 142, 1101 [NASA ADS] [CrossRef] [Google Scholar]
 Cernicharo, J., & Bujarrabal, V. 1992, ApJ, 401, L109 [NASA ADS] [CrossRef] [Google Scholar]
 Cernicharo, J., Guelin, M., & Askne, J. 1984, A&A, 138, 371 [NASA ADS] [Google Scholar]
 Cernicharo, J., Bujarrabal, V., & Lucas, R. 1991, A&A, 249, L27 [NASA ADS] [Google Scholar]
 da Cunha, E., Walter, F., Decarli, R., et al. 2013, ApJ, 765, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Daniel, F., & Cernicharo, J. 2013, A&A, 553, 70 [Google Scholar]
 Daniel, F., Cernicharo, J., & Dubernet, M.L. 2006, ApJ, 648, 461 [NASA ADS] [CrossRef] [Google Scholar]
 Debout, V., BockeléeMorvan, D., & Zakharov, V. 2016, Icarus, 265, 110 [NASA ADS] [CrossRef] [Google Scholar]
 Dumont, A.M., Collin, S., Paletou, F., et al. 2003, A&A, 407, 13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Elitzur, M. 1990, ApJ, 363, 638 [Google Scholar]
 Elitzur, M. 1992, Astronomical Masers (Dordrecht: Kluwer Academic Publishers) [Google Scholar]
 Elitzur, M., & Asensio Ramos, A. 2006, MNRAS, 365, 779 [Google Scholar]
 Elitzur, M., Asensio Ramos, A., & Ceccarelli, C. 2012, MNRAS, 422, 1394 [NASA ADS] [CrossRef] [Google Scholar]
 Gersch, A. M., & A’Hearn, M. F. 2014, ApJ, 787, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Kwan, J. 1974, ApJ, 189, 441 [NASA ADS] [CrossRef] [Google Scholar]
 GonálezAlfonso, E., & Cernicharo, J. 1993, A&A, 279, 506 [NASA ADS] [Google Scholar]
 GonálezAlfonso, E., Alcolea, J., & Cernicharo, J. 1996, A&A, 313, 13 [NASA ADS] [Google Scholar]
 Guilloteau, S., Lucas, R., & Omont, A. 1981, A&A, 97, 347 [NASA ADS] [Google Scholar]
 Imanishi, M., Nakanishi, K., & Izumi, T. 2017, ApJ, 849, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Ivezic, Z., Nenkova, M., & Elitzur, M. 1999, ArXiv eprints [arXiv: astroph/9910475] [Google Scholar]
 Krolik, J. H., & McKee, C. F. 1978, ApJS, 37, 459 [NASA ADS] [CrossRef] [Google Scholar]
 Lockett, P., & Elitzur, M. 1989, ApJ, 344, 525 [NASA ADS] [CrossRef] [Google Scholar]
 Lockett, P., & Elitzur, M. 2008, ApJ, 677, 985 [NASA ADS] [CrossRef] [Google Scholar]
 Masnavi, M., Szilagyi, J., Parchamy, H., & Richardson, M. C. 2013, Appl. Phys. Lett., 102, 164102 [NASA ADS] [CrossRef] [Google Scholar]
 McEwen, B. C., Pihlström, Y. M., & Sjouwerman, L. O. 2014, ApJ, 793, 133 [NASA ADS] [CrossRef] [Google Scholar]
 McEwen, B. C., Pihlström, Y. M., & Sjouwerman, L. O. 2016, ApJ, 826, 189 [NASA ADS] [CrossRef] [Google Scholar]
 Papadopoulos, P. P., van der Werf, P. P., Xilouris, E. M., et al. 2012, MNRAS, 426, 2601 [NASA ADS] [CrossRef] [Google Scholar]
 Pihlström, Y. M., Sjouwerman, L. O., Frail, D. A., et al. 2014, AJ, 147, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Rezaei, F., Karimi, P., & Tavassoli, S. H. 2013, Appl. Opt., 52, 5088 [NASA ADS] [CrossRef] [Google Scholar]
 Scoville, N. Z., & Solomon, P. M. 1974, ApJ, 187, L67 [NASA ADS] [CrossRef] [Google Scholar]
 Solomon, P. M., & Vanden Bout, P. A. 2005, ARA&A, 43, 677 [NASA ADS] [CrossRef] [Google Scholar]
 van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Van der Vorst, H. A. 1992, SIAM J. Sci. Stat. Comput., 13, 631 [CrossRef] [MathSciNet] [Google Scholar]
 Weiß, A., Downes, D., Walter, F., & Henkel, C. 2005, A&A, 440, L45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weiß, A., Downes, D., Neri, R., et al. 2007, A&A, 467, 955 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wiesemeyer, H., Güsten, R., Heyminck, S., et al. 2016, A&A, 585, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yang, B., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2010, ApJ, 718, 1062 [NASA ADS] [CrossRef] [Google Scholar]
 Yun, Y. J., & Park, Y.S. 2012, A&A, 545, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yun, Y. J., Park, Y.S., & Lee, S. H. 2009, A&A, 507, 1785 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Figures
Fig. 1.
Left: optical depth variation of β (Eq. (6)), plotted in green for the Doppler profile; this is the fundamental CEP function describing the average probability for escape from a uniform slab with optical depth τ across its full width (see text). Also plotted are escape probability approximations for a uniform, static sphere (Eq. (13)) and LVG approximations for spherical (Eq. (14)) and planeparallel (Eq. (15)) geometries. In the LVG cases, the appropriate velocity gradient replaces the ratio Δυ/ℓ in the definition of optical depth. Right: the ratio of each escape probability approximation and the fundamental slab function. 

In the text 
Fig. 2.
MOLPOPCEP results for the CO SLEDs of uniform slabs, computed with different methods (the line cooling coefficient j is defined in Eq. (11)). The yaxis scale of each panel is listed at its topleft corner. Columns have the same temperature, listed at the top, rows the same H_{2} density, marked in each leftmost panel. Optical depths are set from [CO] ℓ/Δυ = 10^{−4} pc (km s^{−1})^{−1}, where [CO] is the CO abundance. Full CEP calculations using grids with 10 and 20 zones are plotted with solid lines of different colors, as marked. When the 20zones line (blue) is missing, it is because it coincides with the 10zones result (red), demonstrating CEP convergence in practically all cases. Dashed lines show results for different escape probabilities (see Fig. 1). Note that the slab escape probability calculation is the same as singlezone CEP. 

In the text 
Fig. 3.
Same as Fig. 2, only [CO] ℓ/Δυ = 10^{−5} pc (km s^{−1})^{−1}. 

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.