The modelling of line emission from axially symmetric structures

We have written and tested a computer modelling program to solve line-radiative transfer in any cloud with axial symmetry using a lambda interation method. It can accept any axially symmetric velocity field, thus enabling the modelling of such objects as rotating discs, oblate or prolate spheroids, filaments and outflows. We present here a full description of the method and, by way of illustration of its use, the benefits that proposed new large millimetre telescopes would have in observing rotating circumstellar molecular discs embedded in collapsing envelopes.


I N T R O D U C T I O N
A primary aim of observing star formation regions with radio telescopes is to deduce the internal composition and distribution of material within the objects of interest. Since a telescope can produce only a 2D map of the integrated emission along its line of sight through the cloud, it is necessary to deduce a 3D model of the material that is causing the emission. Although there are several analytical methods that can give some information on the structure (e.g., rotation diagrams), most of them rely on the cloud being optically thin. The most comprehensive approach is to perform radiative transfer modelling. This generally involves assuming a physical cloud structure (i.e., such parameters as H 2 distribution, temperature distribution, etc. must be assumed) and using a computer program to determine the radiation emitted by such a structure and how this would be viewed by a telescope on (or off) Earth. By varying the initial assumptions, it is then possible to match the resultant output with that observed by a real telescope and thus determine the best-fitting cloud model.
Although radiative transfer modelling using computers for astrophysical purposes was undertaken as early as the 1950s (e.g. Koelbloed 1956; see Rybicki & Hummer 1967 for a summary of early work), it was not until well into the 1970s that the main methods that are in use today for modelling spectra line formation in star formation were developed. Monte Carlo modelling (developed in its modern guise by Ulam, von Neumann and Fermi in the mid-1940s) was introduced to molecular cloud radiation transfer by Bernes (1979). At around the same time a lambda iteration method was described by Stenholm (1977) using a simplification (core saturation) developed by Rybicki (1971). The lambda iteration method employs an initial guess at the excitation temperatures to calculate the radiation field. This radiation field is then used to calculate a new set of excitation temperatures, from which in turn another set of excitation temperatures is derived. This iterative process is continued until a converged solution is obtained.
Due to the limitations of computers at the time of these developments, the initial versions of these programs generally permitted modelling only of spherically symmetric clouds (although some early Monte Carlo models used cylindrical symmetry; e.g. Stockhausen 1968or Avery 1969. Whilst for some objects this restriction is not particularly severe, for others it renders realistic modelling impossible. For this reason we have developed a much more flexible lambda iteration radiative transfer program. A generalization of a spherically symmetric lambda iteration program written by Matthews (1986), it removes most of these restrictions by allowing any axially symmetric model structure. Within this structure any desired velocity field can be specified (as long as it remains axially symmetric). This enables the modelling of most structures of astronomical relevance, e.g., oblate and prolate spheroids, discs, filaments and outflows. It differs significantly in method from the multidimensional programs that have been written by Park & Hong (1995), Juvela (1997), Srivastava (1999) and Hogerheijde & van der Tak (in preparation). These all use the Monte Carlo approach to solve for the radiation field in the cloud.
The method presented here has no restrictions on the type of cloud that can be modelled, other than the requirement for axial symmetry. We are thus able to apply it to types of astronomical objects for which it has not previously been possible to conduct accurate radiative transfer modelling.

The radiative transfer equation
We are required to solve the radiative transfer equation where f n is the function that describes the line profile, n is the frequency, and V is a solid angle. The line profile function is defined to be normalized so that f n dn 1Y and for a Gaussian profile is given by f n xe 2 n 2 n 2 0 with x 1 p p n0 X The emissivity function, E n , is defined by where A ul is the Einstein A coefficient (and the Bs are the Einstein B coefficients), and N u is the number of molecules in the upper level. The Cs are the collision coefficients between the molecule of interest and H 2 (and/or He). In general, E n will vary with distance, l, along the line of sight. Therefore, in general, there will be no analytical solution to equation (1); however, it is possible to have a computer solve the equation by splitting the integral up into sections. The number of sections is chosen to be sufficiently large such that within each the value of E n ,i is effectively constant. If this condition holds, then it is possible to integrate the individual segments to give where l i is the length of segment i which should also satisfy t i ! 1X This is now in a format that a computer program can easily handle. Using this equation, it is possible to calculate the total radiation intensity, J Å , on any point in the cloud due to all the other material in the cloud. This is done by integrating over lines of sight evenly scattered over 4p steradians as in equation (3). This can then be integrated in segments to give where DV s is the velocity step, DV d is the solid angle step size, n is the number of lines of sight (counted by d), and m is the number of velocity steps across the line profile (counted by s). The solid angle step size is dependent on the number of lines of sight, and is simply 4p/n. Due to the normalization of f n , it is possible to rewrite this as J DV s n p p n 0 n This is the total flux on a particular point in the cloud due to all the rest of the material in the cloud.

Axially symmetric geometry
The major steps in the implementation of a program to solve the radiative transfer equations in Section 2.1 for a modelled object are shown in Fig. 1. We have chosen a cylindrical coordinate system for our model [using the usual (r,f,z) convention], as this enables an easy description of most of the types of cloud that we wish to model. Our geometry units are then a series of stacked concentric rings. We consider a series of planes and cylinders given by: plane : Thus, by choosing a series of values for r c and z p , we are able to fully define the rings that make up our model cloud. The values of r c and z p are chosen by the user before the program is run, and may have any desired distribution. We have generally chosen a distribution that is weighted towards the centre of the cloud, as this is where the highest opacities occur in our models. In order to calculate the radiation field on each ring due to all other rings in the cloud, we take a series of lines of sight, calculate the radiation caused by molecules along each line of sight, and then integrate over all lines of sight to get the total radiation field. Unlike the spherically symmetric case, there are no symmetries that can be used here, and it is necessary to calculate the radiation field over the full 4p steradians.
We can describe the locus of any point on any line of sight as a combination of two components, as shown in Fig. 2: The x 1 component defines a point at the centre of the cross-section of the ring, in the x±z plane. From there the lines of sight emanate as x 2 . Each line of sight is defined simply by its value of f los and h los . This has the advantage that the x 2 are the same for each ring in the cloud and therefore need to be calculated only once. It is then necessary to find all the intersections between each line of sight and all the cylinders and planes between the start point of the line of sight and where it exits the cloud. This can be done simply by working along the line of sight and setting X c X L and X p X L X Manipulating equations (8) and (9) shows that the plane intersections are at b p z p 2 z mid ah los Y and the cylinder intersections are at b c 2r mid cos f los^ r 2 c 2 r 2 mid sin 2 f los p Y where the b p and b c are the positions along the line of sight (relative to its start position at the centre of the ring defined by r mid and z mid ) of the intersection with the appropriate plane or cylinder (as defined by z p or r c ). This can be easily converted into an actual distance along the line of sight by determining the magnitude of x 2 , i.e., jx 2 j 1 1 h 2 los p X Hence the lengths of the line segments in each ring are known.
It is also necessary to know the velocity component along the line of sight at each of the cylinder/plane intersections in order to be able to calculate the Doppler shift in the radiation emitted. The velocity field in the cloud is defined in terms of three orthogonal components.
(2) Rotational velocity ± this defines the rotation of the cloud (positive is clockwise rotation when the cloud is viewed from above' ± positive z direction).
(3) Vertical velocity ± this defines velocity along the z-axis (and can be used to define the motion in bipolar outflows, for example).
These three components can be considered as vectors V r , V f and V z . It is necessary to calculate the components of each velocity vector along the line of sight at each intersection with a disc or a cylinder (as well as at the start of the line of sight). In order to do this, the angle between the line of sight and each component of velocity must be known. This is best calculated using the dot product for vectors on each velocity vector component and the line of sight x 2 as given by For a general position along the line of sight the vector components may be written as where r mid is the starting position for the line of sight, and f los is the angle of the line of sight. These are the cylindrical coordinate system components, but their components are given in Cartesians (the dot product needs orthogonal axes). Now, by picking the required value of b, the angle between the line of sight and the velocity components at any point on the line of sight can be calculated. Using equations (10) and (11) The velocity component along the line of sight is then given by V jV r j cos u r 1 jV f j cos u f 1 jV z j cos u z X 13 The difference between the value of V at the start of the line of sight (i.e., when b 0 and its value at the intersection of interest is the relative velocity between the start position of the line of sight and the intersection of interest. This now means that each segment on the line of sight has its velocity relative to the starting point of the line of sight defined at both ends. These values can be averaged or interpolated to yield a relative velocity for the entire segment (however, as will be shown in the next section, this interpolation is not always as simple as it at first seems).

Solving the radiative transfer equation
For the first iteration it is necessary to calculate, for each ring and each transition, the new radiation field, given the assumed level populations. This is done by calculating, for each transition in each ring, the radiation intensity contributed by each line of sight.
In each case this simply involves adding up the contributions from each ring segment that the line of sight intersects. Integrating the intensity around 4p steradians gives the total radiation intensity at any one point in that ring. However, as the rings may be moving relative to one another, it is also necessary to take into account the Doppler shifting of the emitted frequencies. This is done by considering n different velocities centred around the line centre (we generally use n 11X The calculations are performed on each segment (m) on each line of sight (k) for each velocity step (l) for each transition ( j) in each ring (i). The first step is to adjust the optical depth to take into account the velocity shift due to the Doppler effect. The optical depth per unit length for each ring was calculated in the first pass through the convergence check section, and is adjusted to give: where t i,k,m is the optical depth, and v m is the velocity of the segment relative to the start of the line of sight (it is given by the average of the two relative velocities at the start and end of the segment). DV i is the velocity width of the ring in which the line of sight starts, and DV i,k,m is the velocity width of the ring in which segment lies. n is the number of velocity steps that are to be used (and must be odd), and the [ ] brackets mean`the integer part of'. So k m represents the optical depth per unit length of segment m as seen by an observer in ring i at a particular velocity shift, indexed by l. By multiplying this by the length of the segment the total optical depth for the segment at the required velocity shift is produced. Now, applying the radiative transfer equation yields the radiation intensity due to this segment on the starting position of the line of sight: where L m is the length of the segment. So I i,j,k,l,m represents the radiation intensity on ring i due to the molecules in the ring containing segment m that are radiating at velocity offset l after taking into account the absorption of this radiation by the material in segment m and all other intervening segments.
There is unfortunately a problem with this when the optical depth for a particular segment gets too large. Consider equation (15) for m 1 (i.e., the first segment on the line of sight); then if k 1 is large, the term e 2knÂLn will become very small, and hence the intensity for that segment will be lower than it should be. In an extreme case the intensity calculated would be virtually zero. The solution for this is to split the segment up into sufficient shorter subsections such that the optical depth for each segment is small enough to reduce this effect to an insignificant level. The program therefore checks at this stage to see what the optical depth for the segment is. If it is greater than 0.2, the segment is split into n knÂLn 0X2 1 1 segments. It is also necessary to split a segment up if the velocity change across the segment becomes too large.
Unacceptably large errors will occur if jv 2 2 v 1 j . 0X3 Â DVY where v 1 and v 2 are the velocities at the beginning and end of the segment, and DV is the ring velocity width. If this condition occurs, the segment is split into jv22v1j 0X3ÂDV 1 1 subsegments. 1 Some thought needs to be given to situations where this interpolation does not work well. Fig. 3 demonstrates the problem by considering a sample line of sight (the thick line) passing offcentre through the cloud (assume that it travels parallel to the planes, and thus only the cylinder intersections are of concern). At each intersection the velocity is defined. When the line-of-sight segments are split up into subsegments as described above, the velocities at the subsegment boundaries are defined by simple linear interpolation between these values. The graph on the right shows the value for the velocity at various positions in the cloud for a r 2 1 2 velocity distribution (e.g., a Keplerian disc) as a solid line. The velocity changes only slowly in the outer part of the cloud, but then rises dramatically towards the centre. The velocities for the intersections (i.e., 1, 2, 3, 5, 6, 7) are indicated by the dotted lines. The dashed line in the velocity graph joins the velocity values at these six points, and represents the velocity values that the program uses for its calculations. It is immediately obvious that this method fails badly in the central part of the cloud (i.e., between intersections 3 and 5). This is due to the intersections at 3 and 5 having identical velocities. The effect this has on the final output can be seen in Fig. 4. This is the result for a beam with 5 Â 5 25 gridding points (see Section 3.2) slightly offset from the centre of a cloud with an r 2 1 2 velocity distribution and a constant H 2 and abundance distribution. The offset is sufficient so that none of the beam reaches to the centre of the cloud (this is why only one`wing' is seen). It can be clearly seen that there is an artificial looking blip at the 26 km s 21 position. This is caused by one of the offcentre lines in the grid experiencing the effect described above; all the material in the innermost section emits at the same velocity rather than being spread out over a range of velocities, and thus there is too much emission at that velocity. The program deals with this by searching for the position along the line of sight that is closest to the central line of symmetry in the cloud (i.e., the r 0 line). At that point the innermost segment is split into two, and the correct velocity for that central point is calculated (i.e., effectively point 4 in Fig. 3 is inserted). The interpolation then functions as shown by the dashed/dotted line in Fig. 3, which is a much closer representation of the actual velocity field. The result of this modification can be seen in Fig. 4. Further improvement will come only by adding more cylinders to the model. This system is not perfect, but it deals quite well with an r 2 1 2 velocity distribution. However, there are certainly other velocity distributions where it will not work as well.
Using the above method for calculating the emission from each segment, the program then proceeds to add up all the contributions from all the segments along the line of sight. This then gives the total radiation intensity from that line of sight, i.e., The second term here adds the radiation contribution, B(n j , 2.73), from the cosmic background radiation (at 2.73 K), after taking into account its passage through the cloud. It is now possible to calculate the radiation intensity on the ring at this particular frequency by summing all the lines of sight and integrating over 4p steradians:  and hence, from the I n f n in equation (3), The radiation now needs to be integrated over the velocity range in order to finally yield the total radiation intensity on a ring due to all other rings (and the cosmic background) for each transition: Once the radiation field has been established, the new population levels based on this radiation field can be calculated through a simple matrix inversion.
In the preceding discussion it is clear that in order to deduce the radiation field correctly, sufficient lines of sight must be used to sample the cloud fully. As this will vary for different cloud structures, we have left the number of lines of sight as a variable, and as many or as few as desired can be used. The lines of sight are distibuted approximately evenly over 4p steradians. Details of the method used can be found in Phillips (1999). As we show in the appendix, provided that sufficient rings are used to cover the rapidly varying parts of the cloud, not many lines of sight are needed. For all the results presented in this paper we have used 22 lines of sight. Similarly, sufficient velocity sampling must be used. We have retained the original 11 velocity steps used in the spherically symmetric program of Matthews (1986). However, it should be emphasized that this is per position in the cloud; due to the comoving frame method employed, the number of steps across the entire velocity field of the cloud will be much higher (the actual number will depend on the cloud velocity structure).

Convergence check
The convergence check starts by calculating the new emissivity function (using equation 4) and storing the old value of the emissivity function. The convergence check is then done by comparing the old and new values of the emissivity function using for all shells k and all transitions i 3 jY where DE must be less than the convergence criteria specified in the input files for that transition to have converged. Convergence is defined as having occurred when all transitions in all shells converge twice in a row. Typical values for DE that we use are between 0.01 and 0.1, depending on the accuracy of the solution desired. Note that this only defines convergence for any given model cloud ± for the entire problem to be considered converged the solution must also be invariant to the number of rings and the number of lines of sight per ring used. Thus the values we have chosen to use for the model clouds studied here may not be appropriate for other types of structure.
Once convergence has taken place, the final step is to calculate the line profiles.

C A L C U L AT I N G T H E L I N E P R O F I L E S
There are several parts to this problem. Since the cloud is a 3D object as viewed from outside, it is necessary to describe the exact position from which it is being viewed (i.e., the output line profile will depend on the viewing position). With the position specified, the internal geometry for the lines of sight that pass through the cloud is calculated. Then, within the correct velocity window the line profile is calculated. These steps will be described in turn.

Define the line of sight from the telescope through the cloud
The telescope is located at some angle relative to the centre line of symmetry in the cloud. When observing with the telescope, observations of the cloud will be made not just at the central position but also at various offset positions, often forming a grid of observations. In Fig. 5, for example, nine positions are shown. Offset positions are also required when performing the gridded convolution necessary to simulate the telescope beam (see Section 3.2). These offset positions are referenced with respect to the telescope's coordinate system (generally RA and Dec.). This is, however, not the same coordinate system as used in the cloud.
In order to calculate the line profile that a telescope would see, it is necessary to know all the positions in the cloud that the line of sight passes through. In order to do this, the vector components for the line of sight need to be translated from the coordinate sytem used by the telescope into the coordinate system used in the cloud.
The position of the observer relative to the cloud is defined in the input files with the parameters u e and u r . This situation is shown in a side view in Fig. 5. u r is the rotation of the cloud about the y-axis of the telescope's coordinate system as shown in the diagram on the left. u e is then the rotation around the x-axis (again the telescope's coordinate system) as shown in the middle diagram. It is now necessary to find the relationship between the coordinate system drawn on the page (i.e., the x 2 z system shown in Fig. 5) and the now-tilted cloud coordinate system. Stated formally, the problem here is that we have a vector x [ R 3 (i.e., the line of sight) stated with respect to the basis B 1 (i.e., our coordinate system as viewed by the telescope), but we wish to have the vector x stated with respect to the basis B 2 (i.e., the coordinate system in the cloud). We choose B 1 1Y 0Y 0Y 0Y 1Y 0Y 0Y 0Y 1 for simplicity. If we consider this as a matrix, then it is possible to use a`change matrix' (X) to convert vectors from one basis to the other using

A XBX 17
The change matrix between the two bases is given by (using the angle notation as shown in Fig. 6) So we now have a method for converting the telescope coordinates into cloud coordinates. The line of sight from the telescope will be a vector that is parallel to the y 1 axis. It can thus be described as g., the telescope is pointed off the centre of the cloud by (15 HH ,15 HH )]. Thus by picking the required value of b, any point along the line of sight can be described by L 1 o a Y bY o b X Now, using this as matrix B in equation (17) with X as derived in equation (18) which are the components for the line of sight in the cloud coordinate system as required.
It is now possible to proceed as before and find the intersections between this line of sight and the cylinders and planes. An intersection with a cylinder occurs when the three components of X c in equation (8) equal the three components of L 2 as given above (and similarly for X p ).

Gridded convolution
The aim of the output routine is to simulate what a given telescope would see if it were pointed at the cloud. Since all real telescopes have a non-point-like beam, it is necessary to simulate the beam. This is done by placing a grid of points over the beam and calculating the emission at each of those points. These values are then added together with an appropriate weighting. This is demonstrated in Fig. 7, where the circle has a radius equal to the beam's FWHM (i.e., it is two beamwidths in diameter). The weighting factor for each position is derived by assuming the beam is Gaussian (i.e., it follows e 2d 2 , where d is the distance from the centre of the beam) and is given by v e 2 u 2 u 2 0 X 20 The absolute emission temperature is then obtained from where the T i are the emission temperatures at the ith position in the grid (with a total of n positions), and the w i are the appropriate weightings for each position. The exact number of gridding positions needed varies, depending on the type of cloud being observed (e.g., a steep power-law H 2 distribution will require more positions to cover the rapid changes at the cloud centre than, say, a Gaussian distribution). Our strategy is to start with a low number of positions (say 11 Â 11 121 and to increase this until the line profiles stop changing. Generally 21 Â 21 441 positions are sufficient, although if the inner cloud is in Keplerian rotation, more may be needed.

S O M E R E S U LT S U S I N G T H E P R O G R A M
The resultant program consists of over 4500 lines of fortran code and an additional 2000 lines of Tcl/TK code for an interactive front end. We have been running this on a 400-MHz Pentium II PC using g77 as the compiler under Linux. 2 The model results that follow took between 30 minutes and several hours to run on this system, depending on the precise problem.
We have extensively tested the program against an earlier 1D version of this program (see Matthews 1986or Heaton et al. 1993, as well as against a completely independent 1D radiative transfer program described in Kru Ègel & Chini (1994). To test the new rotational aspect of the program, we reproduced the solid body 1D rotation results presented by Buckley (1997), which are based a modified version of the 1D version of this program. We also successfully reproduced the position±velocity diagrams  presented by Richer & Padman (1991). Extensive details of the testing, a summary of which is included in the appendix, are presented in Phillips (1999).

The observability of rotation features
Thus far there are only a few examples of protostars known that exhibit features that could be caused by rotating discs. The main reason for this is that the discs tend to be quite small (100s±1000s of au), and hence are at the very limit of what present-day telescopes can observe. However, there are several new telescope projects in progress at present that should deliver dramatically higher resolutions and sensitivities in the first decade of the 21st century. It will therefore be instructive to simulate what features telescopes of differing resolutions can distinguish.
To do this, we will consider a model cloud similar to that used by Buckley (1997), which contains many of the essential features of a cloud in the early stages of star formation, but we modified it to include inner rotation ± see Table 1.
A disc is now placed in the centre of the cloud by making a disc-shaped segment rotate with a Keplerian velocity distribution ± i.e., proportional to c 2 1 2 Y where c is the distance from the axis of symmetry in the cloud. There will be no additional density enhancement in this disc (i.e., the parameters in the disc are exactly those that the non-rotating cloud has, except for the addition of rotation). The Keplerian velocity is defined by v t 1 10 c 2 1 2 km s 21 0 , c , 0X01 pcY 0 , jdj , 0X001 pc where c is the distance (in pc) perpendicular to the central axis of symmetry, and d is the distance (in pc) from the central point of the cloud along the axis of symmetry. The other parameters are as in Table 1. For the simulated observations that are presented we assume the cloud is at a distance of 140 pc as though it were in the Taurus molecular cloud. The collision coefficients for the HCO 1 modelling were provided to us by T. S. Monteiro, but are based on work published in Monteiro (1985). For CO we used rates from Green & Chapman (1978). In both cases the rates are defined at a series of temperatures; where the temperatures of our rings fall between the defined temperatures, we linearly interpolated the rates to the appropriate temperature. For the HCO 1 rates we were forced to extrapolate the rates down below the lowest tabulated temperature (by continuing the linear gradient between the lowest two temperatures). We checked that this was acceptable by running comparisons for low-J transitions with the rates published in Green (1975) which cover the low tempertaures but not the high-J transitions we require. Details of the comparison are given in Phillips (1999).
In addition, the velocity is capped at a maximum of 8 km s 21 . This may seem rather artificial, but it is necessary to stop the velocity climbing towards infinity at the centre of the cloud. We have tested different values for this cut-off, and for the resolution of the results that we present this value is sufficient to make negligible difference to the final result, since only a tiny part of the central cloud has velocities higher than this. This would not be the case for an extremely large telescope with a tiny beam.

Observing using the JCMT
It is of interest to know how the line profiles of clouds similar to the one described above are affected by looking at different transitions with various different telescopes. HCO 1 is a commonly used tracer to help deduce the structure of star formation regions. The JCMT has receivers that can detect the HCO 1 J 3 3 2 and J 4 3 3 transitions at resolutions of respectively 18 and 14 arcsec. There is also the potential for it to operate with a receiver that could detect the J 9 3 8 transition at 802 GHz, at which frequency it has a beamsize of 6 arcsec. The position±velocity diagrams for this model at these three transitions are shown in Fig. 8. Also shown for comparison in Fig. 8 is the 4 3 3 transition with a pencil beam. These diagrams show two problems associated with attempting to see such features with the present generation of telescopes. With a 15-m dish such as that on the JCMT, only at very high frequencies does the beamsize become small enough to start revealing the signature of the disc rotation. However, at these high frequencies the emission from HCO 1 is restricted to only a small area in the centre of the cloud, which in turn restricts the amount of detail about the structure that can be seen simply, because most of the disc has no emission. This problem is clearly visible in the 9 3 8 diagram ± the emission here is weak compared to that at the low-frequency transitions, due to the very high densities needed for excitation of high-J transitions of HCO 1 . It is clear from the lower right diagram in Fig. 8 that the situation is much improved by using a telescope capable of high resolution at low frequencies. This is only really feasible with an array. As most of the features of interest are fairly compact, this problem is well suited to being studied with an interferometer, since the problems of large structures should be minimized. However, as we are not yet able to model correct sampling of the u±v plane, we have not presented any simulated observations from interferometers. The 1-arcsec diagram in Fig. 8 can be considered as being representative of interferometer resolution but without any missing flux. By far the most commonly observed tracer molecule is carbon monoxide and its isotopes. Due to the numerous observable isotopes of differing abundances it is usually possible to choose an isotope that probes to the required depth in any given situation. CO also has a much lower dipole moment than HCO 1 , so higher order transitions are more easily excited. Figs 9 and 10 show the position±velocity diagrams for the same cloud as used previously, except with abundances of 1 Â 10 24 and 1X1 Â 10 26 for CO and 13 CO respectively.
The line profiles for a non-rotating cloud show the standard double-peaked signature of infall with the outer, cooler, cloud layers showing in absorption against the inner hotter regions. The symmetrical bulge in the position±velocity diagram is also a signature of infall caused by the inner cloud having high-velocity material causing the line wings.
For all the diagrams, each line profile was calculated for 101 velocity points. The number of positions where the line profile was calculated varies according to the beamsize. In most diagrams a line profile was calculated every half beamwidth; however, in some cases, in order to deal with interpolation problems with our plotting routine, it was necessary to calculate profiles at smaller spacings. For each line profile calculated the number of gridding points used to simulate the beam was either 51 Â 51 2601 or 101 Â 101 10201Y depending on the problem. This very large number is needed to ensure smooth contours around the cloud centre, where the rapid changes in velocity are taking place.
The position±velocity diagrams in general show that the features that could be caused only by rotation are difficult to observe. The line wings in these diagrams are in fact fairly easy to observe; however, they can also be caused by infall. It is asymmetric line profiles that point to rotation. Although the line wings here are significantly stronger than those in a non-rotating cloud, increasing the infall velocity would also cause the wings to become larger. Thus, unless some other means exists to determine the infall velocity, the large line wings alone do not point to rotation. The asymmetry caused by rotation is easily distinguishable from that of infall, and is caused by the opposite velocities being in the plane of the sky rather than along the line of sight as for infall. Beyond detecting the asymmetry, being able to determine the rate at which the line wing decreases away from the centre of the cloud is necessary to be able to deduce the rotation curve, which in turn yields the central mass of the object. This requires significantly better observations than just being able to detect asymmetry in the line shape. Assymetric position± velocity diagrams are unfortunately also a sign of outflow and, whilst we have not investigated this here, it is an important potential confusion. The confusion would depend quite strongly on the opening angle of any outflow, since it is expected to be perpendicular to the disc. An outflow with a narrow opening angle would thus have a negligible velocity component towards the telescope. For discs that are not edge-on this would, of course, not apply.
In the CO diagrams (Fig. 9) the rectangular shape in the lower transition diagrams here is caused by the emission being in the optically thick LTE regime. This causes flat-topped line profiles across virtually all the cloud. The only exception is close to the centre, where the high rotation velocities cause the emission to be spread out sufficiently so that it is no longer optically thick at the higher velocities, leading to the line wings. This occurs only for those beams that cover the central few arcseconds of the cloud. This effect is reduced for the higher transitions, and by the 6 3 5 transition it has disappeared altogether. These higher J transitions are promising candidates for observing, as they are producing bright lines that would be observable with the inevitably highnoise receivers at these frequencies, and yet are not so optically thick that all the information from the central regions of the cloud is lost. Additionally, of course, the higher frequencies mean a much higher resolution. CO is also unique in being the only  isotope to have the 7 3 6 transition observable in the 800± 900 GHz atmospheric window.
For the 13 CO diagrams (Fig. 10) the 2 3 1 transition (not shown) is the only transition for this isotope that suffers from the excessive optical depth problem that many of the CO transitions have. The self-absorption visible in nearly all the CO transitions is significant only in the 2 3 1 and 3 3 2 transitions here. The 6 3 5 transition especially shows a significantly different line shape near the cloud centre for this isotope. The main differences compared with a non-rotating cloud are a slightly weaker maximum and more pronounced asymmetric line wings, although the difference is detectable only at fairly weak levels.

The potential usefulness of proposed large telescopes
One such proposed telescope is the Large Millimetre Telescope (LMT). This is a 50-m telescope that is to be built through a collaboration between the Mexican Institute of Astronomy and the University of Massachusetts. The surface accuracy of the dish is to be sufficient to work at frequencies up to around 300 GHz, which we assume will enable the first three transitions of CO and its isotopes to be observed. Fig. 11 presents the predicted position± velocity diagrams for two of the observable transitions of 13 CO.
It is immediately clear, by comparing the position±velocity diagrams that the LMT would produce with those from the JCMT, that the improved resolution makes a significant difference to the amount of detail that can be gleaned from the data. The main difference is in how well the rotation curve can be determined. In the JCMT diagrams, whilst it is clear that rotation is taking place due to the obvious asymmetric wings, the variation of that speed with position in the cloud is lost due to insufficient resolution. For the higher transitions the resolution is good enough to start to see the rotation curve, but the emission is at such a low level that it would be very difficult to detect. For the LMT, however, the factor of 3 improvement in resolution makes the curve quite easily detectable even in the lower transitions. It is also noticeable that the large dish significantly increases the temperatures seen, making detection of some of the weaker features much easier.
Of course, even the LMT does not represent the best resolution that is possible. There already exist interferometers capable of working near the necessary frequencies (such as the Nobeyama Array, BIMA, Plateau de Bure and the Caltech Array), which have already produced significant results in these areas of research (e.g. Saito et al. 1996;Hogerheijde et al. 1999). The Smithsonian Array will be able to obtain vastly improved resolution images in the submillimetre spectrum. In order to show how differing resolution affects the position±velocity diagrams, Fig. 12 shows the C 18 O, 6 3 5 transition at a variety of different resolutions (as observed with a fully filled aperture). Note that this cloud was assumed to be 140 pc away. This is the closest that such objects are likely to be found. The majority of such objects will, of course, be further away, so these position±velocity diagrams can also be viewed as showing the effect that distance has on the resolution that an particular telescope may have. In Fig. 12 the diagrams are labelled as being caused by a telescope of a given resolution; Table 2 shows the equivalent distance that would produce each diagram for the JCMT, LMT, a 1-km diameter telescope and a 10-km diameter telescope. These last two are, of course, only likely to be interferometers for the near future. Note that some standard distances are: Taurus molecular cloud 140 pc, Orion molecular cloud 450 pc, Galactic Centre 8.5 kpc, Large Magellanic Cloud 50 kpc, Andromeda galaxy 675 kpc (Pasachoff 1995). The proposed ALMA will have a maximum baseline of around 10 km, which is sufficient to produce diagrams similar to that labelled 1 arcsec for objects anywhere in our Galaxy. The massive improvement in the amount of detail that can be observed using such a telescope is immediately obvious.

Effects of the viewing angle of the disc
All the position±velocity diagrams presented so far have been for edge-on discs, and thus all the conclusions about the observability of this model object rely on the object being edge-on from the Earth's point of view. In reality, of course, these objects will be arranged at random inclinations relative to the Earth. We have therefore also produced position±velocity diagrams of our model at different viewing angles. The signature for rotation is the asymmetry in the position±velocity diagram (as explained previously). The problem is in distinguishing between the radial infall which produces the bulge that can be seen in the later Figure 11. Position±velocity diagrams for 13 CO with transitions for 1 3 0 (top) and 3 3 2 (bottom) for the LMT 50-m telescope. diagrams and the almost identical shape seen in a non-rotating cloud. Surveys for rotating discs like that described in this paper would obviously be limited by the inclination effect.

C O N C L U S I O N S
In order to help draw new results from the ever higher quality data being obtained by telescopes, it will increasingly become necessary to use more sophisticated methods of interpreting these data. We have described here one such method and, by way of example of its use, have applied it to predicting the usefulness of various proposed or actual telescopes for observing rotating discs in star-forming molecular clouds.
Our major conclusion is that while higher frequency instruments, used to exploit improved resolution on present telescopes, will offer substantial further insights into the structure of star forming cores, the new generation of telescopes, in combination with better modelling such as we have developed here, should revolutionize them. Although the predicted emission is substantially decreased at these higher frequencies, the CO 7 3 6 transition none the less has considerable potential for probing the inner motions of the cores (partly due to its fortuitous location in an observable frequency band). For the lower frequencies we confirm what might initially have been expected, namely that it is necessary to select carefully an appropriate isotope that will enable the inner structure to be observed, rather than just the outer parts. For future large single-dish telescopes this is important, as they will probably not have surfaces capable of observing at very high frequencies. We suggest that the 13 CO 3 3 2 is likely to be a good choice for a LMT-sized telescope. The huge increase in detailed structure visible using ALMA class interferometers (assuming that sufficient baselines are used so that no structure is lost) shows that such instruments would prove invaluable in studying early star formation regions.
We will show elsewhere how this program can be applied to interpret actual data (Phillips & Little 2000), and expect it will help to give a better understanding of the conditions in molecular clouds as they collapse to form new stars.

A P P E N D I X A : T E S T I N G T H E P R O G R A M
Clearly, this program is of no use if the user cannot be certain that it is giving the correct results. It is therefore necessary to test the results in as many cases as possible to be as certain as possible that it is working as intended. This task is rather difficult; it is not easy to find a foolproof method of testing such a complex program, since there is no way of knowing (in the general case) whether the line profiles it produces are correct or not. We have used two main methods to test the program. First, there are a few special situations for which spherically symmetric clouds have an analytical solution. For example, the optically thick LTE case produces saturated line profiles whose intensity can easily be calculated, thus confirming correct operation in this case. A more complex test for the program is the LVG case, which also has an analytical solution. Fig. A1 presents output at different cloud offsets for a spherically symmetric cloud 3 with constant temperature, n H 2 and V infall krY where r is the normalized radius, and k 2 km s 21 X This clearly shows the correct functioning of those sections of the program that deal with velocity (especially since, in order to simulate spherically symmetric infall, our 3D velocity field has to deal with velocity components in all directions).
Beyond these special cases the only way of testing the more general cases is by comparison with the output from other programs. Since the astra program is based on the spherically symmetric sten program (Matthews 1986), which has been rigorously tested over many years, it is possible to compare results with this program by simulating spherically symmetric clouds. This we did by reproducing the results presented by Buckley (1997). His work is chosen, as the models he ran are well documented and his version of the sten program has been well Figure A1. Four positions across a model cloud which satisfies the LVG conditions. Output is shown at the centre, 1 2 Y 9 10 and 99 100 of the way to the edge of the cloud (respectively from left to right).  1997) for different abundances of CS 5 3 4 solid 1 Â 10 29 Y dashed 2 Â 10 29 Y dotted 5 Â 10 9 Y dot ÿ dashed 1 Â 10 28 X Bottom: Line profiles for the same abundances of CS 5 3 4 from the astra program.
tested against other programs that are not related to the sten program. We show just one example of this (in Fig. A2) by comparing his output with ours for the CS 5 3 4 transition in a cloud that has parameters as in Table 1, except for the relative abundance which is as in the figure caption. The cloud is assumed to be 200 pc away, and the simulations are made assuming the JCMT as the telescope (i.e., beamsize 20 arcsec). The collision coefficients are from Green & Chapman (1978) for CS, and as in Section 4.1 for HCO 1 . The slight differences in the two diagrams are probably due to the use of slightly differing collision coefficients (e.g., allowing for He collisions). We also show in Fig. A3 a test to determine the correct number of lines of sight to use. This is for HCO1 3 3 2 for the same cloud as for the CS line (i.e., infall but no rotation). The run time was shortest for the case with 22 lines of sight at 19 min, rising to 5.75 h for the case with 332 lines of sight 4 (note that the six lines of sight run time is longer than that for 22 lines of sight as, although there are fewer lines of sight, many more iterations were needed to achieve convergence due to the weak coupling of the radiation field from one part of the cloud to another).
Finally, there has been some theoretical work done by Richer & Padman (1991) on the rotation of protostellar discs. A test of differential rotation is possible by reproducing their modelled results. Fig. A4 shows the position±velocity diagram of an optically thin rotating disc with emission in a Gaussian ring. This can be directly compared with fig. 3 in Richer & Padman (1991), and represents the type of disc seen around high-mass stars such as Cepheus A and NGC 6334I.
The successful reproduction of all these results gives us confidence that the program is indeed working as intended. It should be noted that these tests are only a small (but representative) subset of all the tests performed. A much more detailed description can be found in Phillips (1999). This paper has been typeset from a 7tex;/L A T E X file prepared by the author.  Figure A4. astra reproduction of an optically thin rotating disc. 4 The algorithm used to distribute the lines of sight evenly over 4p steradians allows only certain numbers of lines of sight ± hence the strange choice of numbers.