The STENHOLM Method

The method used to solve the radiative transfer equation in the programs described in this thesis was first described in detail by L.G. Stenholm [30]^3.3. What has come to be referred to as the STENHOLM program described in section 4.2 is the latest product of several decades of advance in the writing of computer programs to simulate line generation in spherical clouds. The enormous increase in the speed of computers over the last two decades has enabled these unavoidably large and complex programs to be run on cheap computers with acceptable run times.

Stenholm described a method by which a spherical cloud or other geometrical shape could have line generation simulated by calculating the radiation intensity using either equation 3.12 or 3.13 as appropriate in what is referred to as a co-moving^3.4 frame. This is so called because the frame of reference used is different for each line of sight. It will be described below for a spherically symmetric cloud consisting of concentric shells but the method applies just as well to any other geometry.

Co-moving Frame   [l]   With reference to figure 3.1 consider a spherical cloud with many concentric shells which are collapsing inward at different speeds (labelled $V_1$, $V_2$, $V_3$, $V_4$ in figure 3.1). In order to calculate the radiation field on each shell due to all other shells a series of lines of sight must be drawn through the cloud. The radiation along these lines of sight at different frequencies can then be calculated (using equations 3.12 & 3.13), however, since the different shells are moving relative to one another the radiation will be Doppler shifted and thus it will be necessary to calculate the velocities of the various segments along the line of sight relative to one another. The question then arises as to where the velocity should be defined to be zero, ie. which frame of reference should be used. The co-moving frame always places the frame of reference at the start of the line of sight, ie. all velocities are relative to the start of the line of sight which then of course has a velocity of zero. This is shown in figure 3.1 where the line of sight starts at the small filled circle, this is given a velocity of zero. A program would then calculate the component of the shell velocities that are aligned with this line of sight at each line of sight/shell intersection.

During the 1970s as the first serious attempts to write programs to model radiative transfer were made the primary concern was keeping the number of calculations down to a level where computers of the time could easily run the models. To this end Stenholm used a technique developed by Rybicki several years earlier which he had termed the core-saturation technique [28]. This relies on the simplification that when the optical depth along a line of sight becomes large (ie. $\tau \rightarrow \infty$) the $\frac{dI_\nu}{d\tau_\nu}$ term in the equation of transfer as given in equation 3.3 tends to zero. The physical explanation for this is that at very high optical depths the radiation intensity is dependent only on the local surroundings and does not vary with position in the cloud since the photons are effectively immediately reabsorbed before they have a chance to travel any distance through the cloud. This means that equation 3.3 can be simplified to

$$I_\nu=S_\nu$$ (3.12)

Core Saturation     [r] The condition of large optical depth occurs in what Rybicki referred to as the 'core' of a line as shown in figure 3.2. Within this region it is possible to calculate the radiation intensity very quickly thus saving considerable computing time. The width of the region is defined by a parameter which Stenholm calls $\lambda$. The approximation is used when the optical depth for a particular line rises above $\lambda$. Outside of this core-saturation region it is necessary to solve equation 3.3 for $I_\nu$. Note that Stenholm uses a slightly different version of 3.3, namely $\frac{dI_\nu}{d\tau_\nu}=I_\nu-S_\nu$. The difference is in the choice of the positive direction, in equation 3.3 the direction is taken as positive from the edge of the universe to the observer. In this case it is from the observer to the edge of the universe which has the effect of reversing the signs. Therefore the integrating factor is $e^{-\tau}$ and thus the integrated equation becomes

$$\int^{\tau_\nu}_{0} \frac{d}{d\tau_{\nu}}\left(I_\nu e^{-\tau... ...=-\int^{\tau_\nu}_{0} S_{\nu}(\tau)e^{-\tau_{\nu}} d\tau_{\nu}$$

Since $S$ is now a function of $\tau$, only the left side can be easily integrated to give

$$I_\nu e^{-\tau}-I_{\nu_0} =-\int^{\tau_\nu}_{0} S_{\nu}(\tau)e^{-\tau_{\nu}} d\tau_{\nu}$$

Now since the emission at the edge of the universe can be considered to be zero the value of $I_\nu$ is zero and thus this can be written as

$$I_{\nu_0}=\int^{\tau_\nu}_{0} S_{\nu}(\tau)e^{-\tau_{\nu}} d\tau_{\nu}$$ (3.13)

This integral represents the emission along a line of sight between the observer and the edge of the universe. It is clearly not possible to find an analytical solution when $S=f(\tau)$ (unless it is a very simple function, which in reality is rarely the case). The aim of a radiative transfer program is to model the movement of radiation through a cloud as governed by this equation.

The programs described in this thesis are based on this co-moving frame method. Figure 3.3 provides an overview of the method employed. These programs use lambda iteration to find a solution to the radiative transfer equation. The first step is to give the program information on the physical conditions in the cloud (ie. H$_2$ density, molecular abundance, etc.). The cloud geometry is then calculated. For simple geometries this means calculating the intersections of the lines of sight that will be used later and storing the positions of the intersections and the relative velocities between different sections on each line of sight. For more complex geometries (this generally means non-spherically symmetric models) the size of the arrays needed to store all this information becomes prohibitive and so this data is calculated only as needed. The lambda iteration section, really the heart of the program, calculates the radiation field within the cloud by considering the radiation falling onto each geometry unit (eg. shell for a spherically symmetric cloud) due to all other geometry units in the cloud. This radiation field will be determined by the level populations in each geometry unit (which was determined by the initial input parameters). Once the radiation field is known it is possible to use the statistical equilibrium equations to calculate the new population levels. Statistical equilibrium assumes that there are an equal number of molecules leaving a particular level as there are entering it - ie. the system is in equilibrium. Therefore, given a radiation field, knowledge about how collisions between molecules cause level changes and the Einstein coefficients it is possible to calculate the level populations necessary to keep the system in equilibrium. Once this has been done the program loops back to re-calculate the radiation field in the cloud using the new level populations. This is the first loop of the lambda iteration. The idea is that the differences between the level populations at the end of each iteration gets smaller and smaller eventually causing no further change in the level populations between iterations. When this situation is reached the lambda iteration section of the program has been solved. The final stage is to determine the line shapes that would be seen by a telescope viewing the cloud. This again involves solving either equation 3.12 or 3.13 for lines of sight from the observer through the cloud.

 

In order to avoid repetition a more detailed description of the exact method and equations used in the programs described in this thesis will be given in chapter 4 (specifically section 4.7) where the ASTRA program is described.