LTR - Radiative transfer in molecular lines

Version 1.40 (6/17/1996)

Dr. V. Ossenkopf, University Observatory Jena

Users Manual

Abstract

LTR is a FORTRAN code to compute the profiles of molecular rotational lines for linear molecules in spherically symmetric clouds with arbitrary density, temperature and velocity structure. Turbulence and clumping effects are treated in a local statistical approximation combined with a radial dependence of the correlation parameters. The code consists of two parts: the self-consistent solution of the balance equations for all level populations and energy densities at all radial points and the computation of the emergent line profiles observed from a telescope with finite beam width and arbitrary offset. All numerical discretizations are done in an adaptive way, the user has full error control and the line profiles may be computed in almost arbitrary accuracy. An accelerated lambda iteration is used to solve the equation system. The optical depths in the lines may vary between 0 and about 3000. The input of physical cloud parameters, telescope parameters, and numerical control parameters may be done either interactively or as file input.

1 The basic theory

1.1 The solution of the radiative transfer problem

The general formulation of the frequency dependent radiative transfer equation

\frac{dI(ν, r⃗, n⃗)}{ds} = −κ(ν, r⃗, n⃗)I(ν, r⃗, n⃗) + ϵ(ν, r⃗, n⃗)  (1)

may be reduced in the spherically symmetric case down to three independent coordinates. It is sufficient to consider radiation in the direction of the z axis and to neglect the ϕ dependence of the radius vector. Then, the radiative transfer equation can be written as

\frac{dI_{z}(ν, p, z)}{dz} = −κ(ν, p, z)I_{z}(ν, p, z) + ϵ(ν, p, z)  (2)

where p denotes the displacement variable (p = \sqrt{x^{2} + y^{2}}). The z axis is chosen in the direction towards the observer. Since the spatial grid is adjusted at distinct radial points, the program does not use the form of the radiative transfer equation with source function and derivative on the optical depth, but with emission and absorption coefficients.

For the numerical solution of the radiative transfer equation it is convenient to use stepwise the solution of the integral equation allowing a reduction of the error to third order

I_{z}(ν, p, z_{i}) = exp (− ∫_{z_{i−1}}^{z_{i}}κ(ν, p, z)dz) [I_{z}(ν, p, z_{i−1}) + ∫_{z_{i−1}}^{z_{i}}ϵ(ν, p, z) exp (∫_{z_{i−1}}^{z} κ(ν, p, z′)dz′) dz] (3)

As incident radiation at the outer boundary of the cloud, the cosmic background radiation

I_{bg}(ν) = \frac{2hν^{3}}{c^{2}} [exp (\frac{hν}{k × 2.7K}) − 1]^{−1}  (4)

will be taken.

The program assumes complete redistribution of energy between the molecules at a given level, so that one can define a unique set of level populations n_{j} at any point r⃗ which is independent from frequency or the actual velocity of a molecule. Furthermore, this means that the line profiles for emission and absorption are the same. In molecular clouds it is sufficient to consider pure Doppler broadening so that the profiles are given by

Φ(ν) = \frac{1}{\sqrt{π}σ} exp (−\frac{(ν − ν_{0}(1 + v_{z}/c))^{2}}{σ^{2}} )  (5)

where v_{z} is the systematic velocity in z direction, σ is the local line width, and ν is the frequency in the observers frame. (Without loss of generality, the program assumes that the centre of the cloud is in rest within the observers frame.).

With the assumtion of complete redistribution, the absorption and emission coefficients for a given transition are determined by

κ_{j}(ν) = \frac{hν_{0}}{c} (n_{j}B_{j,j+1} − n_{j+1}B_{j+1,j})Φ(ν)  (6)
ϵ_{j}(ν) = \frac{hν_{0}}{4π} n_{j+1}A_{j+1,j}Φ(ν)  (7)

where stimulated emission is treated as negative absorption.

The level populations at a given radius r may be computed from the linear system of balance equations

n_{j}(A_{j,j−1}+B_{j,j−1}u_{j}+B_{j,j+1}u_{j+1}+∑_{i≠j} c_{j,i}) = n_{j+1}(A_{j+1,j}+B_{j+1,j}u_{j+1})+n_{j−1}B_{j−1,j}uj+∑_{i≠j} n_{i}c_{i,j} (8)

where the radiative energy density u_{j} for a transition j → j − 1 at a given point is computed from the integral

u_{j}(p, z) = \frac{1}{c} ∫_{0}^{2π} dϕ ∫_{0}^{π} dθ sinθ ∫_{−∞}^{∞} dν Φ_{ϕ,θ}(ν, v⃗)I(ν, p, z, ϕ, θ) .  (9)

Exploiting the spherical symmetry, this may be transformed into

u_{j}(r) = \frac{2π}{rc} ∫_{−r}^{r} dz′ ∫_{−∞}^{∞} Φ(ν)I_{z}(ν, p′ = \sqrt{r^{2} − z^{′2}}, z′)dν  (10)

where the radius r is given as r = \sqrt{p^{2} + z^{2}}. Due to the symmetry, only radiation parallel to the z axis is encountered. For the central line frequency ν_{0} within the expression for Φ (Eq. 5), the frequency of the transition j → j − 1 has to be taken when computing u_{j}.

The equation system is truncated whenever the excitation of a level falls below the chosen accuracy limit or at the maximum level for which the collision rates are known.

1.2 The local turbulence approximation

The turbulence description uses two parameter for each point r⃗: the total width of the velocity distribution σ giving the local profile for optically thin lines and the correlation length r_{corr}.

The width of the velocity distribution σ is composed of a turbulent and a thermal contribution

σ = \frac{ν_{0}}{c} \sqrt{\frac{2kT_{kin}}{m} + \frac{2}{3}⟨v_{turb}^{2}⟩}  (11)

where a Maxwellian distribution of turbulent velocities is assumed. In the data input, the turbulent velocities are given as FWHM of the velocity distribution (FWHM(v_{turb}) = \sqrt{8/3 × ln2 ⟨v_{turb}^{2}⟩}). The thermal velocities are automatically computed from the kinetic gas temperature.

