A&A 485, 127-136 (2008)
DOI: 10.1051/0004-6361:200809440

Radiative transfer revisited for emission lines in photon dominated regions

M. Gonzalez Garcia - J. Le Bourlot - F. Le Petit - E. Roueff

LUTH, Observatoire de Paris and Université Paris 7 Place Jansen, 92190 Meudon, France

Received 23 January 2008 / Accepted 6 March 2008

Abstract
Context. Transfer in lines controls the gas cooling of photon dominated regions (PDR) and provides many of the observational constraints that are available for their modelling.
Aims. The interpretation of infrared and radio observations by the new generation of instruments, such as Herschel, requires sophisticated line radiative-transfer methods. The effect of dust emission on the excitation of molecular species in molecular regions is investigated in detail to explicitly show the origin of various approximations used in the literature. Application to H_2O is emphasised.
Methods. The standard 1D radiative transfer equation is written as a function of the space variable (as opposed to the usual optical depth). This permits to simultaneously consider all pumping contributions to a multi-level species in a non-uniform slab of dust and gas. This treatment is included in the Meudon PDR Code (available at http://aristote.obspm.fr/MIS/).
Results. Infrared emission from hot grains at the edge of the PDR may penetrate deep inside the cloud, providing an efficient radiation source to excite some species at a location where cold grains no longer emit. This leads to non-negligible differences with classical escape probability methods for some lines, e.g. water. Cooling efficiency does not follow directly from line emissivities. The infrared pumping contribution leads to a higher excitation that enhances collisional de-excitation and reduces cooling efficiency.

Key words: ISM: general - line: formation - radiative transfer - methods: numerical - ISM: molecules

  
1 Introduction

The derivation of molecular column densities from molecular emissions requires detailed radiative transfer treatment of the various emitted transitions. A summary of the various radiative transfer methods can be found in van Zadelhoff et al. (2002). The large velocity gradient technique or escape probability method provides first-order estimates (van der Tak et al. 2007). A full radiative transfer treatment within a given assumed geometry provides significantly improved results, and various numerical methods have been developed such as the accelerated $\Lambda$-iteration method (ALI) (Hubeny 2001; Chevallier et al. 2003) and, more recently, the coupled escape probability method (CEP) (Elitzur & Asensio Ramos 2006), which aims to solve the radiative transfer equations exactly. Monte-Carlo techniques are also highly efficient and allow taking various geometries into account (Juvela 2005; Bernes 1979). The major difficulties are the computation time and the control of the resulting errors.

The aim of this paper is to compute the excitation of atomic and molecular species by considering excitation by collisions and pumping due to external radiation fields, as well as internal sources such as dust and lines.

This paper presents the full derivation of the detailed balance equations describing the steady-state level populations of any atomic or molecular system subject to collisions, external and local (dust and gas) radiation fields. We computed the exact expression of the mean radiation intensity $J_{\nu}$ in Sect. 2, where chemical reactions, inelastic collisions, and local emission/absorption effects are taken into account. We outline different approximations corresponding to three different contributions of the monochromatic intensity in Sect. 3. The energy released in line emissions following collisional excitation is one of the major cooling processes of the gas. However, the efficiency of the cooling depends on the ability of the radiation to escape, as described in Sect. 4.2. We apply the different methods and approximations described above to the case of the water molecule in Sect. 5.

As already discussed by other authors (van der Tak et al. 2006; Poelman & Spaans 2005; Poelman et al. 2007; Poelman & Spaans 2006), the complicated level structure of H_2O combined with its high permanent dipole moment requires a particularly careful analysis of its excitation and of its subsequent cooling emission transitions. The corresponding implementation in our PDR (photon dominated region) code has been explicitly achieved, complementing the upgrade of ultraviolet radiative transfer by Goicoechea & Le Bourlot (2007).

The case of S140, a well-documented PDR, is specifically addressed as a benchmark of our calculations and permits a detailed comparison with the previous work of Poelman & Spaans (2005, hereafter PS05). Finally, we summarise our results and display our conclusions in Sect. 6. Useful mathematical subtleties are outlined in separate appendices.

  
2 Line excitation

In a medium far from LTE, excitation of any given species can only be computed by solving detailed balance equations. Le Petit et al. (2006) have shown that these detailed balance equations should be written with chemical formation and destruction terms, thus leading to a fully determined set of equations that does not require the introduction of a conservation equation. For a level i, detailed balance gives

 

$$% n_{i}~\left[~ \sum_{j<i}~ A_{ij} +\sum_{j \neq i}~ \left( B_{ij}~\bar{J}_{ij} + \sum_{X}~k_{ij}^{X}~n^{X} \right)+D \right]=$$

$$\sum_{j>i}~ n_{j}~ A_{ji} + \sum_{j \neq i}~ n_{j}~\left( B_{ji}~\bar{J}_{ji} + \sum_{X}~k_{ji}^{X}~ n^{X} \right) + F_{i}$$ (1)

where n_i is the population of level i, A_ij and B_ijthe Einstein coefficients, k_ji^X the collision rate with species X, n^X the abundance of collision partner X, D the total destruction rate, and F_i the formation rate on level i. Both D and F_i are known from the steady state solution to chemical balance equations, where elemental conservation is incorporated.

To solve this system, the most difficult term to compute is the mean intensity  $\bar{J}_{ij}$, which controls radiative excitation. In the following, we deal with a purely 1D plane parallel cloud, as described in Le Petit et al. (2006). At a position $\tau$ within a cloud, $\bar{J}_{ij}$ is computed from

$$% J_{\nu}(\tau)=\frac{1}{2}~\int_{-1}^{+1}I_{\nu}(\tau,\mu)~ {\rm d}\mu$$ (2)

$$\begin{eqnarray*}\bar{J}_{ij}(\tau)=\int_{-\infty}^{+\infty}\phi_{\nu}(\tau)~ J_{\nu}(\tau)~ {\rm d}\nu \end{eqnarray*}$$

$$% \bar{J}_{ij}(\tau)=\frac{1}{2}~\int_{-\infty}^{+\infty}\phi... ...tau)~ \int_{-1}^{+1}I_{\nu}(\tau,\mu)~ {\rm d}\mu~{\rm d}\nu,$$ (3)

where $I_{\nu}(\tau,\mu)$ is the specific intensity at frequency $\nu$ in the direction $\mu = \cos (\theta )$, and $\phi_{\nu}(\tau)$ is the line profile. See Fig. 1 for the definitions, which are exactly the same as Goicoechea & Le Bourlot (2007) and reproduced here for convenience.

	Figure 1: Cloud geometry and sign conventions from Goicoechea & Le Bourlot (2007).
Open with DEXTER

Since the system of coupled detailed balance equations must be solved on a common grid, it is mandatory to compute the mean intensity for all lines on that same grid. Therefore, it is not possible to use the optical depth in the line as the independent variable as is usually done for a single line radiative-transfer problem. All relevant expressions have to be derived using a single common independent variable for all lines. In all that follows, we use the usual optical depth in the visible $\tau_{\rm V}$ as our independent variable ( $\tau_{\rm V} = 2.5 \log_{10}({\rm e})\; A_{\rm V}$).

  
2.1 Line profile

The line profile $\phi_\nu$ is often expressed using a non-dimensional centred variable $x = \frac{\nu - \nu_{\rm c}}{\Delta \nu_{\rm D}}$ such that $\phi_\nu~ {\rm d}\nu = \phi_x~ {\rm d}x$, with $\nu_{\rm c}$ the line centre frequency and $\Delta \nu_{\rm D}$ the Doppler width. However, in a medium with variable temperature and (potentially) turbulent velocity, this quantity is not constant as a function of the position in the cloud, and thus it is not possible to simply exchange integration over position in the cloud and integration over the line profile.

