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 -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  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

 (1)

where ni is the population of level i, Aij and Bijthe Einstein coefficients, kjiX the collision rate with species X, nX the abundance of collision partner X, D the total destruction rate, and Fi the formation rate on level i. Both D and Fi 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  , 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  within a cloud, is computed from

 (2)

 (3)

where is the specific intensity at frequency  in the direction , and 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 as our independent variable ( ).

#### 2.1 Line profile

The line profile is often expressed using a non-dimensional centred variable such that , with the line centre frequency and 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 , with T the kinetic temperature, m the molecular mass, and the turbulent velocity, we can define a compound thermal velocity that accounts for all the (s) position dependent quantities. A natural variable is then . The link between profiles is then .

For a Gaussian line profile, we have

 (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 .

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

 (5)

where and , is the dust absorption coefficient, and 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  , we use the following relations:

 (6)

with , , and . With the standard values for our Galaxy,   . The negative sign in the definition of comes from the axis conventions, see Fig. 1.

We introduce the dust absorption cross section , so , with the density of nucleons, , , ; and following Elitzur & Asensio Ramos (2006), we introduce , where ni is the number density of level i and gi 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:

 (7)

#### 2.3 Optical depths

From Eq. (7) we see that the total optical depth  at frequency  has two contributions (Dust and Line), , with

 (8)

 (9)

The first, , is independent of frequency, and reduces to if the dust composition is constant throughout the cloud. The second, depends on through only. At line centre and with Eq. (4), we get

This suggests introducing a new function of position  , which will be used in Sect. 3:

 (10)

#### 2.4 Formal solution

For a cloud of total thickness , the formal solution of Eq. (7) is

The first case corresponds to rays propagating from left to right, with a boundary condition  at , and the second case corresponds to rays propagating from right to left with a boundary condition at .

We see that is the sum of three contributions. The first is the incident external radiation field attenuated by a factor  . By splitting the integral in two, we get a second contribution from the line emission at t, attenuated by the difference in total optical depth between the emission and absorption point, and a third contribution from the dust emission at t, 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

 (11)

From the expression of we see that is the sum of three contributions that we consider now separately:

##### 2.5.1 External field contribution

Let us consider an external radiation field that is isotropic on each side of the cloud  (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

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

In the following, we keep I*- = I*+ = 0 for simplicity. With a change in variable in the integrals, we recognise exponential integrals of the second kind, which leads to

 (12)

where and are the escape probabilities from position  to the left or right side of the cloud, respectively:

 (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 E1 and E2 and escape probabilities  are described in Appendix A.1.

##### 2.5.2 Dust contribution

We introduce a step function such that if and if . With this definition, we can use a single function inside the exponential term, change the sign of the integration over for negative , and write a single sum over t to get

Setting , we recognise in the inmost sum, the exponential integral of the first kind. Thus, after integration over the line profile, we get

 (14)

where we have defined the kernel function  by

 (15)

Here, characterises the transport properties of photons from an emission point at t to the absorption point at  summed over all angles and all frequencies.

##### 2.5.3 Line emission contributions

Using a similar approach, we find

 (16)

where we have defined a kernel function by

 (17)

Note the existence of a second profile function  in the definition of K1, 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 introduced in Sect. 2.1. The factor  is introduced so K1 is a dimensionless quantity. Some properties of the kernel functions L1 and K1 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 L1 and K1 introduced by Mihalas (1978) in chapter 11, which provides the basis for this discussion. For a given , they are functions of position t and the unknown populations xi through the contribution of to .

#### 2.6 Detailed balance

From Eqs. (12), (14), and (16) we see that may be written as

 (18)

where the term Zij is a function of position in the cloud ( ) and the unknowns xi and xj over the whole cloud. This can be used to write a modified expression of Eq. (1):

 (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  at p positions ( , ), we have  ordinary equations in the unknowns. Although cumbersome, the expressions for , L1, K1, and are explicit enough to permit computation of all derivatives with respect to the variables xil, 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 using the exact expression of Eq. (14) and then compare our results to PS05.

### 3 Approximate expressions

We first need to compute L1 and K1 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  . From Eqs. (9) and (4) we get

Introducing this expression into Eqs. (15) and (17) and introducing , we get

and:

with:

Those rather convolved expressions may be slightly simplified if and stay reasonably'' close to one another, which is the case for all species except hydrogen at low to moderate temperatures (say less than . In this case, the  term may be taken out of the inmost integral to give

 (20)

with (using Eq. (10)):

When is unchanged (formally), becomes

Using a Simpson quadrature rule for , this is a simple'' analytical function of xi. Note the crucial difference between L1 and K1: the factor  instead of . 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

Using how K1 decreases as increases, we may evaluate the integration over t in Eq. (17) using for all functions of t their value at . First of all, the total optical depth  between tand becomes

with

By splitting the integration over t at (and using as required), using the definition of E1 and changing the integration order over and t, we get

where xi and are evaluated at . The inmost integration can be carried analytically, leading to

We do a second approximation by assuming that, in the line, is dominated by the line contribution. We may thus neglect the factor  , and simplify one factor of in the integral. This gives precisely Eq. (13) for the two double integrals, and we finally get

where one recovers the usual expression of the source function at . One sees that to get that approximation, we first need to consider that the source function is constant at its value at 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 K1 peaks at helps a lot. The second assumption may be quite crude, too, which people are not always aware of.

#### 3.3 An approximation for

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

 (21)

where is the escape probability and the total dust optical depth. Appendix B shows how to recover that approximation from the full expression of (Eqs. (14) and (15))

However, in the following sections, we use that full expression, using only the approximate expression of of Eq. (20) to simplify the computation of L1. The integration over L1itself is performed using a trapezoidal rule, except in the vicinity of where L1 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  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 212 to 101 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 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:

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

where Tij is the energy of the transition in Kelvin, the kinetic temperature, and the radiation field is expressed as an equivalent black body at a temperature of . The factor has been introduced to recover the classical expression of a two-level atom with spontaneous emission, leading to the definition of the critical density . For a vanishing radiative pumping ( ), 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,

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 . 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):

where is the expression used in PS05, while uses a fit to the dust cross section per nucleon  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 :

where is the slab column density. The ratio of the two approximations is displayed in Fig. 5 and the expressions of  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 at the edge of the cloud, closest to the star, and decreases to at cloud centre. The dust temperature ranges from to . 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 (< ), 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 Fi) leads to a set of equations that does not require the use of the conservation condition . 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 Fi. 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 000 and ortho-level 101account 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 202 (para) and 212 (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 ) 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 of the total) is reduced by 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 on both sides, and a lower value for . Figure 15 is similar to Fig. 14 and shows that H_2O now contributes about to the total cooling at cloud centre.

 Figure 15: Cooling rates (log scale) for a cloud illuminated by a field of at both sides, with and   , 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: .

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

### Appendix A: Properties of various functions

#### A.1 Properties of exponential integrals

Exponential integrals are defined by

They verify the relation

Both E1 and E2 are displayed in Fig. A.1. One can see that E2(0)=1 whether . Close to 0E1 may be expanded as

where is the Euler constant.

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

#### A.2 Properties of the escape probability function

The total escape probability from a slab of gas of total optical depth at line centre   is

where and are computed from Eq. (13). We recall that corresponds to an optical depth  towards one side of the cloud, and to an optical depth  towards the other side.

 Figure A.2: Escape probability . is the optical depth at cloud centre, and the optical depth towards the closest edge. Open with DEXTER

The evolution of as a function of the total optical depth  and position into the cloud  is displayed in Fig. A.2. Three plateaus are visible in that plot: for optically thin lines (small ( ), the escape probability  is close to 1 everywhere. Then, as the line becomes thicker, 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 is close to 0. Transitions from one to the other are smooth, but not trivial.

 Figure A.3: Escape probability from cloud centre . 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 . One can see that it is a fair approximation of . However, this is only true at cloud centre! Closer to the edges for thick lines, goes to 0.5 when the approximation goes to 1.

#### A.3 Properties of functions L1 and K1

 Figure A.4: Kernel functions and . The approximate expansion for small is displayed up to and is excellent up to . Open with DEXTER

Figure A.4 shows the kernel functions  and defined by (see Sect. 3)

and:

Integration is performed using a Gauss-Hermite scheme. There K1 decreases much faster than L1 for large , 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 , we may use the above expansion for E1to get

and

which explicitly shows the mild logarithmic divergence at .

### Appendix B: Analytical approximation to

The approximate value of used in PS05 may be derived from Eqs. (14) and (15). First, all quantities dependent on t are taken at their value at and taken out of the integrals. The sum over t is split at , and integration over and t are swapped to give

 (B.1)

From Eqs. (8) and (9), we see that may be written formally as , which leads to

 (B.2)

We make the change of variable in the first integral and in the second to get

 (B.3)

Now using we have

 (B.4)

And we recognise (cf. Eq. (13)), the expression of the escape probabilities. Now multiplying and dividing by and writing , the total dust optical depth, we have

 (B.5)

The final expression comes from the approximation discussed above and leads to the expression used in PS05:

 (B.6)