The long range correlation of the turbulence spectrum as described by a Kolmogorov exponent can be simulated in 1D by a radial dependence of the turbulent velocity dispersion ⟨v_{turb}^{2}⟩ ∝ r^{γ}. Observations of Fuller & Myers (1992) suggest exponents around γ ≈ 0.7. Several other authors proposed somewhat smaller values around γ ≈ 0.5.

For the local treatment of clumping in real space or velocity space the considered volume element is subdivided into numerous clumps whith a thermal internal velocity dispersion. From a Gaussian distribution of the density of molecules with a certain velocity within the single fragments

n(r) = n_{0} × exp(−r^{2}/r_{cl}^{2})  (12)

we get an effective absorption coefficient for the whole medium given by

κ_{eff} = n_{cl} × πr_{cl}^{2} ∫_{0}^{τ_{cl}} \frac{1 − exp(−τ)}{τ}dτ  (13)

(Martin et al. 1984), where n_{cl} is the number density of clumps which contribute at the considered frequency and τ_{cl} = \sqrt{π} κr_{cl} is their central opacity.

Assuming a Maxwellian turbulent velocity distribution of clumps and a thermal velocity dispersion σ_{th} giving at most 1/3 of the total velocity dispersion σ, we obtain an effective absorption coefficient

κ_{eff}(ν) = n_{ges}πr_{cl}^{2} × A(τ_{cl}) × \frac{σ_{th}}{σ} exp (−\frac{(ν − ν_{0})^{2}}{σ^{2}} )  (14)

with

A(τ) = \frac{1}{\sqrt{π}} ∫_{−∞}^{∞} dv ∫_{0}^{τ exp(−v^{2})} \frac{1 − exp(−τ′)}{τ′}dτ′  (15)

Here, n_{ges} is the total number of clumps. In case of pure turbulence in velocity space, n_{ges} is given by the volume of the single clumps and we obtain

κ_{eff}(ν) = κ(ν) × A(τ_{cl})/τ_{cl} .  (16)

Consequently, A(τ_{cl})/τ_{cl} is a measure of the reduction of the opacity due to turbulence. For small clump sizes, A(τ_{cl}) = τ_{cl} so that we are in the microturbulent limit. For τ_{cl} ≫ 1 the macroturbulent limit is reached. In case of real clumping in space, κ_{eff}(ν) is further reduced by the filling factor. In the program, this is simulated by a corresponding artificial reduction of the molecular abundance.

The clumps size r_{cl} is determined by the correlation length of the velocity or density structure. r_{cl} is the length on which the abundance of molecules within the same thermal velocity profile is reduced by the factor 1/e. Consequently, it can be computed from the turbulence correlation length by r_{cl} = r_{corr} × σ_{th}/σ.

One has to keep in mind that this is a statistical approximation, always assuming a complete Maxwellian velocity distribution within each volume element even if the single clumps become relatively large. Furthermore, the source function S = ϵ(ν)/κ(ν) is not influenced in this local approach.

1.3 The central H II region

For the simulation of the internal heating produced by a central continuum source in the cloud, it is possible to assume an H II region in the cloud core. The H II region is characterized by two parameters, the electron density N_{e} and the kinetic electron temperature T_{e}.

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

κ(ν) = \frac{8}{3\sqrt{2π}} \frac{e^{6}}{(4πϵ_{0}m_{e})^{3}c} (\frac{N_{e}}{ν} )^{2} (\frac{m_{e}}{kT_{e}} )^{3/2} lnΛ  (17)

where it is assumed that the gas is singly ionized and Λ is given by

Λ = (\frac{2kT_{e}}{δm_{e}} )^{3/2} \frac{4πϵ_{0}m}{πδe^{2}ν} ≈ 4.9573 10^{7} (\frac{T}{K})^{3/2} \frac{Hz}{ν}  (18)

for T_{e} < 3.2 10^{5}K. The quantities e and m_{e} denote the electron charge and mass and c is the vacuum light velocity.

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

ϵ(ν) = κ(ν) × B_{ν}(T_{e}) .  (19)

In the radiative transfer computations within the H II region the small frequency dependence of these continuum coefficients within the molecular rotational lines is neglected.

1.4 The Sobolev approximation

As an initial guess to the level populations, the Sobolev approximation may be used. This approximation is strictly valid for clouds with large velocity gradients and velocities which are increasing outwards so that all radiative coupling is restricted to a small spatial region. Then the level populations can be computed from local quantities only.

The radiative energy density at a radius r is determined by the local source function S_{j}(r) = ϵ_{j}(r)/κ_{j}(r)

u_{j}(r) = (1 − β_{j}(r)) S_{j}(r) + β_{j}(r)I_{bg}(ν_{0})  (20)

where β_{j}(r) denotes the probability for the escape of a photon within the line j from the cloud when emitted at r. It can be computed from an integral over all spatial directions from a given point by

β_{j}(r) = \frac{1}{2} ∫_{−1}^{1} \frac{1 − exp(−κ_{j}(r)/Q(r, μ))}{κ_{j}(r)/Q(r, μ)} dμ  (21)

with the line integrated effective absorption coefficient κ_{j} and the velocity gradient in spherical symmetry

Q(r, μ) = μ^{2} \frac{dv_{r}}{dr} + (1 − μ^{2})\frac{v_{r}}{r}  (22)

Hence, the radial velocity gradient is the most important quantity for the local energy density.

In contrast to the iterative solution of Eq. (8) it is possible to apply a more efficient way here, since the energy densities are analytically coupled to the level populations via Eq. (7) and (20-22). Using the derivatives of the matrix coefficients, the balance equations can be solved by a Newton-Raphson approach. The corrections to the level populations within each step will be computed from the matrix equation

∑_{i} (∑_{k} \frac{∂A_{jk}}{∂n_{i}} n_{k} + A_{ji}) Δn_{i} = ∑_{i} A_{ji}n_{i}  (23)