From the expression of $\Delta \nu_{\rm D} =\frac{\nu_{\rm c}}{c}~ \sqrt{\frac{2kT}{m}+v_{\rm t}^{2}}$, with T the kinetic temperature, m the molecular mass, and $v_{\rm t}$ the turbulent velocity, we can define a compound thermal velocity $v_{\rm T}(s)=\sqrt{\frac{2kT(s)}{m}+v_{\rm t}^{2}(s)}$ that accounts for all the (s) position dependent quantities. A natural variable is then $y = \frac{\nu - \nu_{\rm c}}{\nu_{\rm c}}$. The link between profiles is then $\phi_{y} =\nu_{\rm c}~\phi_{\nu}$.

For a Gaussian line profile, we have

$$% \phi_{y} = \frac{1}{\sqrt{\pi}}~ \frac{c}{v_{\rm T}(s)}~ \exp\left[-\left(y~\frac{c}{v_{\rm T}(s)}\right)^{2}\right]\cdot$$ (4)

This could be easily extended to take any velocity differences into account between a reference point (say, the cloud edge) and a position within the cloud, thus permitting also treatment of ``large velocity gradient'' cases, or any other macroscopic velocity distribution. Here, we consider a purely microturbulent contribution $v_{\rm t}$.

  
2.2 Radiative transfer equation

To compute the mean intensity (3), we need an expression for the specific intensity $I_{\nu}$. It follows from a formal solution to the 1D radiative transfer equation:

$$% \mu~ \frac{\partial I_{\nu}}{\partial s} = -\left(\kappa^{\... ...nu}+\eta_{ij}~\phi_{\nu}+\kappa^{\rm D}_{\nu}~ I^{\rm D}_{\nu}$$ (5)

where $\chi_{ji}=\frac{h\nu}{4\pi}~(B_{ji}~ n_{j}-B_{ij}~ n_{i})$and $\eta_{ij}=\frac{h\nu}{4\pi}~ A_{ij}~ n_{i}$, $\kappa^{\rm D}_{\nu}$ is the dust absorption coefficient, and $I^{\rm D}_{\nu}$ the dust emissivity. We assume complete redistribution, so that emission and absorption profiles are the same.

Note that in Eq. (5) the independent variable s is a ``physical'' distance. To introduce the optical depth  $\tau_{\rm V}$, we use the following relations:

$$% {\rm d}s = -C~\frac{1}{n_{\rm H}}~ {\rm d}\tau_{\rm V}$$ (6)

with $C = \frac{C_{\rm D}}{R_{\rm V}}~2.5~\log_{10}{\rm e}$, $C_{\rm D} = \frac{N_{\rm H}}{E_{B-V}}$, and $R_{\rm V}=\frac{A_{\rm V}}{E_{B-V}}$. With the standard values for our Galaxy, $C \simeq 2.0$ $\times$ $10^{21}~{\rm atoms~cm^{-2}}~{\rm mag^{-1}}$. The negative sign in the definition of $\tau_{\rm V}$ comes from the axis conventions, see Fig. 1.

We introduce the dust absorption cross section $\sigma^{\rm D}_{\nu}$, so $\kappa^{\rm D}_{\nu} = \sigma^{\rm D}_{\nu}\; n_{\rm H}$, with $n_{\rm H}$ the density of nucleons, $a^{\rm D}_{\nu}=C~\sigma^{\rm D}_{\nu}$, $D_{ij}=C~\frac{h\nu_{ij}}{4\pi}~ g_{i}~ A_{ij}$, $E_{ij} = \frac{c^{2}}{2h~\nu_{ij}^{3}}~ D_{ij}=C~ \frac{\lambda_{ij}^{2}}{8\pi}~ g_{i}~ A_{ij}$; and following Elitzur & Asensio Ramos (2006), we introduce $x_{i}=\frac{n_{i}}{g_{i}~ n_{\rm H}}$, where n_i is the number density of level i and g_i its statistical weight.

This leads to an expression of the radiative transfer equation where all physical constants and atomic and molecular data are clearly separated from physical variables describing local excitation conditions, i.e. depending explicitly on position within the cloud:

$$% \mu~\frac{\partial I_{\nu}}{\partial\tau_{\rm V}} = \left(a... ...u}-D_{ij}~ x_{i}~ \phi_{\nu}-a^{\rm D}_{\nu}~ I^{\rm D}_{\nu}.$$ (7)

  
2.3 Optical depths

From Eq. (7) we see that the total optical depth  $\tau_{\rm T}^{\nu}$ at frequency $\nu$ has two contributions (Dust and Line), $\tau_{\rm T}(\tau_{\rm V})=\tau_{\rm D}(\tau_{\rm V})+\tau_{\rm L}(\tau_{\rm V})$, with

$$% \tau_{\rm D}(\tau_{\rm V}) = \int_{0}^{\tau_{\rm V}}~ a^{\r... ... = C~\int_{0}^{\tau_{\rm V}}~\sigma^{\rm D}_{\nu}(t)~ {\rm d}t$$ (8)

$$% \tau_{\rm L}(\tau_{\rm V}) = E_{ij}~ \int_{0}^{\tau_{\rm V}} (x_{j}(t)-x_{i}(t))~\phi_{\nu}(t)~ {\rm d}t.$$ (9)

The first, $\tau_{\rm D}$, is independent of frequency, and reduces to $\tau_{\rm D}(\tau_{\rm V})=C~\sigma^{\rm D}_{\nu}~\tau_{\rm V}$ if the dust composition is constant throughout the cloud. The second, $\tau_{\rm L}$ depends on $\nu$ through $\phi_\nu$ only. At line centre and with Eq. (4), we get

$$\begin{eqnarray*}\tau_{\rm L}(\tau_{\rm V}) = \frac{E_{ij}~ c}{\nu_{\rm c}~\sqrt... ...au_{\rm V}}~ \frac{(x_{j}(t)-x_{i}(t))}{v_{\rm T}(t)}~ {\rm d}t. \end{eqnarray*}$$

This suggests introducing a new function of position  $X_{i}(\tau_{\rm V})$, which will be used in Sect. 3:

$$% X_{i}(\tau_{\rm V}) = \int_{0}^{\tau_{\rm V}}~ \frac{x_{i}(t)}{v_{\rm T}(t)}~ {\rm d}t.$$ (10)

  
2.4 Formal solution

For a cloud of total thickness $\tau_{\rm V}^{\rm max}$, the formal solution of Eq. (7) is

• $\mu<0$

  
  

  $$\begin{eqnarray*}&& I_{\nu}(\tau,\mu) = I_{\nu}(0,\mu)~ \exp\left( \frac{\tau_{\... ...m T}(t) -\tau_{\rm T}(\tau)}{\mu} \right)~ \frac{{\rm d}t}{\mu} \end{eqnarray*}$$
  
  

• $\mu>0$

  
  

  $$\begin{eqnarray*}&& I_{\nu}(\tau,\mu) = I_{\nu}(\tau_{\rm V}^{\rm max},\mu)~ \ex... ...(t)-\tau_{\rm T}(\tau)} {\mu} \right)~ \frac{{\rm d}t}{\mu}\cdot \end{eqnarray*}$$
  
  

