Appendix B: SimLine - A one-dimensional radiative transfer code

B.1 The radiative transfer problem

SimLine solves the line radiative transfer problem discussed in Appendix A.1 in a spherically symmetric configuration by means of a $\lambda$-iteration. The code is similar to the concept described by Dickel & Auer (1994) but it contains several extensions and achieves a higher accuracy from an adaptive discretisation of all independent quantities.

SimLine integrates the radiative transfer Eq. (A.2) for a number of rays numerically. In spherical symmetry it is sufficient to consider the propagation of radiation in one arbitrary direction which is taken as z here. The integral is computed stepwise from z_i-1 to z_i

 

$$I_\nu(p,z_i) = \exp \left( -\int_{z_{i-1}}^{z_i} {\rm d}z\; \kapp... ...nu(p,z) \exp \left( \int_{z_{i-1}}^z {\rm d}z' \kappa_\nu(p,z') \right) \right]$$ (B.1)

where p denotes the displacement variable perpendicular to the z direction. To minimise the discretisation error the integral does not use the source function but only the emission and absorption coefficients which are linear in the level populations. They are assumed to change linearly between the grid points and the exact integration formula for a linear behaviour (which is not given here, but can be derived straight forward) is applied. This approach provides a reduction of the integration error to third order. The radial grid is dynamically adjusted to give a maximum variation of the level populations between two neighbouring points below a certain limit. In case of strong velocity gradients additional points are included on the z-scale for a sufficiently dense sampling of changes in the profile function. The incident radiation at the outer boundary of the cloud is assumed to follow a black body spectrum.

In spherical symmetry the spatial integration of the radiative energy density (Eq. (B.12)) can be reduced to

 

$${2\pi\over r c} \int_{-r}^r {\rm d}z' \int_{-\infty}^\infty {\rm d}\nu\; \Phi_\nu(p',z') I_\nu(p',z')$$ u(r) = (B.2)

$${\rm with} p'=\sqrt{r^2-z'^2}.$$

The grid of rays tangential to the radial grid is refined by additional rays at intermediate p' values to guarantee a sufficiently dense sampling on the z' scale. The integration uses a cubic spline interpolation.

With the values of the radiative energy density at each radial point and for each transition, the system of balance equations can be solved providing new level populations. Here, a LU decomposition algorithm with iterative improvement (Press et al. 1992) is used. The new level populations are used in the next iteration as input for the radiative transfer equation. The whole $\lambda$-iteration scheme is solved using the convergence accelerator introduced by Auer (1987).

Depending on the physical situation the initial guess is either the optically thin limit, thermalisation or the solution of the radiative transfer equation using the LVG approximation (Eq. (B.15)). The stability of the local radiation field is used as convergence criterion. The number of iterations required for convergence depends strongly on the optical depth of the model cloud. For the examples discussed in this paper only about a dozen iterations were necessary but other test cases with complex molecules like water, non-monotonic velocity gradients, and high optical depths require several hundred iterations.

B.2 The local turbulence approximation

The turbulence description uses two additional parameters for each spatial point: the width of the velocity distribution $\sigma$ providing the local emission profile for optically thin lines and the correlation length of the macroturbulent density or velocity distribution $l_{\rm corr}$.

The width of the velocity distribution $\sigma$ is composed of a turbulent and a thermal contribution

$$\sigma = {\nu_0 \over c} \,\sqrt{ {2kT_{\rm kin} \over m_{\rm mol}} + {2 \over 3} \langle v_{\rm turb}^2 \rangle}$$ (B.3)

where a Maxwellian distribution of turbulent velocities is assumed. The relation between the FWHM and the variance of the turbulent velocity distribution is given by FWHM $(v_{\rm turb})=\sqrt{8/3\times{\rm ln}2\; \langle v_{\rm turb}^2 \rangle}$). The long range variation of the turbulence spectrum as described by means of a Kolmogorov or Larson exponent is simulated by a radially varying turbulent velocity dispersion $\sqrt{\langle v_{\rm turb}^2\rangle} \propto r^\gamma$. Exponents $\gamma$ between about 0.1 (Goodman et al. 1998) and 0.7 (Fuller & Myers 1992) are observationally justified.

For the local treatment of coherent units in a turbulent medium the considered volume element is subdivided into numerous clumps with a thermal internal velocity dispersion. When each clump is characterized by a Gaussian density distribution of molecules with about the same velocity $n(r)=n_0\times\exp(-r^2/r_{\rm cl}^2)$the effective absorption coefficient at the considered velocity for the whole medium is

$$\kappa_{\rm eff}=n_{\rm cl} \times \pi r_{\rm cl}^2 \int_0^{\tau_{\rm cl}} {1-\exp(-\tau) \over \tau} {\rm d}\tau$$ (B.4)

where $n_{\rm cl}$ is the number density of contributing cells and $\tau_{\rm cl}=\sqrt{\pi}\,\kappa r_{\rm cl}$ is their central opacity (Martin et al. 1984). As the clumps size $r_{\rm cl}$ is the length on which the abundance of molecules within the same thermal velocity profile is reduced by the factor $1/\rm e$, we can compute it from the correlation length of the velocity or density structure by $r_{\rm cl}=l_{\rm corr}\times \sigma_{\rm th} /\sigma$.