where the A_{ij} are the coefficients of the matrix of balance equations (Eq. 8).

1.5 Computation of beam temperatures

When the level populations are known, the beam temperature relative to the background is computed from the emergent intensity by integration over the projection of the telescope beam on the cloud.

T_{beam}(ν) = \frac{c^{2}}{2kν_{0}^{2}} \frac{∫_{0}^{2π} dϕ ∫_{0}^{∞} dp p × (I_{surf}(ν, p) − I_{bg}(ν))f_{beam}(p, ϕ)}{∫_{0}^{2π} dϕ ∫_{0}^{∞} dp p × f_{beam}(p, ϕ)}  (24)

The emergent intensity is the value on the cloud surface I_{surf}(ν, p) = I_{z}(ν, p, \sqrt{R_{cloud}^{2} − p^{2}}). A Gaussian profile for the beam is assumed

f_{beam}(p, ϕ) = \frac{1}{\sqrt{π}σ_{beam}} exp (\frac{−(p − p_{offset})^{2}(1 + ϕ^{2})^{2}}{σ_{beam}^{2}} )  (25)

The projected beam width is computed from the angular width by σ_{beam} = πD^{2}/648000 σ_{beam}[′′] where D is the distance of the cloud. The standard deviation σ_{beam} is coupled to the full width of half power by FWHM = 2\sqrt{ln2} σ_{beam}. The program computes a full map, i.e. the line profiles at a given number of positions on a linear radial scan through the cloud.

Within the program, all frequencies are treated in units of (ν − ν_{0})/ν_{0}, velocities in units of c, and intensities in units of ν_{0}^{3}/c^{2}, i.e. by a factor 4π/h larger than the SI unit.

2 The general code design

The general design of the program is directed towards a high accuracy of the computed line profiles. All of the errors in the different steps of the program are explicitely user controlled by setting thresholds. All discretizations necessary to treat the problem numerically are performed in an adaptive way, i.e. there is no predefined grid at all and all grid parameters may change during the iteration procedure.

Furthermore, the code was pushed towards a high flexibility, i.e. the ability to treat a very broad range of all physical parameters with the same accuracy and without numerical limitations. E.g., the systematic velocities may range from 0 to about 20 times the turbulent velocity which may be further extended by increasing the size of the fields for the frequency points. There is no inherent restriction to a distinct range of optical depths. I have tested it for depths at line centre between 0 and about 3000. However, in case of very high optical depths or large velocity gradients the code will run remarkably slower.

The program is not optimized towards a high velocity of the code. Although I got a factor of about 60 in velocity relative to the precursor program by E. Krügel, other codes with lower inherent accuracy may run another factor 10 faster. Nevertheless, the code is suitable for an interactive work even on my PC-486.

3 User interface

3.1 Programm invocation

The radiative transfer code is started from the command line by:

ltr [options] where options can be:

 -fileinput
 -numdefaults
 -overwrite
 -block or -lines
 -nocomments, -unixcomments or -vmscomments

When the option -fileinput is given, the program assumes that all physical, observational and numerical parameters will be read from files; interactive parameter input is switched off. The option -numdefaults sticks the set of the used numerical parameters to the default values. The other options set the behaviour for the data file output (see Sect. 3.3).

All the options are intended to avoid additional interactive questions during the program execution. If options are omitted, the corresponding parameters will be determined interactively during the program execution.

3.2 Terminal interaction

The program is designed for an interactive computation of line profiles so that the flow is controlled by terminal input. However, all parameters necessary for the computation may be given either interactively or may be read from data files.

There is a strict separation between parameters describing the molecular cloud, parameters of the observation, and numerical parameters which control the accuracy of the different parts of the computation.

A parameter input which was done interactively may be stored immediately in a file with the appropriate structure. This file can be read at the next time of use, and modifications of single values can be done by means of a text editor. In case of errors made when interactively defining the data, the input procedure can be repeated by typing a character which cannot be converted to a numerical value at the next input question.

The two parts of the program may be run in arbitrary order and repeatedly. At first one can compute either a new cloud model, i.e. the level populations at all radial points throughout the cloud, or one can read a file which contains the results of a previous computation. Then, it is possible to compute the observed line profiles for an arbitrary number of lines or telescope parameters. For a hand-made line fit one can change some cloud parameters and repeat the whole procedure.

I recommend to use a graphics software in a parallel session/window to get an immediate impression from the changes on the line profiles.

3.3 Data files

All files treated by the program are plain ASCII files so that they can be inspected by eye, modified with a text editor, and used for further reduction in other programs.

The input files providing the physical, observational and numerical parameters consist of sets of dual lines. The first line in each pair is treated as a comment line by the program and, in general, contains a description of the parameter. The second line contains the parameter in ASCII form. I highly recommend to produce a first set of input files using the LTR program itself by interactive input of the parameters. This guarantees the consistency of the produced parameter files with the structure expected when reading the parameters. Starting from this first set it may be more feasible to change distinct values in these files using a text editor.

Two types of result files are produced. The first one contains the level populations calculated as the result of the lambda iteration between populations and radiative energy densities. There are no comments in the file. The first line contains the molecule description, the number of radial points and the number of rotational levels treated. In the following block, the physical parameters at the different radial points are stored. Each line contains the number of the radial point, its radius, the molecule density, the systematic velocity, the thermal+turbulent velocity σ, and the correlation length over σ . A second block contains the lines with the level populations. For each radial point (denoted by its number) the populations at all computed levels are stored. A third block is created for testing purposes only and will not be used when reading the data for a new computation. It gives the energy densities at each radial point and for all transitions.

The same file structure as produced here (except for the third block) must be given when reading a precomputed cloud model from a file in order to calculate new line profiles.