The first case corresponds to rays propagating from left to right, with a boundary condition  $I_{\nu}(0,\mu)$ at $\tau_{\rm V} = 0$, and the second case corresponds to rays propagating from right to left with a boundary condition $I_{\nu}(\tau_{\rm V}^{\rm max},\mu)$ at $\tau_{\rm V} = \tau_{\rm V}^{\rm max}$.

We see that $I_{\nu}$ is the sum of three contributions. The first is the incident external radiation field attenuated by a factor  $\exp\left( \frac{\tau_{\rm T}}{\mu} \right)$. By splitting the integral in two, we get a second contribution from the line emission at t, $D_{ij}~ x_{i}~ \phi_{\nu}$ attenuated by the difference in total optical depth between the emission and absorption point, and a third contribution from the dust emission at t, $a^{\rm D}_{\nu}~ I^{\rm D}_{\nu}$ attenuated by the same factor. These two contributions are summed over the whole cloud. The arguments of the exponentials are always negative.

  
2.5 Mean intensity

The mean intensity over line profile is computed from Eqs. (2) and (3) and

$$% \bar{J}_{ij}(\tau_{\rm V}) = \frac{1}{2}~ \int_{-\infty}^{+... ..._{-1}^{+1}I_{\nu}~ (\tau_{\rm V},\mu)~ {\rm d}\mu~ {\rm d}\nu.$$ (11)

From the expression of $I_{\nu}$ we see that $\bar{J}_{ij}$ is the sum of three contributions that we consider now separately:

$$\begin{eqnarray*}\bar{J}_{ij}(\tau_{\rm V}) = \bar{J}^{\rm ext} + \bar{J}_{ij}^{\rm int} + \bar{J}_{\rm D}^{\rm int}. \end{eqnarray*}$$

  
2.5.1 External field contribution ${\bar{J}^{ext}}$

Let us consider an external radiation field that is isotropic on each side of the cloud  $I_{\rm ext}$ (although possibly different on each side) plus a contribution perpendicular to the cloud. Such a radiation field keeps the azimutal symmetry, while allowing for a point source (typically a close-by star) on the axis. We have

$$\begin{eqnarray*}&& I_{\nu}(0,\mu) = I_{\rm ext}^{-} + \delta_{-1\mu}~ I_{*}^{-}... ... V}^{\rm max},\mu) = I_{\rm ext}^{+} + \delta_{1\mu}~ I_{*}^{+}. \end{eqnarray*}$$

Inserting those expressions into Eq. (11) and assuming a constant radiation field over the line profile, we get

$$\begin{eqnarray*}\bar{J}^{\rm ext}(\tau_{\rm V}) &=& \frac{I_{*}^{-}}{2}~ \exp\l... ...}(\tau_{\rm V}^{\rm max})}{\mu} \right)~ {\rm d}\mu~ {\rm d}\nu. \end{eqnarray*}$$

In the following, we keep I_*^- = I_*⁺ = 0 for simplicity. With a change in variable $\mu \rightarrow \frac{1}{\alpha}$in the integrals, we recognise exponential integrals of the second kind, which leads to

$$% \bar{J}^{\rm ext}(\tau_{\rm V}) = I_{\rm ext}^{-}~ \beta_{-}(\tau_{\rm V}) + I_{\rm ext}^{+}~ \beta_{+}(\tau_{\rm V})$$ (12)

where $\beta_{-}$ and $\beta_{+}$ are the escape probabilities from position  $\tau_{\rm V}$ to the left or right side of the cloud, respectively:

$$\begin{eqnarray*}\beta_{-}(\tau_{\rm V}) = \frac{1}{2}~ \int_{-\infty}^{+\infty}... ...nu}~ E_{2}\left( \tau_{\rm T}(\tau_{\rm V}) \right)~ {\rm d}\nu \end{eqnarray*}$$

$$% \beta_{+}(\tau_{\rm V}) = \frac{1}{2}~ \int_{-\infty}^{+\in... ... V}^{\rm max})-\tau_{\rm T}(\tau_{\rm V}) \right)~ {\rm d}\nu.$$ (13)

These expressions are exact and show explicitly the symmetry between photons entering or leaving the cloud as long as they do not interact with the medium. Some properties of exponential integrals E₁ and E₂ and escape probabilities $\beta$ are described in Appendix A.1.

  
2.5.2 Dust contribution ${\bar{J}_{D}^{int}}$

We introduce a step function $\epsilon_{\tau_{\rm V}}(t)$ such that $\epsilon_{\tau_{\rm V}}(t)=1$ if $t<\tau_{\rm V}$ and $\epsilon_{\tau_{\rm V}}(t)=-1$if $t>\tau_{\rm V}$. With this definition, we can use a single function inside the exponential term, change the sign of the integration over $\mu$ for negative $\mu$, and write a single sum over t to get

$$\begin{eqnarray*}J_{\rm D}^{\rm int}(\tau_{\rm V})\!=\! \frac{1}{2} \int_{0}^{\t... ...\rm V})}{\mu} \right)~ \frac{{\rm d}\mu}{\mu} \right]~ {\rm d}t. \end{eqnarray*}$$

Setting $\mu \rightarrow \frac{1}{\alpha}$, we recognise in the inmost sum, the exponential integral of the first kind. Thus, after integration over the line profile, we get

$$% \bar{J}_{\rm D}^{\rm int}(\tau_{\rm V}) = \int_{0}^{\tau_{\... ...{\nu}(t)~ I^{\rm D}_{\nu}(t)~ L_{1,\tau_{\rm V}}(t)~ {\rm d}t,$$ (14)

where we have defined the kernel function $L_{1,\tau_{\rm V}}(t)$ by

$$% L_{1,\tau_{\rm V}}(t) = \frac{1}{2}~ \int_{-\infty}^{+\inft... ...au_{\rm T}(t)-\tau_{\rm T}(\tau_{\rm V}) \right))~ {\rm d}\nu.$$ (15)

Here, $L_{1,\tau_{\rm V}}(t)$ characterises the transport properties of photons from an emission point at t to the absorption point at  $\tau_{\rm V}$summed over all angles and all frequencies.

  
2.5.3 Line emission contributions ${\bar{J}_{ij}^{int}}$

Using a similar approach, we find

$$% \bar{J}_{ij}^{\rm int}(\tau_{\rm V}) = D_{ij}~\int_{0}^{\ta... ...i}(t)}{\Delta\nu_{\rm D}(t)}~ K_{1,\tau_{\rm V}}(t)~ {\rm d}t,$$ (16)

where we have defined a kernel function $K_{1,\tau_{\rm V}}(t)$ by

 

$$% K_{1,\tau_{\rm V}}(t) = \frac{\Delta\nu_{\rm D}(t)}{2}~ \int_{-... ... \left( \tau_{\rm T}(t)-\tau_{\rm T}(\tau_{\rm V}) \right) \right)~ {\rm d}\nu.$$ (17)

Note the existence of a second profile function $\phi$ in the definition of K₁, which comes from the fact that both emission and absorption are line processes. Furthermore, the two profile functions are not computed at the same spatial position and need not be identical. This is particularly important for a high kinetic temperature at emission and a low temperature at absorption; hence; the usefulness of the factorization of $\phi$ introduced in Sect. 2.1. The factor  $\Delta\nu_{\rm D}(t)$ is introduced so K₁ is a dimensionless quantity. Some properties of the kernel functions L₁ and K₁ are given in Appendix A.1.