When the turbulent velocity dispersion $\sigma$ is at least three times as large as the thermal velocity dispersion $\sigma_{\rm th}$, we obtain an effective absorption coefficient

$$\kappa_{\rm eff,\nu}=n_{\rm ges} \pi r_{\rm cl}^2 \times A(\t... ...ver \sigma} \exp{}\left(-{(\nu-\nu_0)^2 \over \sigma^2}\right)$$ (B.5)

with

$$A(\tau)={1\over \sqrt{\pi}}\int_{-\infty}^{\infty}{\rm d}v \int_0^{\tau \exp(-v^2)} {1-\exp(-\tau') \over \tau'} {\rm d}\tau'.$$ (B.6)

Here, $n_{\rm ges}$ is the total number density of clumps. In case of incompressible turbulence, i.e. clumping in velocity space, it is equal to the reciprocal cell volume. For small values of the clump opacity, $A(\tau_{\rm cl})$ is identical to $\tau_{\rm cl}$ and we reproduce the microturbulent limit. For large $\tau_{\rm cl}$, the function $A(\tau_{\rm cl})$ saturates and we obtain a significant reduction of the effective absorption coefficient. In case of density clumps, $\kappa_{\rm eff}(\nu)$ is further reduced by the filling factor entering $n_{\rm ges}$. In SimLine, this is simulated by a corresponding artificial reduction of the molecular abundance.

B.3 The central H II region

To simulate the effect of a central continuum source in the cloud, it is possible to assume an H II region in the cloud centre. The H II region is characterised by two parameters, the electron density $n_{\rm e}$ and the kinetic electron temperature $T_{\rm e}$.

The absorption coefficient for electron-ion bremsstrahlung in the Rayleigh-Jeans approximation is given by:

$$\kappa_\nu={8 \over 3\sqrt{2\pi}}\; {{e}^6 \over (4\pi\epsilo... ... \left(m_{\rm e}\over k T_{\rm e}\right)^{3/2} {\rm ln}\Lambda$$ (B.7)

where it is assumed that the gas is singly ionised and $\Lambda$ is given by

$$\Lambda=\left(2k T_{\rm e} \over \gamma m_{\rm e}\right)^{3/2... ...mes10^7 \left(T\over {\rm K}\right)^{3/2} {{\rm Hz} \over \nu}$$ (B.8)

for $T_{\rm e}<3.6\times10^5$ K. The quantities e and $m_{\rm e}$ denote the electron charge and mass, c is the vacuum light velocity, and $\gamma=1.781$ (Lang 1980).

For a thermal plasma, the emission coefficient follows from the Planck function

$$\epsilon_\nu=\kappa_\nu\times B_\nu(T_{\rm e})\cdot$$ (B.9)

In the radiative transfer computations the frequency dependence of these continuum coefficients is neglected within the molecular line width. Within the H II region, we substitute the molecular coefficients in Eq. (B.1) by the quantities from Eqs. (B.7) and (B.9) so that we locally switch to a continuum radiative transfer.  

B.4 Computation of beam temperatures

When the level populations are known, the beam temperature relative to the background is computed by the convolution of the emergent intensity with the telescope beam.

$$T_{\rm mb}={c^2 \over 2 k \nu_0^2} \;{\displaystyle \int_0^{2... ...\rm d}\phi \int_0^\infty {\rm d}p\; p f_{\rm mb}(p,\phi)}\cdot$$ (B.10)

The intensity $I_\nu^{\rm S}(p)$ is the value on the cloud surface $I_\nu^{\rm S}(p)=I_\nu(p,\sqrt{R_{\rm cloud}^2-p^2})$. We assume a Gaussian profile for the telescope beam

$$f_{\rm mb}(p,\phi)= \exp{}\left(-(p-p_{\rm offset})^2(1+\phi^2)^2 \over \sigma_{\rm mb}^2\right)\cdot$$ (B.11)

The projected beam width is computed from the angular width by $\sigma_{\rm mb} = \pi/648\,000\, D\, \sigma_{\rm mb}[''] = 2.912\times10^{-6}\, D$ FWHM[''] where D is the distance of the cloud. The program computes a radial map with arbitrary spacings.

B.5 The general code design

The design of the code is directed towards a high accuracy of the computed line profiles. All errors in the different steps of the program are explicitly user controlled by setting thresholds. All discretisations necessary to treat the problem numerically are performed in an adaptive way, i.e. there is no predefined grid and all grid parameters will change during the iteration procedure. The system of balance equations is truncated whenever the excitation of all higher levels falls below a chosen accuracy limit.

Furthermore, the code was pushed towards a high flexibility, i.e. the ability to treat a very broad range of physical parameters with the same accuracy and without numerical limitations. The systematic velocities e.g. may range from 0 to several times the turbulent velocity and the optical depths may vary from negative values for weak masing to values of several thousands.

The program is not optimised towards a high speed. Other codes with lower inherent accuracy may easily run a factor 10 faster, and further improvements are possible. Nevertheless, the code is suitable for an interactive work even on a small PC with execution times of a few seconds for the models considered in this paper.