The second file type contains the line profiles. These files are designed to be readable by graphics programs as GREG or GNUPLOT. Two general file structures are possible - compact and line oriented. In the first one, the two dimensional field of beam temperatures is stored as a block where the lines represent the different beam offsets and the columns the frequency/velocity points. In two preceding comment lines the offset and frequency values are given. The line-oriented structure requires much more disk space but can be read by more programs. Here, each line contains three numbers - the beam offset, the frequency, and the beam temperature. Consequently the file contains (number of offsets)×(number of frequency points) lines. Preceeding comment lines specify these numbers. The format of comment lines can be chosen to be in VMS or UNIX format, i.e. whether a ‘!’ or ‘#’ shall be used to indicate the beginning of a comment line. Furthermore one can also skip the comment lines.

The default parameters for the structure of the line profile files and the possibility to always overwrite existing files with newer ones can be specified as options on the command line. The options -block and -lines turn on block mode or line mode for the file structure. The options -nocomments, -unixcomments, and -vmscomments define the format of the comment lines. When the additional option -overwrite is given, existing files will be replaced without further notice. Otherwise, every overwrite has to be confirmed interactively. All these options can be given also on compile time via conditional defines (e.g. -Dnocomments -Doverwrite). If no parameters are specified, they will be asked interactively before each file output.

The units in the files are always pc for radii and km/s for velocities. Beam offsets are given in arcseconds and frequency points also in km/s. Frequencies above a line central frequency, i.e. blue shifted parts are counted as negative velocities and red shifted parts positive.

3.4 The input parameters

3.4.1 Physical cloud parameters

 Molecule to be considered [CO=1,CS=2,HCO+=3,SiO=4]:

This parameter determines which ladder shall be treated. At the moment, only CO, CS, HCO^{+}, and SiO are possible. Up to 12 levels can be taken into account for CO, CS, and SiO and up to 8 levels for HCO^{+}. Physically, the code makes no difference between the different isotopes of a given molecule since the transitional frequencies and collision rates are very similar. The only difference is the relative abundance which has to be chosen later.

 Number of shells for the cloud representation:

The cloud may be subdivided into a number of shells which have different laws for the physical parameters. Discontinuities at shell edges are possible in case of file input, however, they should be reserved to the experienced user.

 Radius of the inner core [pc]:

To avoid the singularity of power laws at zero, the user has to define an inner core. Within this core all parameters are assumed to be constant. Furthermore, the core may consist of an H II region.

 Electron density [e/cm^3] and temperature:

In case of an H II region for the inner core, these parameters have to be specified to compute the continuum absorptivity and emissivity. An electron density 0 means that no central HII region will be assumed.

For each shell the following parameters have to be specified. For all parameters except the radius, four numbers have to be given: the value at the inner boundary of the shell, the radial exponent, and possible lower and upper cut-offs. Values of 0.0 for the cuts are interpreted as no cuts.

 Outer radius of the shell [pc]:
 Hydrogen density [H2/cm^3]:
 Temperature [K]:

The kinetic temperature determines the collision rates and provides a contribution to the local line width.

 Relative molecular abundance [X/H2]:

Here, possible molecule depletions should be taken into account.

 FWHM of turbulent velocity [km/s]:

The thermal velocities will be added to get the total kinetic line width.

Turbulence correlation length [pc]:

The scale on which the autocorrelation function of turbulence goes to 1/e.

Systematic radial velocity [km/s]

Infall is counted with negative numbers.

3.4.2 Observational parameters

Upper level of the transition to be observed:

Only one line is considered at a time.

Frequency resolution [km/s]:

This might be restricted by the detector or the purpose of the observation.

Beam width [FWHP in ``]:

The beam is considered to be Gaussian.

Step size for mapping [``]:

One typically choses steps in the order of half a beam width.

Central offset for the first point [``]:

In case of observational offsets or composed maps, one may start the radial map not at the centre of the cloud but with a finite offset.

Number of map points [0=map whole cloud]:

A number 1 means just one observation, not a map. When giving a zero, the program automatically computes the number of map points up to the outer radius of the cloud.

Distance of the cloud [pc]:

3.4.3 Numerical parameters

The same set of numerical parameters can be used to determine the accuracy in the computation of the level populations and in the computation of the emergent line profiles. In a few cases a parameter is interpreted slightly different in the two parts.

Cut-off for Gaussians in the line emission [sigma]:

This one and the next parameter determine at which argument a Gaussian is considered to be zero. For higher arguments the analytic treatment is truncated. The cut-off in line emission determines at which frequencies in radiative transfer the absorptivity/ emissivity is treated as zero, i.e. the background field remains unchanged. The default value implemented in the program is 3.2.

Cut-off for Gaussians in freq. integration [sigma]:

This parameter says where to set the limits in integration over the intensity multiplied with a Gaussian profile. It is applied in the frequency integration to determine the local energy density at each point and in the spatial integration over the telescope beam. The default value is 4.2.

Resolution in scanning the Gaussians [sigma]:

The parameter defines the resolution of the grid for scanning the frequency scale in the integration of the radiative energy density. This frequency grid is locally adapted to be always at least as fine as this resolution factor times the local σ. The default value is 0.55.

Maximum change factor for the level densities between grid points:

The gas densities and level populations are linearily interpolated between the radial grid points. This parameter determines the maximum change of the level densities n_{j} between the grid points. In case of steeper changes, additional grid points are added. The factor has to be greater than 1. The default value is 1.3.

In many cases, the reduction of this parameter provides the strongest improvement of the computational accuracy at the costs of a strong increase in the time consumption. It may, therefore, be considered as the most important numerical parameter.

Hysteresis in the radial grid adjustment [0..1]:

Since the radial grid is dynamically adjusted to the gradient of the level densities n_{j}, grid points may be added or removed during the lambda iteration. The hysteresis gives the relative difference between the thresholds for addition and removal of points. It supresses grid oscillations during the iteration. The default value is 0.2.

Maximum change of the radial velocity between grid points [sigma]:

Beside the density structure, also the velocity structure may require additional radial grid points. The systematic velocity is linearily interpolated between the grid points. The parameter determines the maximum change of the systematic velocity between two radial grid points. In case of stronger changes, additional grid points are included. The default value is 2.5.