These two kernel functions include all the necessary information on all radiative processes within the cloud. They are similar to, but different from, the functions L₁ and K₁ introduced by Mihalas (1978) in chapter 11, which provides the basis for this discussion. For a given $\tau_{\rm V}$, they are functions of position t and the unknown populations x_i through the contribution of $\tau_{\rm L}$ to $\tau_{\rm T}$.

  
2.6 Detailed balance

From Eqs. (12), (14), and (16) we see that $\bar{J}_{ij}$ may be written as

$$% \bar{J}_{ij} = \frac{2h\nu_{ij}^{3}}{c^{2}}~ Z_{ij}$$ (18)

where the term Z_ij is a function of position in the cloud ( $\tau_{\rm V}$) and the unknowns x_i and x_j over the whole cloud. This can be used to write a modified expression of Eq. (1):

 

$$x_{i}~ \left( \sum_{j<i}~ A_{ij}~ g_i~ \left( 1+Z_{ij} \right) + \sum_{j>i}~ A_{ji}~ g_{j}~ Z_{ji} \right)$$

$$\qquad + x_{i}~ \sum_{j\neq i}~ g_i~ k_{ij}^{X}n^{X} + x_{i}~D = \sum_{j\neq i}~ x_{j}~ g_j~ k_{ji}^{X}~ n^{X} +F_{i}$$

$$\qquad + \sum_{j<i}~ x_{j}~ A_{ij}~ g_{i}~ Z_{ij} + \sum_{j>i}~ x_{j}~ A_{ji}~ g_j~ \left( 1+Z_{ji} \right).$$ (19)

This is a set of n coupled non-linear integral equations (with n the number of levels of the considered species). Upon discretisation of the space variable  $\tau_{\rm V}$ at p positions ( $\tau_{\rm V}^{1} = 0$, $\tau_{\rm V}^{p} = \tau_{\rm V}^{\rm max}$), we have $p \times n$ ordinary equations in the $p \times n$ unknowns. Although cumbersome, the expressions for $\tau_{\rm L}$, L₁, K₁, and $\bar{J}_{ij}$are explicit enough to permit computation of all derivatives with respect to the variables x_i^l, where i indexes the level and l the position in the cloud. It also allows us to build a Jacobian matrix. Thus, the system may be solved in a fully consistent way in a single (large) step to provide an exact solution to the line transfer problem. This is in line with the method proposed in Elitzur & Asensio Ramos (2006), but with the further advantage that all physical quantities may be position dependent (density, temperature, turbulent velocity dispersion, dust emissivity, etc.). We have not yet implemented this computation into our code, but are planning to do so soon.

In the following, we focus on the effect of computing $\bar{J}_D^{\rm int}$ using the exact expression of Eq. (14) and then compare our results to PS05.

  
3 Approximate expressions

We first need to compute L₁ and K₁ explicitly. This requires a choice of line profile, and in the rest of this paper we use a Gaussian profile.

3.1 Nearly constant Gaussian wings

We need an explicit expression of the line optical depth  $\tau_{\rm L}$. From Eqs. (9) and (4) we get

$$\begin{eqnarray*}\tau_{\rm L}(\tau_{\rm V}) = \frac{E_{ij}~ c}{\nu_{\rm c}~ \sqr... ...2}~ \left( \frac{c}{v_{\rm T}(t)} \right)^{2} \right]~ {\rm d}t. \end{eqnarray*}$$

Introducing this expression into Eqs. (15) and (17) and introducing $z=y~\frac{c}{v_{\rm T}(\tau_{\rm V})}$, we get

$$\begin{eqnarray*}L_{1,\tau_{\rm V}}(t) = \frac{1}{\sqrt{\pi}}~ \int_{0}^{+\infty... ...\rm V}}(t)~ \Delta\tau_{\rm T}(t,\tau_{\rm V}) \right)~ {\rm d}z \end{eqnarray*}$$

and:

$$\begin{eqnarray*}K_{1,\tau_{\rm V}}(t) &=& \frac{1}{\pi}~ \int_{0}^{+\infty}~ {\... ...\rm V}}(t)~ \Delta\tau_{\rm T}(t,\tau_{\rm V}) \right)~ {\rm d}z \end{eqnarray*}$$

with:

$$\begin{eqnarray*}\Delta\tau_{\rm T}(t,\tau_{\rm V}) &=& \tau_{\rm D}(t)-\tau_{\r... ...(\tau_{\rm V})}{v_{\rm T}(t')} \right)^{2} \right]~ {\rm d}t'. \end{eqnarray*}$$

Those rather convolved expressions may be slightly simplified if $v_{\rm T}(\tau_{\rm V})$ and $v_{\rm T}(t')$ stay ``reasonably'' close to one another, which is the case for all species except hydrogen at low to moderate temperatures (say less than $1000~{\rm K})$. In this case, the ${\rm e}^{-z^2}$ term may be taken out of the inmost integral to give

$$% \Delta\tau_{\rm T}(t,\tau_{\rm V}) = \tau_{\rm D}(t)-\tau_{... ... \frac{\Delta X_{ij}(t,\tau_{\rm V})}{v_{\rm T}(\tau_{\rm V})}$$ (20)

with (using Eq. (10)):

$$\begin{eqnarray*}\Delta X_{ij}(t,\tau_{\rm V}) = X_j(t)-X_i(t)-X_j(\tau_{\rm V})+X_i(\tau_{\rm V}). \end{eqnarray*}$$

When $L_{1,\tau_{\rm V}}(t)$ is unchanged (formally), $K_{1,\tau_{\rm V}}(t)$becomes

$$\begin{eqnarray*}K_{1,\tau_{\rm V}}(t) = \frac{1}{\pi}~ \int_{0}^{+\infty}~ {\rm... ...rm V}}(t)~ \Delta\tau_{\rm T}(t,\tau_{\rm V}) \right)~ {\rm d}z. \end{eqnarray*}$$

Using a Simpson quadrature rule for $\Delta X_{ij}(t,\tau_{\rm V})$, this is a ``simple'' analytical function of x_i. Note the crucial difference between L₁ and K₁: the factor  ${\rm e}^{-2z^2}$instead of ${\rm e}^{-z^2}$. This single factor of 2 leads to huge differences in the behaviour of the two kernels (see Fig. A.4).

3.2 An approximation for ${\bar{J}_{ij}^{int}(\tau_{V})}$

Using how K₁ decreases as $\vert t-\tau_{\rm V}\vert$ increases, we may evaluate the integration over t in Eq. (17) using for all functions of t their value at $\tau_{\rm V}$. First of all, the total optical depth  $\tau_{\rm T}$ between tand $\tau_{\rm V}$ becomes

$$\begin{eqnarray*}\tau_{\rm T}(t)-\tau_{\rm T}(\tau_{\rm V}) \simeq \Lambda_{ij}~ (t-\tau_{\rm V}) \end{eqnarray*}$$

with

$$\begin{eqnarray*}\Lambda_{ij} = a^{\rm D}_{\nu} + E_{ij}~ (x_{j}(\tau_{\rm V})-x_{i}(\tau_{\rm V}))~ \phi_{\nu}(\tau_{\rm V}). \end{eqnarray*}$$

By splitting the integration over t at $\tau_{\rm V}$ (and using $\epsilon_{\tau_{\rm V}} = \pm 1$ as required), using the definition of E₁ and changing the integration order over $\mu$ and t, we get

$$\begin{eqnarray*}\bar{J}_{ij}^{\rm int}(\tau_{\rm V}) = \frac{D_{ij}~ x_{i}}{2} ... ...\tau_{\rm V})\right){\rm d}t~ \frac{{\rm d}\mu}{\mu}~ {\rm d}\nu \end{eqnarray*}$$