Maximum relative distance of points on the spatial integration scale:

This parameter determines the number of displacement values treated in the radiative transfer code. Since the computation of the energy density at a given radial point requires the integration over the z axis (Eq. 10), many points should be given along the z axis. The parameter constrains the maximum relative distance between two points on the z axis within this integration procedure. The default value is 0.25.

Maximum velocity shift between radiative transfer points [sigma]:

The radiative transfer code interpolates the emissivity and absorptivity for each frequency linearily within each step. In case of large velocity gradients, this will be inappropriate since the line profiles are widely shifted between two grid points. Here, additional steps have to be included. The parameter gives the maximum length of such a step on the scale of velocity in z direction v_{z}. It is certainly a good idea to take the same value as for the frequency resolution so that the frequency shift between consequent radiative transfer points is always according to at most one grid spacing at the frequency scale. The default value is 0.55.

Initial populations [1=thermal,2=no radiation,3=Sobolev]:

Here, one can determine the first guess for the level population the lambda iteration starts with. Three choices for the first guess are possible. The thermal distribution is an appropriate approach for clouds which are optically thick within at least one line of the molecule considered. The radiation free distribution should be chosen for optically thin clouds. The Sobolev approximation is good for clouds with large velocity gradients. However, I found that, due to the high efficiency of the accelerated lambda iteration method, the influence of the initial guess on the total number of iterations is relatively small. The default value is 2.

Neglection threshold for the level populations:

The parameter sets the truncation of the equation system for the level populations. Whenever the excitation of one and all higher levels troughout the whole cloud falls below this threshold, the level is neglected in the further computation. According to the machine precision used for the representation of the numbers values below 10^{−14} make no sense. The default value is 10^{−8}.

Relative accuracy as convergence criterion:

This is the convergence criterion. It is tested only after a full lambda acceleration cycle. When even the lambda acceleration does not find a larger step size than this accuracy value, convergence is reached. To obtain accuracies better than 2 10^{−7}, the program has to be compiled with real*8 numbers instead of the normally used real*4. The default value is 10^{−6}.

Due to the strong nonlinearity of the whole line radiative transfer problem it is not possible to give a general error estimate for a given set of numerical parameters. Tests using a broad range of cloud parameters have shown, however, that the error in line intensities using the default numeric parameters falls always below 5–6%.

When the command line option -numdefaults is given at program invocation, the default numerical parameters are used in all computations. This can avoid the annoying confirmation of the usage of these default values before every computation in cases when one does not plan to change them. Then, it is, however, impossible to use changed numerical parameters within the same program run when it turns out that a given problem is not soluble with these default numerical parameters.

3.5 Error messages and warnings

Due to the static declaration of fields in FORTRAN 77 the possible numbers of grid points on each scale is restricted to the maximum size of the corresponding fields. In the current version, there are at most 140 radial points, 280 different impact parameters, 223 frequency points and 12 atomic levels possible. Most errors occur when the combination of cloud properties and numerical parameters requires a grid in one of these scales which exceeds these numbers of grid points. You will find the following messages:

Number of radial points insufficient to treat the problem with the required
accuracy !

The gradient in the level densities requires more than 140 radial points either due to a very steep density gradient or due to strong changes in the level populations, e.g. at the edge of very thick clouds. The user will be brought back to the input routines to change some parameters. One might either somewhat increase the numerical parameter characterizing the maximum change of level densities between grid points or in case of extremly differing densities within the cloud substitute regions with steep density gradients by a discontinuity in the gas density.

Number of ray points insufficient to treat the problem with the required accuracy !

This error may scarely happen when the user has chosen a very small distance of the grid points on the z integration scale in connection with many radial grid points. It may also appear in the computation of the emergent line profiles when the systematic velocities are very large compared with the local random velocities or a large cloud is mapped with a very fine beam. He will be brought back to the input routines. One should slightly increase the numerical parameter characterizing the relative distance of grid points on the spatial integration scale or, in case of the large velocities, increase the numerical parameter giving the resolution in scanning the Gaussians.

Number of frequency points insufficient to treat the problem with the required
accuracy !

Looking at the maximum velocities in the cloud and the frequency resolution required either by the local line width in some points of the cloud (first part of the program) or by the observation (second part), more than 223 frequency points would be required. The user will be brought back to the input routines to change some parameters. When the error appears in the first part, i.e. the computation of the level densities, one should increase the numerical parameter characterizing the resolution in scanning the profiles or slightly reduce the numerical parameter for the cut-off in the frequency integration. However, there is no solution when the systematic velocities exceed the local line width by many order of magnitude. If the error appears in the second part, i.e. the computation of the emergent line profiles, one should consider to reduce either the frequency resolution chosen as observational parameter or to reduce the numerical parameter for the cut-off in the frequency integration.

There is no explicite test whether the maximum number of atomic levels treated by the program is sufficient for a given problem. I case of doubt one should inspect the output file for the level densities looking at the excitations of the last level. The number of levels can be a serious restriction when treating relatively warm and optically very thick clouds.

A small number of consistency checks is done for the parameter input and the following messages may appear:

The chosen transition is practically not excited !

This error occurs when the observational parameters select a transition the excitation of which falls below the neglection threshold of the previous computation of the level densities. A new transition has to be selected.

Only finite beam widths are possible !

The FWHM of the beam has to be > 0.

Nevertheless, most parameters are not checked and strange results or abnormal program termination may occur in case of unphysical or numerically absurd parameters.

4 Hints for the application of the code in line fitting

My own tests and the experiences of some colleagues lead to a few recommentations the user of the program should be aware of:

• To get a first feeling one should always start using the default numerical parameters directly implemented in the program. Only in a late stage, the accuracy of the result should be checked by setting stronger limits to the numerical accuracy.

• When fitting a number of lines, the user should always start with a simple one-shell-model (except there is a detailed physical model in background suggesting a more complicated structure). Modifying the parameters within this one shell will already provide the possibility to produce a very wide variety of line shapes and intensity ratios. In complicated models, one can easily get lost within the parameter space. Remember, in 1-D, there is never only one unique solution to a line fitting problem !

• To get a reasonable initial guess for the cloud parameters, traditional ways for the determination of the gas temperature and the column density, like the comparison of the line strengths from different isotopes of a considered molecule, should always serve as starting point.

• Whenever the cloud model covers a range of gas densities > 10^{3} or the density structure is relatively complex, it may happen that one runs out of the number of points available for the radial grid. For smooth density structures, the user should relax the numerical parameter giving the maximum difference in the level densities between two radial points. In case of strongly variing density profiles it may be more appropriate to substitute steep gradients by artificial jumps in the density (The input of discontinuities in the parameters is only possible using file input).

• At the end of a density fit, the user should check the consistency of the temperature estimate with the computed line temperature of the thick line used for the temperature derivation. I found that even quite thick lines (τ > 10) are often not completely thermalized. At least when ^{12}CO is used, this way can provide a much more accurate temperature determination.

• As a first orientation to constrain the turbulence parameter space, one might follow Miesch & Bally (1994) which got correlation lengths of about 0.1 pc for the best resolved clouds and correlation lengths up to 0.7 pc for more distant clouds. Their values for the Kolmogorov exponent γ giving the radial dependence of the turbulent velocity dispersion range from 0.74 for L1228 down to values around 0.5 for relatively many objects. For a homogeneous isotropic infinite turbulence, there is also a theoretical relation between γ and the correlation length: γ ∝ 1/ ln(r_{c}) showing a general anticorrelation.

• The approximation of effective absorption coefficients for a clumpy medium is valid both for clumping in space and velocity space. The used description by the correlation length is according to a homogeneous medium with pure turbulence in velocity space. There is no additional parameter for the filling factor in space in case of a clumpy medium, since this parameter is not independent. Instead, one has to use the clump hydrogen density for the hydrogen density parameter,the molecular abundance has to be reduced by the filling factor and the correlation length has to be increased by one over the filling factor.

5 System requirements

The program is written very close to standard FORTRAN 77 so that most f77 compilers should be able to compile it. It has been tested as it is on DEC Alpha running OSF/1, DEC VAX running Ultrix, SUN Sparc running SunOS, IBM RISC running AIX, HP with HP-UX, and Intel machines running Linux (with f77 and g77). On DEC VAX machines running VMS, four lines in ltrio.f have to be changed due to a different record treatment of this operating system compared to the Unixes above. On none of these platforms the compiler or linker should give any warnings. I did not yet succeed to run the code under MS-DOS due to a strange memory management of the Microsoft-Fortran-Compiler/Linker.

Although all fields are restricted to a relatively small number of grid points (see Sect. 3.5), the field of energy densities at all radial points has to consist of 140 × 560 × 11 real*4 numbers which makes in total about 3.5MB alone. All other fields use another 500kB for storage. Adding 500kB up to 3MB for the code itself and the dynamically linked libraries (depending on the machine and compiler) leads to a total need in RAM/Swapspace of 4 to 7 MB.

6 Limitations and bugs

At the moment the program only knows four types of molecules: CO, CS, HCO^{+}, and SiO with different isotopic constitutions. Further molecules will be added as soon as I find appropriate sets of temperature dependent collision rates for them.

Unfortunately, there is only a limited temperature range where collision coefficients were published. At higher temperatures the rate coefficients have to be extrapolated. Consequently, the computations are quite uncertain at temperatures above 300 K for CO, CS and SiO and above 70 K for HCO^{+}. Moreover, one should also have in mind that the maximum number of levels is restricted in the code so that higher temperatures might also provide a non-vanishing excitation at levels which are neglected in the computation. The user has to check for his own, that the maximum level number is sufficient.

The number of shells with different density or velocity laws is restricted to 5 in this version.

The program is not designed to fit the situation of a strongly variing turbulence. The gradient in both the absolute value of the turbulent velocity dispersion and in the turbulence correlation length should not exceed the steepest gradient of either the systematic velocity or the density.

Due to the limited number of frequency points the systematic velocities may fall at most one order of magnitude above the random velocities.

At the moment, I know of no bugs in the code.

7 Future plans

New versions of the program are planned to contain:

• a more user-friendly frontend for the parameter input

• improvement of the treatment of scale-invariant turbulence expilitely considering the back reaction on the source function

• treatment of nonlinear molecules

• …

8 Internals - Structure of the code

Here, I can only give a short overview over the general structure of the code. It shall help as a guide through the source code in case of a deeper inspection of the numerical internals or in case of error debugging. The overview explains the function of the different routines and how they are concatenated within the program. It cannot substitute for a study of the sources which I have tried to comment relatively detailed.

Please, do not perform any essential changes to the sources without contacting me and send me all your suggestions to improve the code or to make it more user-friendly.

8.1 Main program

LINETRANS is the main program. It provides the loops for calling MICRO and LINE, the routines for the computation of the molecular level populations and the observed line profiles.

Both MICRO and LINE start with reading the parameters specified by the user. MICRO then initializes the grids and starts the loop made of the radiative transfer (STPMIC), the lambda iteration (LAMBDA, LEVMIC) and the grid adjustment (ADAPTGRID). In case of convergence the results may be stored (WRITEDENS).

LINE continues with setting up the grids for beam offset and frequency, calls the radiative transfer in one line (STPLINE), the beam integration (INTBEAM) and finally writes the results (WRITELINE).

PRESENT shows the opening and closing screen of the program.

8.2 Input - Output

The routines for the parameter input READPHYS, READNUM, READOBS have the same general structure. One part runs the file input and another one the interactive input. In case of interactive input an inverse version of the file read can be used to produce a parameter file. For the input of numerical parameters (READNUM), there is an additional part so that it is possible to use the built-in default values or to use the parameters from a previous computation if called more than once.