where x_i and $\phi_\nu$ are evaluated at $\tau_{\rm V}$. The inmost integration can be carried analytically, leading to

$$\begin{eqnarray*}\bar{J}_{ij}^{\rm int}(\tau_{\rm V}) = \frac{D_{ij}~ x_{i}}{2} ... ...}^{\rm max}-\tau_{\rm V})\right)\right]~ {\rm d}\mu~ {\rm d}\nu. \end{eqnarray*}$$

We do a second approximation by assuming that, in the line, $\Lambda_{ij}$ is dominated by the line contribution. We may thus neglect the factor  $a^{\rm D}_{\nu}$, and simplify one factor of $\phi_\nu$ in the integral. This gives precisely Eq. (13) for the two double integrals, and we finally get

$$\begin{eqnarray*}\bar{J}_{ij}^{\rm int}(\tau_{\rm V}) = \frac{2h~ \nu_{ij}^{3}}{... ...ft(1-\beta_{\rm L}(\tau_{\rm V})-\beta_{R}(\tau_{\rm V})\right), \end{eqnarray*}$$

where one recovers the usual expression of the source function at $\tau_{\rm V}$. One sees that to get that approximation, we first need to consider that the source function is constant at its value at $\tau_{\rm V}$ everywhere in the cloud, and that the dust contribution to opacity over the whole line profile is negligible. The first assumption is crude, but that the kernel function K₁ peaks at $\tau_{\rm V}$ helps a lot. The second assumption may be quite crude, too, which people are not always aware of.

3.3 An approximation for ${\bar{J}_{D}^{int}(\tau_{V})}$

PS05 (Eqs. (3) and (4)) use

$$% \bar{J}_{\rm D}^{\rm int}(\tau_{\rm V}) = \beta~ \tau^{\rm D}~ I^{\rm D},$$ (21)

where $\beta$ is the escape probability and $\tau^{\rm D}$ the total dust optical depth. Appendix B shows how to recover that approximation from the full expression of $\bar{J}_{\rm D}^{\rm int}$ (Eqs. (14) and (15))

However, in the following sections, we use that full expression, using only the approximate expression of $\Delta\tau_{\rm T}$ of Eq. (20) to simplify the computation of L₁. The integration over L₁itself is performed using a trapezoidal rule, except in the vicinity of $\tau_{\rm V}$ where L₁ has a logarithmic divergence. There we use an analytic integration of an expansion for the short argument.

  
4 Line processes

  
4.1 Line formation

Line profiles can be easily computed from a knowledge of the structure of the cloud - as computed by the code - after the whole iteration has converged. The process is two-stepped: first compute the optical depth at a specific frequency from the ``far'' side of the slab along a ray up to the front. This can be done with Eqs. (8) and (9). Then the formal solution to the transfer equation is computed by using a simple trapezoidal rule that includes three contributions: attenuated radiation from the background (which may reduce to the CMB, but may also include some illuminating source), emission from the dust along the whole path, and emission from the molecules themselves.

The intensity in the line is the contribution from that single source, integrated over frequency. It is important to note that potential emission from matter (mainly dust behind the zone in the cloud where the line forms) is screened by absorption in the line (followed either by collisional deexcitation or emission in $4\pi$ steradian) and will thus not contribute to the underlying background. Therefore, it is not always possible to compute an integrated line intensity from observations by simply removing a base line extending beyond the line. One must have at least a fair knowledge of the continuum emissivity to compute the full line emissivity.

This effect can lead to underestimating the line intensity by at least a factor of two for the most intense lines of H_2O. Figure 2 shows the contributions to the total emissivity of the 2₁₂ to 1₀₁ line in the model of Sect. 5.1 of dust (upper panel) and molecules (lower panel). The dust contribution is high in a tiny slab at the cloud edge where the dust is warm. Deeper into the cloud, its emissivity is lower but remains non vanishing. However, those photons are blocked at line centre (where, in this case, the line optical depth is 1.13) as can be seen from the black area for frequencies close to line centre (up to about two Doppler width). Line emissivity is high on the front side of that high optical depth zone. The decreasing emissivity in the wing is partially offset by a longer optical path, leading to the wing features seen in Fig. 3.

	Figure 2: Local emissivity of the 202 to 111 transition of H_2O as a function of position in the line (in units of Doppler width) and location on the slab. Upper panel, dust contribution, lower panel line contribution.
Open with DEXTER

	Figure 3: H_2O $2_{02} \rightarrow 1_{11}$ line profile from model PDR-c of Sect. 5. Note the explicitly computed contribution from the dust.
Open with DEXTER

It can be seen that the continuum level is lower under the line than far into the wings, so that a simple linear baseline drawn from the observed continuum misses half of the line signal. Only the strong lines with great optical depth are subject to this bias. However, these are the most likely detectable lines by forthcoming facilities such as the Herschel HIFI instrument. This is shown quantitatively in Fig. 4.

	Figure 4: Ratio of real to apparent line emissivity of the lowest lying lines of H_2O as a function of line frequency for model PDR-c. See Sect. 4.1 for details.
Open with DEXTER

  
4.2 Line cooling

It is customary to compute the cooling contribution of any species by a sum of the energy released by all radiative transitions from an upper level i to a lower level j. However, this approximation is correct only in the limit where excitation is dominated by collisions. This is not the case as soon as the lines become thick, and is particularly erroneous if radiative pumping by infrared dust emission leads to significant excitation. A simple sum would include in the cooling the energy at a given position brought there by a photon originating from somewhere else in the cloud, which never had a chance to release its energy in a kinetic energy form.

The only way to accurately take into account the net energy balance of the gas is to consider the gains and losses directly through collisions. This can be illustrated by the example of a two-level atom. Consider the balance including induced transitions, but without chemical formation and destruction terms:

$$\begin{eqnarray*}x_1 \left( k^{X}_{12}~ n^{X} + B_{12} \bar{J}_{12}~ \right) = x... ...eft( k^{X}_{21}~ n^{X} + A_{21} + B_{21}~ \bar{J}_{12} \right). \end{eqnarray*}$$

The net cooling through this transition is the difference between kinetic energy removed from the gas, proportional to x₁  k₁₂  nminus what is injected into the gas x₂  k₂₁  n (per particle). From the detailed balance equation, we have

$$\begin{eqnarray*}\frac{x_2}{x_1} &=& \frac{g_2}{g_1} \exp\left( -\frac{T_{12}}{T... ...right)}{ \exp \left( \frac{T_{12}}{T_{\rm e}}\right)-1} \right), \end{eqnarray*}$$

where T_ij is the energy of the transition in Kelvin, $T_{\rm k}$ the kinetic temperature, and the radiation field is expressed as an equivalent black body at a temperature of $T_{\rm e}$. The factor $1+ \frac{ A_{21}}{ k_{21}~ n}$has been introduced to recover the classical expression of a two-level atom with spontaneous emission, leading to the definition of the critical density $n_{\rm c}$. For a vanishing radiative pumping ( $T_{\rm e} \rightarrow 0$), one recovers the usual two-level atom.

For transitions well below that critical density, the correction factor (second line of the above expression) simplifies to,