WRITEDENS is the routine to store the level populations in a file with the structure described in Sect. 2.2. READDENS reads the populations from such a file. Finally, WRITELINE is used to store the integrated beam temperatures in a file with variable structure. The default parameters for the structure of the files written are set in PREDEFINE which reads the command line.

8.3 Radiative transfer for all lines

These routines perform the radiative transfer computation necessary to obtain the radiation energy density at each point and within each molecular transition.

STPMIC contains the loops for all displacement parameters and all radial points. For each “ray” LEFTKPL is called to initialize the fields for the intensity (with the background radiation - UBACK) and the absorption coefficient at the edge of the cloud.

To carry out the transfer between two radial grid points, at first FIRSTSTEP is called. This subroutine determines the number of real radiation transfer steps necessary between these points to have sufficiently small shifts in the line profiles. If more than one step is required NEXTSTEP is called in a following loop to compute the size of further steps.

TRANSFER is the real radiative transfer routine. It constrains the frequency range where the line profile is non-vanishing and computes the new intensity in this range. In case of small optical depths a linear approximation of the exponential functions (Eq. 3) is done, otherwise the full function will be treated assuming that κ and ϵ show a linear dependence within the integration range. CTRANSFER computes the continuum line transfer within a central H II region. It assumes, that the opacity and source function remain constant within this core.

At each radial point STPMIC calls the integration routines UNYINT and UZINT to compute the energy density.

8.4 Radiative transfer for a single line

The subroutines STPLINE, FIRSTLSTP, NEXTLSTP, LEFTKPL, TRANSFL, and CTRANSFL are reduced versions of STPMIC, FIRSTSTEP, NEXTSTEP, LEFTKP, TRANSFER, and CTRANSFER in the sense that they compute the radiative transfer only in a single line instead of all transitions important for a molecule. They are used to compute the emergent line profiles which are to be observed. No intensity integration is called. The additional routine TAUCENTRAL is used here to give a control output of the central opacity in the line centre.

8.5 Initialization and adjustment of the grids

INITRAD builds the initial radial grid. It checks the gradients in density and velocity and takes as many points as necessary to fulfil the conditions for the maximum changes specified by the numerical parameters. Additional points are included in case of discontinuities detected in the density structure. To assign the physical parameters density, temperature, systematic velocity … to a radial point VALUES is called which computes these values from the physical input parameters.

INITN initializes the level populations at the different radial grid points according to the chosen first guess approximation. In case of a central H II region, it calls UFROMHII which computes the integrated intensity at a given point due to the continuum radiation from the H II region.

INITFREQ scans through all points and determines the local line width. From this quantity the maximum frequency/velocity step size for each point is computed.

ADAPTGRID is responsible for changing the radial grid so that the differences between the points fulfil the conditions for the maximum change of the level population between two points and the maximum velocity difference. It calls ADDPOINTS which checks for the need of additional points and includes them and REMOVEPOINTS which seeks for superfluid points which will be removed. In case of adding new grid points VALUES is called which computes the gas density and velocities at the new radii.

The routines MAKEPS and PSFORLINE build the grid of displacement parameters providing the “rays” on which the radiative transfer is calculated. MAKEPS is responsible for the grid used in STPMIC, PSFORLINE builds the grid for STPLINE. Both grids differ. MAKEPS controls that for a given radial point the relative distance of the displacment points on the z scale is below the numerical parameter chosen in NUMREAD (see Sect. 2.3.3). The grid is made as fine as necessary to fulfil this condition. PSFORLINE starts with the radial grid and makes it denser for the displacement grid to have the same accuracy in covering the beam. In contrast to MAKEPS all radial points are used.

8.6 Solution of the balance equations

MATRIXUP is the routine to compute the constants used in the coefficients of the transition matrix which is build up by AMATRIX. Together with the normalization equation the transition matrix forms the linear system of balance equations (Eq. 8). LUSOLVER solves this system using an LU decomposition method. LEVMIC stores the energy densities at all radial points for three cycles and controls the full lambda iteration. After three ordinary lambda iteration cycles LAMBDA is called which computes a second order lambda acceleration step to the energy densities which is immediately applied. In case of a sufficiently small acceleration step (convergence criterion) LEVMIC terminates the iteration. HIISMTH corrects the obtained level populations within the HII region to enable a uniform treatment for all grid points.

KAPSRC is a small routine to compute the new absorption coefficients and the source functions from the level populations using the Einstein coefficients from AIJ,BIJ, and BJI. EFFKAP is the routine which obtains the effective absorption coefficients for the turbulent medium from the correlation length and the clump absorptivities. ATAU is the corresponding reduction function. In case of a central HII region, HIIKAP corrects the opacities for the most inner point to get the continuum values.

8.7 The Sobolev approximation

SOBOLEV contains the loop for all radial points and all iterations to determine the local energy densities and level populations in the Sobolev approximation. It computes the velocity gradients which are used by ESCAPE to determine the escape probability in a given direction and its first derivative with respect to n_{j}. GAUSSOB integrates the escape probabilities over all directions using a Gaussian quadrature, the grid points of which are initialized at the beginning by GAULEG. From the integrated escape probabilities SOBOLEV computes the local energy densities which are used by SOBMATRIX to obtain the level populations.

SOBMATRIX is the routine to set up the differential form of the balance equation system used for the Newton-Raphson approach to the local level populations (Eq. 23). It calls AMATRIX to get the original form of the matrix (used on the right hand side of the system) and DMATRIX for the differential form (for the left hand side). The system is resolved by LUSOLVER. SOBKAPSRC is the routine to get the absorptivities and the source function and their derivatives from the level populations. SOBEFFKAP computes the effective absorption coefficient and its derivatives for a clumpy turbulent medium.

8.8 Intensity integration

UNYINT sets the limits for the frequency integration at each radial point and multiplies the intensities with the local line profile (Eq. 10). Then, it calls QUINTGAUSS, a simple equidistant integration routine, which carries out the frequency integration. The line profile is computed by PHI.

UZINT prepares the fields for the z integration (Eq. 10). In ranges with very small steps it calls the linear integration routine DINTLIN, for larger steps it calls the cubic spline integration routine DINTCUB. DINTCUB itself makes use of the cubic spline interpolation routine SPLINE.

INTBEAM is responsible for the intensity integration over the telescope beam. It computes the range of the cloud covered by a Gaussian beam and calls SPECERR for the ϕ integration and finally DINTCUB, the cubic integration program, for the p integration. SPECERR computes an numerical approximation to the ϕ integral within Eq. (24) which is independent from the intensity.

TRANSLAT is a small routine translating the angular beam width of the telescope into the projected beam width on the cloud.

8.9 Molecular constants

The following routines contain the full “physics” of the problem. All physical quantities are specified within them. All previous routines only contain “numerics”.

INITEINSTEIN looks which molecule shall be considered and initializes the fields containing the molecule specific quantities, especially a_{j,j−1}, b_{j,j−1}, and ν_{0}. The functions AIJ, BIJ, BJI and CFREQ simply read these tables and provide the results as function values.

PROBTHERM, VTHERM, UBACK also base on the tables produced by INITEINSTEIN. These functions provide the population of a distinct level in case of thermal distribution, compute the mean thermal velocity of the molecules at a given temperature, and determine the intensity of the cosmic background radiation at the frequencies of the molecule transitions. HIIINIT computes the continuum opacity and source function for all transitions of the considered molecule from the table of frequencies and the electron density and temperature within a given HII region.

CIJ calls CIJCOH2, CIJCSH2 or CIJHCO to get the collision rate coefficients at a given temperature for the transition of interest. Before CIJ is called, the ratio of ortho- to para-H_{2} has to be determined by ORTHOPARA so that CIJ can combine the rate coefficients for collisions with both species. CIJCOH2, CIJCSH2, and CIJHCO contain the tabulated rate coefficients, use ILOCAT to find the appropriate table entries, and interpolate the values for the given temperature.

9 Source files

The following files contain the sources for the different routines:

ltr.f: LINETRANS, MICRO, LINE, PHI, TRANSLAT, PRESENT

ltrio.f: READPHYS, READNUM, READOBS, WRITEDENS, READDENS, WRITELINE, PREDEFINE

initial.f: INITRAD, INITN, INITFREQ, ADAPTGRID, ADDPOINTS, REMOVEPOINTS, VALUES, UFROMHII

stpdiff.f: STPMIC, FIRSTSTEP, NEXTSTEP, LEFTKP, TRANSFER, MAKEPS

stpline.f: STPLINE, FIRSTLSTP, NEXTLSTP, LEFTKPL, TRANSFL, PSFORLINE, TAUCEN-TRAL

matrix.f: LEVMIC, LAMBDA, KAPSRC, EFFKAP, ATAU, HIIKAP, HIISMTH, MATRIXUP, AMATRIX, LUSOLVER

sobolev.f: SOBOLEV, ESCAPE, GAUSSOB, GAULEG, SOBMATRIX, DMATRIX, SOBKAPSRC, SOBEFFKAP

einstein.f: INITEINSTEIN, AIJ, BIJ, BJI, CFREQ, PROBTHERM, VTHERM, UBACK, HIIINIT

collrate.f: CIJ, ORTHOPARA, CIJCOH2, CIJCSH2, CIJHCO, ILOCAT

ltrinteg.f: INTBEAM, SPECERR, UNYINT, QUINTGAUSS, UZINT, DINTCUB, DINTLIN, SPLINE

10 Auxiliary programs

I have written some small auxiliary programs, mainly for testing purposes around the LTR code. The program TESTCLOUD reads the level populations and energy densities from a data file produced in the first part of LTR and computes the excitation temperatures, radiation temperatures and few additional quantities. The input file must have the structure as produced by WRITEDENS, the output file is self explanatory.

The program CLINE reduces the two dimensional data set as produced by WRITELINE and extracts the profile for zero beam offset, i.e. the resulting file contains only one line profile. This may be suitable for users which do not have a sophisticated graphics front end as GREG or IDL.

11 Acknowledgements

A direct precursor of this program was written by Endrik Krügel (MPIfR Bonn). Basing on his code of about 1000 lines, I could set up this more comprehensive and flexible program. The subroutines CIJCOH2 and ILOCAT are still very similar to the original code. The temperature dependent collision rates for CO used there were tabulated by D.R.Flower & J.M.Launay (1985). The collision rates for CS and SiO, computed by Sheldon Green, were tabulated by B.E. Turner, K. Chan, S. Green & D.A. Lubowich (1992). The collision rates for HCO^{+} are based on computations of T.S Monteiro (1985) and S. Green (1975).

The local statistical treatment of clumping or turbulence roughly follows the formalism derived by H.M. Martin, D.B. Sanders & R.E. Hills (1985). The accelerated lambda method was published by Lawrence Auer in 1987. The LU decomposition method, the Gaussian quadrature, and the spline interpolation algorithm were taken from “Numerical Recipies” (Press, Flannery, Teukolsky & Vetterling 1986).

12 References

Auer L.H. 1987, In: Kalkofen W. (Ed.), Numerical Radiative Transfer, Cambridge, p. 101

Flower D.R., Launay J.M. 1985, MNRAS 214, 271

Fuller G.A., Myers P.C. 1992, ApJ 384, 523

Green S. 1975, ApJ 201, 366

Green S., Chapman S. 1978, ApJSS 37, 169

Krügel E., Chini R. 1994, A&A 287, 947

Le Bourlot J. 1991, A&A 242, 235

Martin H.M., Sanders D.B., Hills R.E. 1984, MNRAS 208, 35

Miesch M.S., Bally J. 1994, ApJ 429, 645

Monteiro T.S. 1985, MNRAS 214, 419

Press W.H., Flannery B.P., Teukolsky S.A., Vetterling W.T. 1986 Numerical Recipies (Fortran Version), Cambridge

Turner B.E., Chan K., Green S., Lubowich D.A., ApJ 399, 114