$$\begin{eqnarray*}1-\exp\left( -\frac{T_{12}}{T_{\rm e}}\right) + \frac{A_{21}}{k... ...\left( \frac{T_{12}}{T_{\rm k}}-\frac{T_{12}}{T_{\rm e}}\right) \end{eqnarray*}$$

which is always larger more 1. Thus, radiative pumping by infrared radiation tends to decrease cooling by enhancing collisional deexcitation.

  
5 Results and discussion

  
5.1 Comparison with S140

In this section, we compare the formalism proposed here to the results of PS05 on S140, using the same model parameters to the extent possible. We use four different models, all based on the Meudon PDR code as described in Le Petit et al. (2006), which differ only by the expression used for $J_{\rm int}^{\rm D}$. PDR-a does not take IR pumping into account, PDR-b1, and PDR-b2 use Eq. (21) as used in PS05, and PDR-c uses the formalism of Sect. 2. PDR-b1 and PDR-b2 differ by the expression of the total dust optical depth used in Eq. (21):

$$\begin{eqnarray*}&& \tau^{\rm D}_1 = 10^{-3}\; \frac{100~\mu {\rm m}}{\lambda} \\ && \tau^{\rm D}_2 = N_{{\rm H}}\; \kappa_{\rm abs}, \end{eqnarray*}$$

where $\tau^{\rm D}_1$ is the expression used in PS05, while $\tau^{\rm D}_2$uses a fit to the dust cross section per nucleon  $\kappa_{\rm abs}$ from the standard galactic model of Weingartner & Draine (2001) as updated by Draine (2003a,b) and available from http://www.astro.princeton.edu/draine/dust/dustmix.html. For an energy T in Kelvin in the range $[3~{\rm K}, 400~{\rm K}]$:

$$\begin{eqnarray*}\log_{10}(\kappa_{\rm abs}) = 0.194601\; \log_{10}^2(T) + 1.3307\; \log_{10}(T) -28.0681 \end{eqnarray*}$$

where $N_{{\rm H}}$ is the slab column density. The ratio of the two approximations $\tau^{\rm D}_2 / \tau^{\rm D}_1$ is displayed in Fig. 5 and the expressions of  $J_{\rm int}^{\rm D}$ summarised in Table 1.

	Figure 5: Ratio of dust opacities used in models PDR-b2 over PDR-b1.
Open with DEXTER

Table 1: Infrared pumping approximations.

That comparison allows us to understand the effect of IR pumping on physical processes computed in the model. The main input parameters are displayed in Table 2. The chemical network involves about 100 species linked by 1500 reactions. Gas phase elemental abundances are given in Table 3.

Table 2: Physical parameters used to model the PDR of the molecular cloud S140.

Table 3: Gas phase abundances relative to H nuclei.

  
5.1.1 Physical properties of the cloud

Figure 6 shows the gas temperature profile resulting from thermal balance. It ranges from $\sim$ $480~{\rm K}$ at the edge of the cloud, closest to the star, and decreases to $\sim$ $10~{\rm K}$ at cloud centre. The dust temperature ranges from $\sim$ $50~{\rm K}$ to $\sim$ $25~{\rm K}$. These temperatures are the same in all models, showing that, in this case, the overall structure is not affected by cooling in infrared lines.

Our temperature range is within an order of magnitude of PS05. Roellig et al. (2007) have shown that the gas temperature is a very sensitive output of PDR codes. Thus, given the many differences in the computations, our agreement is satisfactory.

	Figure 6: Gas and dust temperature in our S140 models (log scale). The star is to the left of the figure. Higher dust temperature at the edge leads to infrared pumping deeper into the cloud.
Open with DEXTER

	Figure 7: The H_2O abundance in our S140 models (PDR-abc).
Open with DEXTER

	Figure 8: Ortho to para ratio of H_2O in our S140 models. The ratio is set by the balance of chemical formation and destruction alone, with formation proportional to the Boltzman factor.
Open with DEXTER

Figure 7 shows the H_2O abundance. At low temperatures (< $300~{\rm K}$), water is formed by a chain of ion-neutral reactions initiated by cosmic ray ionization. Closer to the edges, water is destroyed by photodissociation. PS05 give a maximum relative abundance of H_2O of a few 10^-6, quite close to our result.

  
5.1.2 Excitation of water molecules

A detailed balance of H_2O is computed using Eq. (19). Collision rates are taken from Green et al. (1993) for collisions with He and from Phillips et al. (1996) for collisions with H_2, as found on http://basecol.obspm.fr/. The explicit inclusion of chemical formation and destruction terms (D and F_i) leads to a set of equations that does not require the use of the conservation condition $\sum_i~x_{i}=1$. Thus, conservation can be checked a posteriori as a numerical test. We also computed explicitly the ortho to para ratio (O/P) of H_2O as a function of position, by taking all processes into account leading to the formation of one or the other form. Here, the only term governing that ratio is the balance of chemistry through a temperature-dependent formation branching ratio F_i. To our knowledge, there are no reactive collision rates with H or H_2 that have been published and could be used. The resulting O/P ratio is displayed in Fig. 8. The lowest para-level 0₀₀ and ortho-level 1₀₁account for most of H_2O and are almost not affected by the different pumping schemes.

This is not true for higher lying levels, as can be seen in Figs. 9 and 10 which show the populations of levels 2₀₂ (para) and 2₁₂ (ortho). Model PDR-c, with its higher pumping rate, populates upper levels more effectively. Neglecting infrared pumping completely (model PDR-a) may lead to huge differences and will not be considered further. The difference between models PDR-b1 and PDR-b2 shows the importance of using the best possible grain model, although this is still a very uncertain field. The same dust properties are used in models PDR-b2 and PDR-c, and it can be seen that the resulting populations are in fair agreement. However, exact computation of the pumping rate results in central abundances about five times higher, which translates in the line emissivities. The discrepancy increases with level energy as can be seen in Figs. 11.

	Figure 9: Relative abundance of para level 202 of H_2O for models PDR-abc (log scale).
Open with DEXTER

	Figure 10: Relative abundance of ortho level 212 of H_2O for models PDR-abc (log scale).
Open with DEXTER

	Figure 11: Relative abundance of ortho level 220 of H_2O for models PDR-ab (log scale).
Open with DEXTER

	Figure 12: Relative level populations of the first six ortho levels of H_2O (log scale). Model PDR-c: exact infrared pumping computation.
Open with DEXTER

	Figure 13: Relative level populations of the first five para levels of H_2O (log scale). Model PDR-c: exact infrared pumping computation.
Open with DEXTER

Figures 12 and 13 show the lowest ortho and para levels, respectively, for model PDR-c, and Table 4 gives the average mean intensities of the first few ortho and para transitions. Agreement is fair with PS05 for the strongest transitions.

We checked that CO is not affected by infrared pumping (variations are less than a few percent), as could be expected from its low permanent dipole moment. In the other hand CS, with a much higher dipole, may be a better target; however, transitions with the right energy are much too high to involve populated levels significantly, so that they are barely modified.

Table 4: Mean intensities of the first rotational transitions of ortho and para  H_2O (in ${\rm erg~cm^{-2}~s^{-1}~sr^{-1}}$) and (last column) optical depth at line centre for model PDR-c.

  
5.1.3 Cooling in S140

Cooling is not affected much here, because the main coolants are insensitive to infrared pumping, as could be guessed from the identical temperature profiles (see Fig. 14). As explained in Sect. 4.2, cooling by H_2O (which accounts for less than $1 \%$ of the total) is reduced by $15 \%$ when infrared pumping is taken into account.

The reduced water cooling comes from radiative energy transfer from the warm outer parts of the cloud to the central water-rich core through trapping of the infrared radiation. The higher populated excited levels transfer their energy by collisional deexcitation to the gas, which reduces the amount of kinetic energy removed by collisional excitation.

	Figure 14: Model PDR-c, cooling rates (log scale). Cooling in the molecular core comes mainly from CO and its isotopes.
Open with DEXTER

We compared the local cooling rate of water to the empirical proxi found in Pineau des Forêts et al. (1986). We found that, where this term matters (high abundance of H_2O), the formula is within a factor of two from our results, which is a tribute to the physical insight of those authors.

We have also simulated a denser cloud where H_2O is expected to contribute more to cooling. We choose a symmetric radiation field of $\chi =100$ on both sides, and a lower value for $A_{\rm V}=10$. Figure 15 is similar to Fig. 14 and shows that H_2O now contributes about $\sim$$20\%$ to the total cooling at cloud centre.

	Figure 15: Cooling rates (log scale) for a cloud illuminated by a field of $\chi =100$ at both sides, with $A_{\rm V}=10$ and $n_{\rm H}=2$ $\times$ $10^{6}~{\rm cm}^{-3}$, exact infrared pumping. H_2O cannot be neglected any more.
Open with DEXTER

  
6 Conclusions

This paper presents a careful analysis of the cooling radiation within a 1D interstellar cloud model containing gas and dust. The various contributions involved in the excitation of atomic and molecular species are discussed, including possible radiative pumping due to internal infrared sources such as dust emission. The main result is that infrared emission coming from hot grains at the edge of the cloud may penetrate deep inside the cloud, providing an efficient radiation source at a location where the cold grains can no longer emit infrared radiation. This procedure is necessary for abundant light species that have their fundamental transitions in the tera-Hertz region, such as  H_2O.

We specifically address the emission of H_2O submillimetric transitions, but the present formalism may be extended to other molecules such as HCN, HNC, floppy carbon chains that present the appropriate level structure, and the like. We applied our model to the S 140 PDR region, which may be taken as a benchmark case, and compared our results with the approach of PS05. We recovered the approximations used by PS05 when we introduced the analytical approximation of the internal mean intensity due to dust, expressed in the escape formalism, and make the additional assumption: $\beta \sim \frac{1}{1+\tau}$.

We have discussed the various modifications in the line intensities with various levels of approximations. Including the exact treatment of the infrared pumping consistently leads to an increase in the line emissivity values, and factors up to two orders of magnitude may result. As shown by PS05 inhomogeneities also impact the integrated line intensities by large factors. It is probable that both effects should be taken into account.

The cooling efficiency does not follow the line emissivities directly. Indeed, adding the infrared pumping contribution implies a higher excitation and also an additional heating contribution. The high excited levels may also be de-excited by collisions, in addition to radiative emission, which reduces the cooling efficiency. The present paper assumes a constant grain size distribution over the whole cloud. This assumption is a shortcoming of the present model and will be removed in subsequent studies.

Despite these inherent model simplifications, including the geometrical assumptions, we think that future observations with Herschel will benefit from this significant improvement in the theoretical treatment of line emissivities. Extension to arbitrary geometry of our semi-analytical computation is not possible since several key equation transformations depend on the 1D assumption. However, other simple geometries are probably tractable. For instance, a spherical cloud with a density profile that depends only on the radial coordinate, but with illumination coming both from an isotropic interstellar radiation field and a close-by star is a good candidate. We expect that the equations are cumbersome, have made no attempt to derive them.

This geometry would lead to smaller optical depths to some cloud boundaries from the bulk of the gas, thus probably further enhancing the line intensities. However, exposure to external ultraviolet photons would also be enhanced, reducing the abundance of molecular species. Thus the overall effect is not easy to anticipate.

Acknowledgements

Part of this work was supported by PCMI national programme of the CNRS. Manuel Gonzalez Garcia acknowledges a grant from the European Marie Curie program ``The Molecular Universe''.

References

• Bakes, E. L. O., & Tielens, A. G. G. M. 1994, ApJ, 427, 822 [NASA ADS] [CrossRef]
• Bernes, C. 1979, A&A, 73, 67 [NASA ADS]
• Chevalier, L., Paletou, F., & Rutily, B. 2003, A&A, 411, 221 [NASA ADS] [CrossRef] [EDP Sciences]
• Draine, B. T. 2003a, ARA&A, 41, 241 [NASA ADS] [CrossRef]
• Draine, B. T. 2003b, ApJ, 598, 1017 [NASA ADS] [CrossRef]
• Elitzur, M., & Asensio Ramos, A. 2006, MNRAS, 365, 779 [NASA ADS] [CrossRef] (In the text)
• Goicoechea, J. R., & Le Bourlot, J. 2007, A&A, 467, 1 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Green, S., Maluendes, S., & McLean, A. D. 1993, ApJS, 85, 181G [NASA ADS] [CrossRef] (In the text)
• Hubeny, I. 2001, ASPCS, 247, 197 [NASA ADS]
• Juvela, M. 2005, A&A, 440, 531 [NASA ADS] [CrossRef] [EDP Sciences]
• Le Petit, F., Nehmé, C., Le Bourlot, J., & Roueff, E. 2006, ApJS, 164, 506 [NASA ADS] [CrossRef] (In the text)
• Mathis, J., Rumpl, W., & Nordsieck, K. 1977, ApJ, 217, 425 [NASA ADS] [CrossRef]
• Mihalas, D. 1978, Stellar Atmospheres (Freeman and Co) (In the text)
• Phillips, T. R., Maluendes, S., & Green, S. 1996, ApJSS, 107, 467 [NASA ADS] [CrossRef] (In the text)
• Pineau des Forêts, G., Flower, D. R., Hartquist, T. W., & Dalgarno, A. 1986, MNRAS, 220, 801 [NASA ADS] (In the text)
• Poelman, D. R., & Spaans, M. 2005, A&A, 440, 559 [NASA ADS] [CrossRef] [EDP Sciences]
• Poelman, D. R., & Spaans, M. 2006, A&A, 453, 615 [NASA ADS] [CrossRef] [EDP Sciences]
• Poelman, H., Spaans, M., & Tielens, A. G. G. M. 2007, A&A, 464, 1023 [NASA ADS] [CrossRef] [EDP Sciences]
• Röllig, M., Abel, N. P., Bell, T., et al. 2007, A&A, 467, 187 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• van der Tak, F. F. S., Walmsley, C. M., Herpin, F., & Ceccarelli, C. 2006, A&A, 447, 1011 [NASA ADS] [CrossRef] [EDP Sciences]
• van der Tak, F. F. S., Black, J. H., Schoier, F., Janssen, D., & van Dishoeck, E. 2007, A&A, 468, 627 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• van Zadelhoff, G.-J., Dullemond, C. P., van der Tak, F. F. S., et al. 2002, A&A, 395, 373 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296 [NASA ADS] [CrossRef] (In the text)

 

  
Online Material

  
Appendix A: Properties of various functions

  
A.1 Properties of exponential integrals

Exponential integrals are defined by

$$\begin{eqnarray*}E_{n}(x) = \int_{1}^{\infty}~ \frac{{\rm e}^{-xt}}{t^{n}}\;{\rm d}t. \end{eqnarray*}$$

They verify the relation

$$\begin{eqnarray*}\frac{{\rm d}E_{n}}{{\rm d}x} = -E_{n-1}. \end{eqnarray*}$$

Both E₁ and E₂ are displayed in Fig. A.1. One can see that E₂(0)=1 whether $E_{1}(0) \rightarrow \infty$. Close to 0E₁ may be expanded as

$$\begin{eqnarray*}E_{1}(x) \simeq x -\gamma -\ln(x) \end{eqnarray*}$$

where $\gamma$ is the Euler constant.

	Figure A.1: Integral exponentials E1 and E2.
Open with DEXTER

  
A.2 Properties of the escape probability function $\beta$

The total escape probability $\beta(\tau)$ from a slab of gas of total optical depth at line centre  $\tau_{\rm max}$ is

$$\begin{eqnarray*}\beta(\tau)= \beta_{+}(\tau) + \beta_{-}(\tau) \end{eqnarray*}$$

where $\beta_{+}$ and $\beta_{-}$ are computed from Eq. (13). We recall that $\beta_{-}(\tau)$ corresponds to an optical depth $\tau$ towards one side of the cloud, and $\beta_{+}(\tau)$ to an optical depth  $\tau_{\rm max}-\tau$towards the other side.

	Figure A.2: Escape probability $\beta$. $\tau _{\rm m}$ is the optical depth at cloud centre, and $\tau$ the optical depth towards the closest edge.
Open with DEXTER

The evolution of $\beta(\tau)$ as a function of the total optical depth  $\tau_{\rm max}$and position into the cloud $\tau$ is displayed in Fig. A.2. Three plateaus are visible in that plot: for optically thin lines (small ( $\tau_{\rm max}$), the escape probability $\beta$ is close to 1 everywhere. Then, as the line becomes thicker, $\beta$decreases, and two cases exists: either we are close to one edge, and the escape probability is close to 0.5, or we are deep into the cloud, and $\beta$is close to 0. Transitions from one to the other are smooth, but not trivial.

	Figure A.3: Escape probability from cloud centre $\beta\left(\frac{\tau_{\rm max}}{2}\right)$.
Open with DEXTER

The escape probability from cloud centre as a function of the cloud size is shown in Fig. A.3. Overplotted is also the function $f(\tau) = \frac{1}{1+x}$. One can see that it is a fair approximation of $\beta$. However, this is only true at cloud centre! Closer to the edges for thick lines, $\beta$ goes to 0.5 when the approximation goes to 1.

  
A.3 Properties of functions L₁ and K₁

	Figure A.4: Kernel functions $L_1(\tau )$ and $K_1(\tau )$. The approximate expansion for small $\tau$ is displayed up to $\tau =2$ and is excellent up to $\tau \simeq 0.3$.
Open with DEXTER

Figure A.4 shows the kernel functions $L_1(\tau )$ and $K_1(\tau )$defined by (see Sect. 3)

$$\begin{eqnarray*}L_{1}(\tau) = \frac{1}{\sqrt{\pi}}~ \int_{0}^{+\infty}~ {\rm e}... ...eft( \frac{{\rm e}^{-z^{2}}}{\sqrt{\pi}}~ \tau \right)~ {\rm d}z \end{eqnarray*}$$

and:

$$\begin{eqnarray*}K_{1}(\tau) = \frac{1}{\pi}~ \int_{0}^{+\infty}~ {\rm e}^{-2~z^... ...ft( \frac{{\rm e}^{-z^{2}}}{\sqrt{\pi}}~ \tau \right)~ {\rm d}z. \end{eqnarray*}$$

Integration is performed using a Gauss-Hermite scheme. There K₁ decreases much faster than L₁ for large $\tau$, which shows that coupling through continuum extends farther than the simple line coupling. Thus infrared radiation from a warm region may affect significantly colder, more deeply shielded molecular regions of a cloud. For small $\tau$, we may use the above expansion for E₁to get

$$\begin{eqnarray*}L_{1}(\tau) = \frac{1}{2}~ \left[ \frac{1}{2}-\gamma +\frac{ \t... ...sqrt{2\pi}} - \ln\left( \frac{\tau}{\sqrt{\pi}}\right) \right] \end{eqnarray*}$$

and

$$\begin{eqnarray*}K_{1}(\tau) = \frac{1}{2\sqrt{2\pi}} \left[ \frac{1}{4} -\gamma... ...\pi}}~\tau -\ln\left( \frac{ \tau}{\sqrt{\pi}} \right) \right], \end{eqnarray*}$$

which explicitly shows the mild logarithmic divergence at $\tau=0$.

  
Appendix B: Analytical approximation to ${\bar{J}_{D}^{int}}$

The approximate value of $\bar{J}_{\rm D}^{\rm int}$ used in PS05 may be derived from Eqs. (14) and (15). First, all quantities dependent on t are taken at their value at $\tau_{\rm V}$and taken out of the integrals. The sum over t is split at $\tau_{\rm V}$, and integration over $\nu$ and t are swapped to give

 

$$\bar{J}_{\rm D}^{\rm int}(\tau_{\rm V}) = \frac{1}{2}a^{\rm D}_{\... ...}} \!\! E_{1}(\tau_{\rm T}(t)\!-\!\tau_{\rm T}(\tau_{\rm V}))~{\rm d}t \right].$$ (B.1)

From Eqs. (8) and (9), we see that $\left( \tau_{\rm T}(t) -\tau_{\rm T}(\tau_{\rm V}) \right)$ may be written formally as $\alpha~(t-\tau_{\rm V})$, which leads to

 

$$% \bar{J}_{\rm D}^{\rm int}(\tau_{\rm V}) = \frac{1}{2}a^{\rm D}_... ...m V}}^{\tau_{\rm V}^{\rm max}}~ E_{1}(\alpha(t-\tau_{\rm V}))~{\rm d}t \right].$$ (B.2)

We make the change of variable $x=-\alpha~(t-\tau_{\rm V})$ in the first integral and $x=\alpha~(t-\tau_{\rm V})$ in the second to get

 

$$% \bar{J}_{\rm D}^{\rm int}(\tau_{\rm V}) = \frac{a^{\rm D}_{\nu}... ...t_{0}^{\alpha,(\tau_{\rm V}^{\rm max}-\tau_{\rm V})}~ E_{1}(x)~{\rm d}x \right]$$ (B.3)

Now using $E_{1}(x)=-\frac{{\rm d}}{{\rm d}x}E_{2}(x)$ we have

 

$$% \bar{J}_{\rm D}^{\rm int}(\tau_{\rm V}) = \frac{a^{\rm D}_{\nu}... ...lpha~\tau_{\rm V}) - E_{2}(\alpha~(\tau_{\rm V}^{\rm max}-\tau_{\rm V})\right).$$ (B.4)

And we recognise (cf. Eq. (13)), the expression of the escape probabilities. Now multiplying and dividing by $\tau_{\rm V}^{\rm max}$ and writing $a^{\rm D}~\tau_{\rm V}^{\rm max} = \tau^{\rm D}_{\rm max}$, the total dust optical depth, we have

$$% \bar{J}_{\rm D}^{\rm int}(\tau_{\rm V}) = \tau^{\rm D}_{\rm... ...{+}(\tau_{\rm V})}{\tau_{\rm T}(\tau_{\rm V}^{\rm max})}\cdot$$ (B.5)

The final expression comes from the approximation $\beta \sim \frac{1}{1+\tau}$ discussed above and leads to the expression used in PS05:

$$% \bar{J}_{\rm D}^{\rm int}(\tau_{\rm V}) = \tau^{\rm D}_{\rm max}~I^{\rm D}_{\nu}\; \beta(\tau_{\rm V}).$$ (B.6)