Modelling submillimetre spectra of the protostellar infall candidates NGC 1333 – IRAS 2 and Serpens SMM4

We present a radiative transfer model, which is applicable to the study of submillimetre spectral line observations of protostellar envelopes. The model uses an exact, non-LTE, spherically symmetric radiative transfer ‘Stenholm’ method, which numerically solves the radiative transfer problem by the process of ‘ L -iteration’. We also present submillimetre spectral line data of the Class 0 protostars NGC 1333–IRAS 2 and Serpens SMM4. We model the data using the Stenholm code. We examine the physical constraints which can be used to limit the number and range of parameters used in protostellar envelope models, and identify the turbulent velocity and tracer molecule abundance as the principal sources of uncertainty in the radiative transfer modelling. We explore the trends in the appearance of the predicted line proﬁles as key parameters in the models are varied, such as infall velocity proﬁle, turbulence and rotation. The formation of the characteristic asymmetric double-peaked line proﬁle in infalling envelopes is discussed. We ﬁnd that the separation of the two peaks of a typical infall proﬁle is dependent not on the evolutionary status of the collapsing protostar, but on the turbulent velocity dispersion in the envelope . We also ﬁnd that the line shapes can be signiﬁcantly altered by rotation. Fits are found for the observed line proﬁles of IRAS 2 and SMM4 using plausible infall model parameters. The density and velocity proﬁles in our best-ﬁtting models are inconsistent with a singular isothermal sphere model (SIS), since for both objects modelled, the infall velocities appear further advanced than a SIS model would predict, given the density proﬁle. We ﬁnd better agreement with a form of collapse which assumes non-static initial conditions, in agreement with other recent ﬁndings. We also ﬁnd some evidence that the infall velocities are retarded from free-fall towards the centre of the cloud, probably by rotation, and that the envelope of SMM4 is rotationally ﬂattened. an outer radius of 6684au and is assumed to lie at a distance of 450pc. The adopted hydrogen number density and velocity proﬁles are of the form n H 2 (cid:133) r (cid:134)/ r 2 3 = 2 and v (cid:133) r (cid:134)/ r 2 1 = 2 respectively, with n H 2 (cid:136) 1 : 1 (cid:2) 10 5 and v (cid:136) 2 0 : 237kms 2 1 at the outer edge of the cloud. The kinetic temperature is set to 10K up to the half-point radius, rising linearly to a value of 20K at the outer edge of the cloud. The microturbulent velocity width and the CS abundance are assumed to be constant throughout the cloud, with values of 0.5kms 2 1 (FWHM) and 2 (cid:2) 10 2 9 respectively.


I N T R O D U C T I O N
The process of star formation is currently not well understood. However, the main protostellar collapse phases of low-mass stars (,0:5-2 M ( Þ have been identified observationally and labelled as Class 0 (André, Ward-Thompson & Barsony 1993) and Class I (Lada & Wilking 1984;Lada 1987) protostars. These are believed to represent the phases during which the circumstellar envelope accretes on to the central protostar and disc (e.g. André 1994;Ward-Thompson 1996). The final pre-main-sequence stages of Classes II and III (Lada & Wilking 1984;Lada 1987) correspond to the classical T Tauri (CTT) and weak-line T Tauri (WTT) stages respectively (André & Montmerle 1994). These protostellar stages are beginning to be understood, at least in outline (for a review, see André, Ward-Thompson & Barsony 2000).
Protostellar infall has been reported in Class 0 sources by a number of authors (e.g. Zhou et al. 1993;Ward-Thompson et al. 1996). However, the manner of the collapse remains a matter for debate. The ideas of static initial conditions for collapse (e.g. Shu 1977) have been disputed by many authors (e.g. Foster & Chevalier 1993;Whitworth et al. 1996). There is now a growing body of evidence which suggests that the collapse occurs from non-static initial conditions, and at a non-constant accretion rate which decreases with time (e.g. Kenyon & Hartmann 1995;Henriksen, André & Bontemps 1997;Safier, McKee & Stahler 1997;Whitworth & Ward-Thompson 2001). Infall has even been detected in some starless cores, such as L1544 (Tafalla et al. 1998).
The densities and temperatures in the gas envelopes surrounding the youngest protostars are favourable for exciting a number of rotational molecular transitions, observable in the submillimetre waveband. The line profiles of these transitions contain information about both the physical state and dynamics of the envelope gas, and may potentially be used to test theoretical models of star formation, as many workers have shown (e.g. Bernes 1979;Rybicki & Hummer 1991;Choi et al. 1995;Juvela 1997;Park & Hong 1998;Hogerheijde & van der Tak 2000).
In this paper, observations are presented of the protostellar candidates NGC 1333-IRAS 2 and Serpens SMM4, in transitions of HCO 1 , H 13 CO 1 , CS, CO, 13 CO and C 18 O. The HCO 1 and CS transitions preferentially trace high-density gas, whereas CO traces a much wider range of gas densities. The paper is laid out as follows: Section 2 introduces the L-iteration method of numerical radiative transfer; Section 3 describes our approach to the modelling; Section 4 explores the sensitivity of the model to the various free parameters; Section 5 describes our observations and data reduction; Sections 6 and 7 present the results of our observations for NGC 1333 -IRAS 2 and Serpens SMM4 respectively; Section 8 compares the observations with the model predictions and finds the best fits to the data; Section 9 presents a brief summary of our main findings. The non-expert reader may prefer to read the second half of the paper first, starting from Section 5.

N U M E R I C A L R A D I AT I V E T R A N S F E R
Consider a cloud of gas with a specified distribution of density, kinetic temperature and composition, which may have internal turbulent and systematic motions. Let any radiation sources not forming part of the cloud itself also be specified. For each molecular species, there exists a steady state solution for the distribution of rotational energy level populations as a function of position in the cloud. In the modelling of this paper we numerically calculate this distribution for the idealized case of a spherically symmetric model cloud, and predict the observed line profiles. This is a complex problem, due to the fact that well-separated parts of the cloud can interact radiatively with each other. The L-iteration method, described below, is conceptually one of the simplest techniques for solving this kind of problem.

The L-iteration method
The method is begun by choosing an initial radiation field in a more or less arbitrary manner. From this a 'false' run of level populations may be calculated, by assuming the validity of the steady state rate equations. Radiative transitions between these level populations will generally produce a radiation field which departs from the one originally assumed. If this radiation field is used to calculate a new set of level populations, and the procedure is repeated a sufficiently large number of times, the radiation field and level populations should eventually converge on a mutually consistent solution. Literation is a kind of diffusion process, where imbalances in the radiation field are smoothed out over a length-scale corresponding to approximately one optical depth at each iteration step.
A number of radiative transfer models have been published (e.g. Rawlings et al. 1992;Zhou 1992;Walker, Narayanan & Boss 1994). However, the modelling in this paper was carried out using a modified version of the Stenholm L-iteration code developed by Stenholm and subsequently expanded by Little and co-workers at the University of Kent (Stenholm 1977;Matthews 1986;Heaton et al. 1993). The code uses the above method to solve the spectral line radiative transfer problem for the rotational transitions of linear molecules in a spherically symmetric model cloud. Radial profiles of systematic velocity, temperature, density, tracer molecule abundance and microturbulent velocity dispersion may be specified.
The model cloud is discretized using a number of spherical shells, and the level populations in each shell are determined in the iteration process from the calculated mean radiation intensity in the comoving shell frame. Once the level populations have converged, a calculation is made of a simulated spectral line observation on the cloud. This part of the code carries out straightforward numerical integrations of the equation of transfer along parallel lines of sight through the cloud to find the emergent intensities for a grid of impact parameters, which are then weighted according to a twodimensional Gaussian beam profile function.

Testing for convergence
It is important to find a reliable test for the convergence of the model. One possible condition is that the maximum fractional change h(i, J) N , of the level population n(i, J) between iterations N 2 1 and N, is less than some specified value a: hði; JÞ N ¼ nði; JÞ N 2 nði; JÞ N21 nði; JÞ N , a; for all i and J, where the subscripts specify the iteration numbers. In many cases, the value of h(i, J) is a good estimate of the absolute fractional error in the corresponding level population, in which case convergence is ensured by setting a suitably low value of a (a value of 0.02 or less is usually found to be sufficient). However, when the rate of convergence is particularly slow, this test may give a misleading result, since the fractional change between two successive iterations may then significantly underestimate the total change produced over a large number of subsequent iterations. The rate of decrease of h should somehow be taken into account. Dickel & Auer (1994) use the following procedure to estimate the total fractional error h tot (i, J) N of the level population n(i, J) after the Nth iteration: h tot ði; JÞ N .
hði; JÞ N ½hði; JÞ N / hði; JÞ N21 2 1 : The convergence criterion is then h tot ði; JÞ , a for all i and J, for some specified a much less than 1. Note that if hði; JÞ N . hði; JÞ N21 ; then h tot (i, J) N is much larger than h(i, J) N -i.e., the total error is much larger than the fractional change between successive iterations. The rate of convergence is fastest when the optical depth of the cloud is small, or when the hydrogen number density is much greater than the critical thermalization density. When the optical depth is small, radiation can travel across the whole cloud in one iteration, allowing the level populations in different parts of the cloud to adjust to each other quickly. When the density is close to, or above, the critical density for the principal thermally accessible transitions of the molecule, the level populations are determined mainly by molecular collisions, and are insensitive to changes in the radiation field.
For optically thick regions where the density is below the critical density, the convergence rate is generally much slower. The low density results in low collision rates, and the level populations therefore depend more sensitively on the radiation field, while the high optical depth means that the radiation field will take many iterations to adjust to the conditions in distant parts of the cloud. The code converges most rapidly near to LTE and when optically thin.
We have made several alterations to the code since it was made available to us. These have included changes to the internal structure of the code by, for example, improving the shell subdivision algorithm and incorporating rotation into the line profile calculation. However, the former did not affect the conclusions -it merely gave the code more resolution. Some of the more significant changes made are now discussed in detail.

The Rybicki approximation
An important modification made to the radiative transfer calculation was the removal of the approximation of core saturation used in the original code, known as the Rybicki approximation (Rybicki 1972(Rybicki , 1984Stenholm 1977), which may be summarized as follows: suppose that at some point in the cloud the following relation for the optical depth t n at frequency n is satisfied for all directions r 0 and solid angles kˆ: t n ðr 0 ; z b ; n;kÞ . g; where z b is the distance to the boundary of the model, and g is some prescribed constant. The intensity I in all directions at that point is defined in the Rybicki approximation to be equal to the local source function S, i.e., I n ðr 0 ; nÞ ¼ S n ðr 0 ; nÞ: This equation is only applied for the particular values of n for which t n . g. For increasing values of g, the volume of the cloud to which the approximation applies becomes smaller, and disappears altogether when g increases beyond the maximum optical depth of the cloud. For the problems we consider, the source function cannot be assumed to be uniform over distances corresponding to a few optical depths, and the Rybicki approximation cannot be used. We therefore do not use the Rybicki approximation in the modelling in this paper.

Solid-body rotation
We have incorporated rotation into the code, in the form of 'solidbody' rotation -i.e., the angular velocity of all points in the cloud about some fixed axis is a constant, but systematic radial motions in the rotating frame are still allowed. The radiative transfer effects of solid-body rotation were fully accounted for by changes to the section of the code which calculates the emergent line profile from the cloud, after the L-iteration solution has been found. No change needs to be made to the L-iteration code itself, since the solution for the level populations in the cloud does not change when solidbody rotation is introduced, as is now explained. The radiation field in a frame comoving with the gas at a given point in the model can be affected by a superimposed velocity field only if the latter introduces changes in the relative line-of-sight velocities of other parts of the cloud (ignoring relativistic effects of order v 2 / c 2 or less). Solid-body rotation does not introduce changes in relative line-of-sight velocities from one part of the cloud to another. For example, consider a non-rotating cloud which has systematic radial velocities v r (r 1 ) and v r (r 2 ) at the positions r 1 and r 2 respectively. Writing r 12 ¼ jr 2 2 r 1 j for the distance separating the two points, the component of relative velocity, v los , along the line of sight between r 1 and r 2 is given by v los ¼ ðr 2 2 r 1 Þ : ½v r ðr 2 Þ 2 v r ðr 1 Þ/ r 12 : Now assume the cloud is made to rotate in solid-body fashion with angular velocity V. Let the new systematic velocities at positions r 1 and r 2 be v(r 1 ) and v(r 2 ), where Substituting into equation (5), the terms in V disappear, and the final expression is unchanged. This implies that the solid-body rotation has no effect on the relative line-of-sight velocities between any two points in the cloud, and hence the radiation field in the comoving frame at each point in the cloud is not affected by solid-body rotation. This result applies only to the radiation emitted by the cloud itself. Radiation entering the model from outside will be Doppler-shifted in the gas frame by the rotation, but the effect of this is completely negligible for realistic external radiation fields. However, solid-body rotation does make a difference to the prediction for the observed line profile, since the rotation does introduce changes in the line-of-sight velocities of different parts of the cloud relative to an observer outside the cloud. If ẑ is the unit vector parallel to the line of sight from the observer to the cloud, and v z (r) is the line-of-sight velocity (relative to the observer) of the gas at the position vector r (relative to the centre of the cloud), then v z ðrÞ ¼ẑ : v r ðrÞ 1ẑ : ðV Â rÞ ¼ v r;z ðrÞ 1 yV x 2 xV y ; where v r,z (r) is the component along the line of sight of the radial component of the gas velocity at r, and x and y are, as before, spatial coordinates in the plane perpendicular to the line of sight at the distance of the cloud. V x and V y are the components of the angular velocity vector in this plane. The component of angular velocity parallel to the line of sight, V z , does not appear in equation (8) and has no effect on the observed spectrum from the cloud.

LTE comparison
Analytical solutions of the equation of transfer are in general unobtainable for all but the simplest of problems. However, in the limit that the level populations are completely thermalised by collisions -i.e., in local thermodynamic equilibrium (LTE) -the spectral line calculations from the numerical model may be compared directly with analytical predictions. In LTE, the excitation temperature for each transition approaches the kinetic temperature T kin of the gas, and the source function for the transition J ! J 2 1 becomes S n (r) ¼ B n [n J , T kin (r)]. If T kin is constant over the cloud, then the solution to the equation of transfer is: where t n (v z ) is the optical depth of the cloud at the frequency corresponding to the velocity bin v z . For simplicity, we consider the intensity emerging from a uniform density isothermal cloud layer of thickness L, with no systematic motions, and with a constant velocity dispersion (s v ) and total tracer molecule number density (n tot ). The optical depth through the cloud layer may then be written as: where n(J) and n(J21) are the number densities of tracer molecules in the upper and lower levels of the transition, and other symbols take their usual meanings. In the last equation we have used the fact that in LTE nðJ 2 1Þ nðJÞ ¼ gðJ 2 1Þ gðJÞ expðhn J /kT kin Þ; ð13Þ and also The explicit equation for f(v z ) is given by and the number density n(J) of tracer molecules in the level J may be evaluated by: The sum in the denominator is the partition function, and usually converges rapidly after about 10 terms (depending on T kin ), and can be easily evaluated to the required accuracy. Substituting these expressions into equation (9), with I ext n ðJÞ ¼ B n ðn J ; 2:73 KÞ, gives an expression for the predicted intensity as a function of v z . To compare these predictions with the Stenholm code, we used a static model cloud with a uniform kinetic temperature, hydrogen number density, velocity dispersion, and abundance. The hydrogen number density was set several orders of magnitude above the critical density to enforce LTE. Fig. 1 shows comparisons of the analytical and numerical predictions for the CSðJ ¼ 5 ! 4Þ line, from LTE models with progressively increasing CS abundance, covering a wide range of cloud optical depths. The excellent agreement between the predicted spectra in all parts of the line, and across the large range of optical depths, is a good indication that the integration of the equation of transfer along the line of sight is being carried out correctly.

Non-LTE comparison
The LTE comparison shown above, while instructive, does not test the ability of the code to carry out non-LTE radiative transfer calculations reliably. We have not found any exact analytical solutions for non-LTE radiative transfer problems in spherical geometry, to compare with the code in a similar manner as shown above. Instead, we attempt to reproduce the line profiles obtained in a previous non-LTE radiative transfer study (Krügel & Chini 1994, hereafter KC), which used a somewhat different numerical method. Fig. 2 shows a comparison of profiles taken from KC with the corresponding profiles calculated by our Stenholm code. The figure caption gives details of the model cloud used, which has strong velocity and density gradients towards the centre, and is far from LTE. Both calculations use the CS collision rates of Green & Chapman (1978). Overall, the agreement between the two sets of line profiles is very good. The intensities of the line profiles predicted by Stenholm are up to 10 per cent weaker than the KC profiles, and the disagreement appears to be greatest for the smallest beam sizes. Comparison of analytical and numerical predictions of the CSðJ ¼ 5 ! 4Þ spectra from four uniform LTE models with progressively increasing CS abundance. The solid and dashed lines show the numerical and analytical predictions respectively. Each model has a radius of 10 000 au, a uniform kinetic temperature of 20 K, a uniform hydrogen number density of 10 10 cm 23 , and a uniform FWHM velocity width of 1.0 km s 21 . The uniform CS number density in the first model is 10 25 cm 23 , and increases by factors of 10 in successive models, to 10 22 cm 23 in the last model. The microwave background intensity has been subtracted from both sets of spectra.
We checked our calculations, and found that this was not the result of inadequate sampling of the lines of sight used in the beam convolution, by recalculating the profiles for a range of different beam sampling densities. The small disagreement in the line profile strengths is not serious, and may be caused by differences in the spatial, angular or frequency discretization schemes used in the models, or small differences in the precise values of some of the parameters used. Particularly encouraging is the very good agreement in the shapes of the corresponding line profiles.

A P P R O AC H T O T H E M O D E L L I N G
In this section we discuss our approach to the modelling of infall line profiles, and investigate the qualitative dependence of the predicted line profiles on a number of different model parameters. The approach to the modelling is necessarily a compromise between the desire to represent realistically the object being modelled, and the practical need to limit the parameter space and running time of the numerical simulations. The observations are modelled using the spherically symmetric microturbulent Literation radiative transfer code described above. Even with the considerable simplification of spherical symmetry the number of model parameters is formidable.
Radial profiles of density, kinetic temperature, molecule abundance, systematic velocity, and microturbulent velocity width must be specified. Additional parameters include the choice of inner and outer boundaries of the model, the incident radiation fields at the inner and outer boundaries, and the distance to the object. Some assumptions must be made about the forms taken by these profiles, so that they may be described by specifying only a few parameters. Some of the profiles are assumed to have a powerlaw dependence on radius, which requires only the power-law exponent and a normalization factor to be specified. Wherever possible, the choice of parametric forms used to describe the radial The model cloud has an outer radius of 6684 au and is assumed to lie at a distance of 450 pc. The adopted hydrogen number density and velocity profiles are of the form n H2 ðrÞ/r 23=2 and vðrÞ/r 21=2 respectively, with n H2 ¼ 1:1 Â 10 5 and v ¼ 20:237 km s 21 at the outer edge of the cloud. The kinetic temperature is set to 10 K up to the half-point radius, rising linearly to a value of 20 K at the outer edge of the cloud. The microturbulent velocity width and the CS abundance are assumed to be constant throughout the cloud, with values of 0.5 km s 21 (FWHM) and 2 Â 10 29 respectively. profiles is guided by theoretical considerations. We now discuss the profile of each physical parameter in turn.

The systematic velocity profile
Following the onset of fully dynamical collapse, the radial velocity profile should tend towards a freefall profile (e.g. Larson 1969;Hunter 1977;Shu 1977;Zhou 1992;Foster & Chevalier 1993). The pure free-fall velocity profile may be written as The infall velocity profiles in non-rotating clouds which collapse from different initial conditions should approximately converge, even when the mass accretion rates are very different, as long as the comparison is made at times when the central protostellar masses are equal (Hunter 1977). Hence we use the singular isothermal sphere (SIS) model velocity profiles to parametrize the infall velocity field in the radiative transfer modelling, which allows the velocity field to be described in terms of only two quantities: the effective sound speed a eff and the infall radius r inf . We stress that these are merely used as parameters characterizing the velocity field, and their actual physical significance is model-dependent. We have found the following empirical fit to the exact velocity profile predicted by the standard SIS model, for r , r inf : This equation fits the SIS collapse field to within a few per cent. If the collapsing cloud has some initial rotation, then centrifugal support may become important at small radii, leading to a decrease in the radial velocity close to the centre of collapse. With an initial angular velocity V, then we find the inner radius, r d , at which centrifugal forces begin to dominate depends on the mass already accreted, M acc , and the radius R enclosing this mass (Terebey, Shu & Cassen 1984) through the following relation: 3.2 Microturbulent velocity profile As discussed above, the radiative transfer code used to model the data in this paper makes use of the microturbulent approximation. This assumes that random 'turbulent' motions can be treated in the same way as thermal molecular motions, by incorporating them into the local line profile function. Comparisons between radiative transfer codes which model macroscopic and microscopic turbulence suggest that the main effect of relaxing the microturbulent approximation is to reduce the strength of selfabsorption features (e.g. Park & Hong 1995).
Although observations appear to indicate that random motions of some kind must be present in protostellar envelopes, very little theoretical work has been done on the role of turbulence and turbulent support in protostellar collapse. Lizano & Shu (1989) derived a 'turbulent equation of state' using the well-known 'Larson's Laws' (Larson 1981), which relate the observed velocity dispersion (s), size (R), and average density (r¯) in molecular clouds and cores: These correlations apply to ensembles of clouds and also to observations of individual clouds in different molecular tracers, and are associated with the turbulent motions in the clouds (e.g. Myers 1983). Combining the two relations above gives Zhou (1992) used this relation to infer the variation of turbulent velocity with density in his radiative transfer calculations of infall line profiles. Vázquez-Semadeni, Cantó & Lizano (1998) have recently argued that Larson's Laws may be inappropriate to dynamically collapsing clouds, because they were established from observations of clouds which were in equilibrium configurations. These authors carried out numerical hydrodynamic and magnetohydrodynamic simulations to test the behaviour of the turbulence during the collapse, and found that, in contrast to equation (22), the turbulent velocity dispersion actually increases with density during the collapse.
For purely hydrodynamic and weakly magnetic collapse models the relation between the turbulent velocity dispersion and the density was seen to approach a power law of the form s / r 1=2 (i.e., P tb / r 2 Þ. In the case of strongly magnetically inhibited collapse, the power-law index was observed to change to s / r 1=4 ðP tb / r 3=2 Þ, consistent with the predicted behaviour for the slow compression of Alfvén waves (McKee & Zweibel 1995). The behaviour of the turbulence in infalling protostellar envelopes is thus theoretically rather uncertain. Turbulence nevertheless plays a very important part in the formation of line profiles, and is one of the principal sources of uncertainty in the infall line profile modelling.

Density profile
In a region where the infall velocity has a free-fall profile there is often a strong tendency for the density profile to adopt an r 23/2 distribution, since if the rate of injection of material into the accretion stream is not changing rapidly with time, the quantity 4pr 2 rv r should be approximately independent of radius. The prediction for the density profile deep inside the infalling region may be written as where, as before, r inf ¼ a eff t is the infall radius, and a eff is the effective sound speed. The exact model density profile inside the whole infalling region is very well approximated (within ,1 per cent) by the following empirical formula: Outside the infall radius, the density profile is given by There are several competing theoretical models which predict different forms of the density profiles at various stages in the contraction and collapse of a dense cloud core (e.g. Foster & Chevalier 1993;Basu & Mouschovias 1994;Whitworth & Ward-Thompson 2001). However, once a central protostar has formed, there is a tendency in the theoretical models for the forms of the density and velocity profiles to converge towards the SIS model forms, in particular towards the centre of the collapse (e.g. Foster & Chevalier 1993). We therefore initially search for fits to our observations using density profiles of the above form, although we do not require that the adopted density profile uses the same values for the parameters r inf and a eff as used in the velocity profile.

Kinetic temperature profile
The temperature of the gas in an infalling protostellar envelope is determined by the balance of heating and cooling processes. Ceccarelli, Hollenbach & Tielens (1996) modelled these processes to derive the temperature profile in protostellar envelopes with central embedded heating sources. They found that in all the cases studied, the gas temperature profile stays within ,30 per cent of the dust temperature, usually lying slightly below the dust temperature, T d . In the following it is assumed that the gas temperature is well coupled to the dust temperature, although this assumption may be significantly in error for hydrogen number densities below ,10 4 cm 23 .
If the emitted radiation reaching the optically thin part of the envelope peaks at a wavelength l * , then the temperature of the dust at a radius r will be given by The value of l * depends on the details of the dust radiative transfer through the optically thick inner region (i.e., how the material in this region is distributed). Kenyon et al. (1993) used the diffusion approximation in the optically thick region to estimate the value of l * . For a density profile of the form rðrÞ ¼ r 0 r 2q this gives Substituting this expression into the above equation and setting b ¼ 1:5 and q ¼ 1:5 (appropriate for a hydrodynamically collapsing region), we find The temperature profile in the optically thin part of the envelope is therefore only weakly dependent on the density normalization and the luminosity (see Kenyon, Calvet & Hartmann 1993 for further discussion). For the objects modelled in this paper, the optically thick region is expected to be very much smaller than the beam size of the observations, and is comparable to the inner shell radius used in the radiative transfer models. At sufficiently large radii, the gas temperature will approach the ambient temperature of the larger scale cloud. We therefore adopt the temperature profile appropriate to the optically thin limit ðT / r s , where 20:4 , s , 20:33Þ in the inner region, and switch to a flat profile at the radius where the power-law inner profile falls below the ambient cloud temperature.

Abundance profile
The profiles of CS and HCO 1 relative abundance in an infalling protostellar envelope are governed by time-dependent chemistry and the transfer of molecules between the gas phase and the ice mantles surrounding dust grains. This problem has been studied by Rawlings et al. (1992) and Willacy, Rawlings & Williams (1994) in the case of isothermal collapsing clouds. The objects studied in this paper contain internal stellar heating sources, and are unlikely to be well described by an isothermal model. The chemistry of internally heated infalling protostellar envelopes has been studied by Ceccarelli et al. (1996). We refer the reader to all of these references for further discussion of the way in which chemistry can affect observed abundance profiles.
Given these uncertainties and the requirement to limit the parameter space, we used flat abundance profiles in our modelling. The fractional abundances of CS and HCO 1 measured in dense molecular gas cores occupy similar ranges, typically between 10 29 and 10 28 (e.g. Irvine, Goldsmith & Hjalmarson 1987;McMullin, Mundy & Blake 1994a;Blake et al. 1995). The measured values of ½ 12 C/ 13 C at the Galactocentric radius of the Sun (which applies to all of our observed objects) lie in the range 50-100 (Wilson & Rood 1994). We therefore restrict the abundances used in the modelling to these ranges of values.

The canonical infall model
It is useful to define a specific infalling envelope model, with plausible physical parameters, which may be used to investigate some general properties of infall spectral line profiles. We use the above density and velocity profiles at a stage when a mass M ( / 2 has already accreted on to the central protostar, and a further .M ( / 2 of envelope gas is infalling towards it. If we choose an effective sound speed of a eff ¼ 0:35 km s 21 , then an infall radius of .3700 au is implied. We therefore use the following model relations for the systemic velocity and hydrogen number density profiles: n H2 2 Â 10 5 / cm 3 ¼ 0:35 r 3700 au 2 3 2 10:65 r 3700 au 20:64 ; and outside the infall radius the systematic velocity is set to zero. The density profile is given by We truncate the density profile at an outer radius of 10 000 au, which encloses a total mass of 2.75 M ( . This corresponds to the typical radius at which dense cores merge with the more diffuse gas in the ambient molecular cloud. The adopted microturbulent velocity dispersion is s tb ¼ 0:3 km s 21 ; which is assumed to be spatially constant. Although there is little empirical or physical justification for a flat turbulent velocity profile (see the above discussion), we make this assumption as a zeroth-order approximation, in the absence of any clear consensus as to the correct description of turbulence in infalling envelopes.
As discussed above, the kinetic temperature profile (assumed to be equal to the dust temperature profile) in the optically thin part of the envelope is expected to be fairly insensitive to the luminosity of the central source. We choose a canonical temperature profile of The normalization of the temperature profile was fixed using the dust temperature profile calculated by Kenyon et al. (1993), for a source luminosity of ,15 L ( . The varying kinetic temperature in the envelope is not strictly consistent with the assumption of constant a eff implied by the SIS model, where a 2 eff ¼ s 2 tb 1 kT/ m. However, our primary concern is not to investigate the SIS model solution in detail, but rather to investigate the essential features of line formation in infalling envelopes, whether the infall is described by the SIS model or not. We set both the CS and HCO 1 relative abundances to 5:0 Â 10 29 , and assume an ½H 12 CO 1 =½H 13 CO 1 isotopic ratio of 60. The set of molecular constants used in the radiative transfer modelling are listed in Table 1, where m denotes the permanent molecular dipole moment, and other symbols take their usual meanings. The data are taken from the JPL spectral line catalogue (Poynter & Pickett 1985).

The origin of the infall profile
We now use the radiative transfer code described above to simulate observed spectral line observations of the canonical infall model, and investigate how the predicted line profiles depend on a number of different model parameters. Fig. 3 shows the predicted HCO 1 , H 13 CO 1 and CS line profiles for the canonical infall model described above. The model cloud was assumed to lie at a distance of 200 pc, and the appropriate beam sizes for each transition were used when carrying out the beam convolution.
Double-peaked, blue-asymmetric line profiles are seen in all of the main isotopomer transitions, apart from the CSðJ ¼ 7 ! 6Þ line, which is nevertheless skewed slightly bluewards of the systemic velocity. The minima between the peaks of the doublepeaked profiles lie close to the systemic velocity. The HCO 1 ðJ ¼ 3 ! 2Þ line shows a stronger self-reversal than the HCO 1 ðJ ¼ 4 ! 3) line, and in both of these lines the selfreversal is stronger than the CSðJ ¼ 5 ! 4Þ line. The two singlepeaked H 13 CO 1 lines peak very close to the systemic velocity, coinciding with the minima in the main-line profiles. Low-level high-velocity wings are visible in all of the main-line profiles, and appear to be stronger in the higher rotational transitions of both molecules.
To illustrate how the asymmetries in line profiles arise, we show in Fig. 4 the separate contributions of the near and far portions of the envelope to the line profiles for the HCO 1 ðJ ¼ 4 ! 3Þ and CSðJ ¼ 5 ! 4Þ transitions. The nearer hemisphere (to the observer) of the envelope (in which the systematic velocity is redshifted) is seen to be mostly responsible for the redshifted peak in the line profile, as expected. However, the more surprising observation is that the nearer hemisphere also dominates the emission at velocities well into the blueshifted side of the line profile, and contributes significantly to the blueshifted peak.
Emission from the far hemisphere becomes dominant only in the most blueshifted part of the line. It is interesting to note that the contribution of the far hemisphere to the blue peak is no greater than the contribution of the near hemisphere to the red peak, and it is the extra contribution of the near hemisphere to the blueshifted emission which produces the blue asymmetry. This contrasts with the explanation usually given for the infall asymmetry, in which it is argued that the optical depth to the warm dense blueshifted gas just beyond the central protostar is less than the optical depth to the redshifted gas just in front of the protostar, due to the Doppler shift of the intervening infalling envelope. If this explanation were the correct one, then we would expect the blueshifted emission from the far hemisphere to be, by itself, much brighter than the redshifted emission from the near hemisphere, which contradicts our findings.    5 shows how the emission from the far half of the envelope is modified by the optical depth of the near half of the envelope on its way to the observer. Self-absorption clearly plays a very important role in shaping the emergent line profile. In Fig. 6 the contribution of the beam-averaged optical depths through both hemispheres of the envelope for both CSðJ ¼ 5 ! 4Þ and HCO 1 ðJ ¼ 4 ! 3Þ are plotted. The total optical depth of the HCO 1 ðJ ¼ 4 ! 3Þ line is much greater than the CSðJ ¼ 5 ! 4Þ line, which accounts for the stronger self-absorption features in the HCO 1 line. As expected, there is a separation in velocity of the peaks of the optical depth profiles from the near and far hemispheres, although this is rather small (,0.2 km s 21 ).
To shed further light on how the asymmetric infall profiles are formed, we show in Fig. 7 how the HCO 1 ðJ ¼ 4 ! 3Þ line profile    changes as successive 'layers' of the model envelope are included, from the far edge of the envelope forwards. For all the 'slice planes' in the rear of the envelope (left-hand panel of Fig. 7) the predicted profiles are single-peaked. This is explained by the fact that in the far half of the envelope, the higher excitation gas always lies in front of the lower excitation temperature gas. Fig. 8 shows the radial dependence of the systematic velocity, kinetic temperature, and HCO 1 ðJ ¼ 4 ! 3Þ and CSðJ ¼ 5 ! 4Þ excitation temperature in the model. The fact that the excitation temperature lies well below the kinetic temperature over most of the envelope is a demonstration of the inapplicability of the LTE approximation to this problem.
As successive layers of the front half of the envelope are included in the line profile calculation (right-hand panel of Fig. 7), there is a reduction in the intensity in the core of the line. This is due to the fact that in the near half of the cloud, lower excitation temperature gas lies in the foreground, and absorbs the emission from the higher excitation gas behind it. Nevertheless, the emission of this gas is also very important in determining the form of the emergent line profile. The emission from gas near the centre, in the near half of the envelope, is seen in Fig. 7 to produce an increase in the line intensity at strongly redshifted velocities. The contribution of this gas to the emission at lower velocities is obscured in the figure by the dominating effect of absorption towards the line centre, but Fig. 7 shows that this emission is significant even on the blueshifted side of the line. In the outermost region of the envelope, where the excitation temperature is lowest, the gas is predominantly absorbing, and simply deepens the self-reversal of the line profile. Fig. 9 illustrates how the predicted HCO 1 ðJ ¼ 4 ! 3Þ and CSðJ ¼ 5 ! 4Þ line profiles depend on the impact parameter of the line of sight. As expected, the strength of the lines diminishes with increasing impact parameter, as the beam samples increasingly less dense and lower temperature gas. The degree of asymmetry in the line profiles decreases with increasing impact parameter, partly because of the lower infall velocities at larger radii, and partly because of the lower column densities, and hence optical depth, along the line of sight. The lines also become narrower with increasing impact parameter, as a result of the decreasing infall velocity with radius.

Impact parameter and beam size
The distances to many nearby protostellar objects are often only known to an accuracy of ,50 per cent, so it is important to examine how sensitively the predicted line profiles depend on the assumed distance. There is, of course, a degeneracy between the assumed distance to the cloud and the beam size.
In Fig. 10 we plot HCO 1 ðJ ¼ 4 ! 3Þ and CSðJ ¼ 5 ! 4Þ line profiles for a number of different assumed distances, whilst keeping the beam sizes fixed. The peak line temperatures increase approximately as the inverse of the assumed distance (or beam size). The increase in the strength of the line wings depends approximately on the inverse square of the assumed distance, which suggests that the line wing flux originates from a region on the source much smaller than the beam size. This is explained by the fact that in the velocity profile we have adopted, the highest infall velocities lie at the smallest radii, and therefore the highvelocity wings of the line profile are formed in a very small region around the centre of the cloud. The different dependence of the line core and line wings on distance produces increasingly broad profiles with decreasing distance and/or beam size.
We conclude that distance uncertainties may contribute significantly to the overall uncertainty in the radiative transfer modelling of protostellar envelopes. Gregersen et al. (1997) analysed HCO 1 observations of a large sample of Class 0 sources using a model cloud at a single 'representative distance'. This  approach will not give reliable results when the scatter in the distances of the objects in the sample is large. Fig. 11 shows line profiles for a number of different infall radii, in which only the velocity profile is varied, as discussed in Section 3.1, and the density profile remains fixed as for the canonical model. For increasing values of the infall radius, larger fractions of the envelope take part in the infall motion, and the collapse speed at a given radius increases approximately as r 1=2

Infall velocity
inf . The line profiles become more strongly skewed towards blue velocities with increasing infall velocity, mainly as a result of the diminishing intensity of the redshifted peak. The actual velocities of the two peaks show remarkably little variation over the range of infall velocities considered. Predictably, the models with the highest infall velocities produce the strongest emission in the line wings. In actual observations this infall signature may often be obscured by emission from the outflow. Fig. 12 shows how the predicted line profiles vary with the assumed relative abundance of CS and HCO 1 in the envelope. The relative abundance of the tracer molecule largely determines the optical depth through the envelope in the observed line, and is therefore a very important parameter in deciding the overall appearance of the line profile. This is apparent from the figure, where the intensities and shapes of the line profiles are seen to vary significantly over the factor of 10 range of relative abundances covered. As expected, the lowest values of the relative abundance (and hence optical depth) produce the weakest infall profiles. The progression of line shapes with increasing optical depth is from: a single-peaked Gaussian; to a single-peaked blue-skewed profile; to a profile with a blueshifted peak and redshifted 'shoulder' (or 'red knee'); to a double-peaked profile with a stronger blueshifted peak. As the abundance and optical depth increase further, the absorption trough deepens, the line peaks become stronger and their separation increases. Fig. 13 shows how the line profiles depend on the magnitude of the turbulent velocity dispersion, s tb , when it is assumed to be uniform throughout the envelope. As the turbulent velocity dispersion increases, the most apparent effect on the line profiles is to increase the velocity separation between the two peaks, as a result of the   broadening of the absorption profile of the foreground envelope. The total integrated flux in the line tends to increase with the turbulent velocity dispersion. This is because the peak optical depth of the foreground gas is reduced as the molecules are spread over a wider range of velocities, and the similar reduction in the optical depth of the strongly emitting gas in the centre of the cloud allows more of the emission to 'escape'.

Turbulence
Comparing the HCO 1 ðJ ¼ 4 ! 3Þ with the CSðJ ¼ 5 ! 4Þ profiles, we see that the effect of changing the value of the turbulent velocity dispersion depends strongly on the peak optical depth in the transition. When the optical depth is small, increasing the turbulence tends to 'smear out' the double-peaked profile, whereas for very optically thick lines the double-peaked structure remains just as promiment as the velocity gap between the two peaks increases.
Note that this finding is in apparent contradiction with some previous models. For example, Myers et al. (1996) used a simple 'two-slab' geometry to model infalling protostars, and claimed that as a protostar evolves, the fraction of the envelope that is infalling increases (under the SIS inside-out collapse assumptions), and that this causes the gap between the two peaks to decrease, and the absorption trough to become narrower. Our modelling indicates that the amount of turbulence is what actually determines the gap between the peaks. The two scenarios could only be reconciled if turbulence were to decrease as a protostar evolves. This may seem unlikely, since the associated outflow is more likely to increase the turbulence as it evolves. However, a recent claim has been made (Jayawardhana, Hartmann and Calvet 2001, hereafter JHC) that in fact Class 0 sources may be formed preferentially in regions of higher turbulence than Class I sources. This is an interesting idea that is consistent with our modelling, although the other ideas of JHC that a significant fraction of a Class 0 envelope is not infalling, but rather static, is not consistent with our modelling of specific Class 0 sources (see below). Fig. 14 illustrates how the line profiles vary with the exponent in the assumed power-law relation between the turbulent velocity and the density. Guided by the discussion on turbulence in Section 3.2, we investigate turbulent velocity laws of the form s tb / r 1=4 and s tb / r 1=2 (as found in the numerical simulations of Vázquez-Semadeni et al. 1998), and s tb / r 21=2 (derived using Larson's Laws). For comparison we also plot the profile produced by the canonical model, which assumes a flat turbulent velocity profile.
The profiles are clearly quite sensitive to the form of the turbulent velocity law. The profiles obtained using the s tb / r 21=2 law (dashed lines) are appreciably narrower than the other profiles. This is particularly evident in the centres of the H 13 CO 1 ðJ ¼  4 ! 3Þ and CSðJ ¼ 7 ! 6Þ lines. This arises because the r 21/2 law assigns the lowest turbulent velocity dispersion to the hot, dense, inner regions of the cloud, where most of the line emission takes place. These profiles also tend to show the strongest emission close to the line centre, arising mainly from the lower optical depth of the outer envelope due to the enhanced velocity dispersion there. However, in the case of the H 13 CO 1 ðJ ¼ 4 ! 3Þ and possibly the CSðJ ¼ 7 ! 6Þ lines, this is due to the increased optical depth in the centre of the cloud, as these lines are not completely optically thick.
As the exponent s in the relation s tb / r s increases from 20.5 to 10.5, the line profiles become much more dominated by broad line wings. For the positive exponents, the largest velocity dispersions are located in the centre of the cloud, which causes the emission from this region to be distributed over a wide velocity range. The low turbulent and systematic velocities in the outer envelope confine the absorption of the gas in this region to line centre velocities, allowing the enhanced high-velocity emission from the core to remain virtually unobscured. Despite the change in overall appearance, the HCO 1 ðJ ¼ 4 ! 3Þ and CSðJ ¼ 5 ! 4Þ profiles retain the characteristic blue-asymmetric, double-peaked structure. The strong wings of these lines are similar to those normally associated with outflow emission. The more optically thin H 13 CO 1 ðJ ¼ 4 ! 3Þ and CSðJ ¼ 7 ! 6Þ lines may provide some means of distinguishing between these two possibilities. If the wings in the main-line profiles are caused by enhanced turbulence near the centre of the cloud, then the optically thin lines should also show very broad profiles. Conversely, if the wings are tracing a bipolar outflow, then the optically thin line profiles, which are expected to be dominated by the envelope, should be narrower. Fig. 15 shows how the predicted line profiles change when a solidbody rotation is superposed on the infall velocity field. This velocity field may not be physically realistic, since angular momentum conservation will tend to cause the angular velocity of the gas to increase as it falls inwards, even in the presence of magnetic braking (e.g. ). We use this example simply to illustrate some of the qualitative effects of rotation on line profiles. A proper treatment of differential rotation would require a three-dimensional radiative transfer analysis, which is computationally a very large step from the spherically symmetric problem presented here.

Solid-body rotation
The on-source line profiles are least affected by rotation, showing only a small decrease in the height of the peaks, and a very slight broadening of the line as a whole. For differential Keplerian ðV sin i , r 23=2 Þ rotation, the on-source profile might be more strongly affected, since the largest rotational velocities would then lie closest to the centre. The off-centre line profiles are more affected by solid-body rotation than the on-source profiles. As expected, the shifts in the centroid velocity of the line profiles follow the rotational velocity gradient. As well as shifting the centroid, the rotation also significantly distorts the shape of the off-centre line profiles. On the blueshifted side of the rotational velocity gradient, the rotation produces an enhanced blue asymmetry in the line profiles, in addition to the blue-asymmetry produced by the infall. At the positions on the redshifted side of the rotational velocity gradient, the blue infall asymmetry is completely reversed by the rotation.
Figs 16 and 17 show HCO 1 ðJ ¼ 4 ! 3Þ and CSðJ ¼ 5 ! 4Þ centroid velocity contour plots from the canonical infall model for projected rotational angular velocities of Vsin i ¼ 0 and 30 km s 21 pc 21 respectively. The non-rotating, infalling model produces circularly symmetric negative centroid velocity contours, which reach a minimum at the central position. The solid-body rotation tends to produce approximately parallel contours of centroid velocity, aligned with the axis of rotation.
However, the effect of the infall is to cause the contours to distort towards bluer velocities close to the centre of the cloud, producing the appearance a 'blue bulge' in the centroid velocity contour plot, which encroaches into the redshifted half of the rotational velocity gradient (Walker et al. 1994;Zhou 1995).
We can use measurements of the centroid velocity gradient across our spectral line maps of protostellar envelopes to place limits on the rotation rate of these envelopes. In Fig. 18 the  calculated HCO 1 ðJ ¼ 4 ! 3Þ and CSðJ ¼ 5 ! 4Þ centroid velocities are plotted, along axes perpendicular to and parallel to the rotation axis, at 1000-au intervals (i.e., 5-arcsec intervals for a distance of 200 pc). A straight-line fit to the centroid velocity along the perpendicular axis is shown in the left-hand panel of the figure. It is apparent that the CS centroid velocities lie behind the overall velocity gradient on both sides of this gradient. This is due to the larger assumed FWHM beam size in the calculation of the CSðJ ¼ 5 ! 4Þ profile (19.3 arcsec compared to 13.6 arcsec), which tends to increase the contribution of the brighter inner regions (which have lower rotational velocities) to the line profile. The centroid velocity along the axis parallel to the rotation axis in the right-hand panel clearly shows the 'blue bulge' signature discussed above, although in this case the magnitude of the associated centroid velocity shift is rather small.
The formal least-squares straight-line fit to the centroid velocity gradient perpendicular to the rotation axis is 23:56 :8 km s 21 pc 21 ; and the intercept is 20.018 km s 21 . The actual angular velocity used in the model was V sin i ¼ 30 km s 21 pc 21 , so the actual and measured values agree to within 1s. The slightly low measurement for the centroid velocity gradient is mainly due to the lagging behind of the CS centroid velocities mentioned above, as can be seen from the figure. Nevertheless, this example shows that the centroid velocity gradient gives a reasonable estimate of the actual projected velocity gradient across a cloud, at least for the case of solid-body rotation.

O B S E RVAT I O N S
The observations were carried out at the James Clerk Maxwell Telescope (JCMT), using receivers RxB3i ) and RxA2 , with the Digital Autocorrelation Spectrometer (DAS) backend. Details of the observing runs are given in Tables 2-4. Pointing and focus checks were performed roughly once per hour during each run. The focus was found to be very stable throughout, and the pointing accuracy was ,2 arcsec. Observations of standard sources were taken several times per night for each transition observed. We estimate the absolute calibration is accurate to within ,25 per cent, although the relative calibration between different observations is much better than this. Throughout the observations the DAS was configured with the optimum frequency resolution of 95 kHz per channel, giving a usable bandwidth of 125 MHz (equivalent to #0.1 km s 21 at the frequencies observed).
Most of the pointed observations were made using positionswitching. Frequency-switching was used only for some of the HCO 1 and H 13 CO 1 observations, and to check for emission-free reference positions during position switching. Beam-switching was used only occasionally, for sources which were known to have very compact emission. Comparison of identical observations using each of these observing modes gives excellent agreement in each case -see Matthews (1996) for discussion of observing modes. Both grid mapping and raster mapping modes were used to make spectral line maps. In some of the maps, the noise is greater towards the edge of the map than at the centre, due to a reduced integration time per map point at these positions. In the individual raster maps, integration times of 4-6 s per cell were used, as recommended by the JCMT Users Guide (Matthews 1996). Longer integration times per point were obtained by repeating the raster. Data from the JCMT data archive were also used to supplement our observations where possible, as indicated in Table 4.
Data reduction was carried out using the SPECX package (Padman 1990). Baseline calibration was performed by subtracting a fit through sections of baseline placed on either side of the spectral feature of interest, carefully chosen to avoid suppression of low-level wing emission. Frequency-switching tends to produce curved or sinusoidal baselines. Curved baselines were removed by subtracting a second-order polynomial fit in the vicinity of the line. This procedure was checked on several occasions by making observations of the same object in both frequency-switched and position-switched modes, giving very good agreement between the line profiles.
The spectra are calibrated in terms of the main-beam temperature scale, apart from the maps, where the corrected antenna temperature, T * A , is used. The integrated intensity and centroid velocity, position-velocity and channel maps were also produced using the standard SPECX map-making facilities. To improve the appearance of the contouring in these maps, the observed grids were first interpolated on to a grid with twice the spatial sampling rate. The relation between centroid velocity v c over a given velocity interval in terms of v i and T * Ai , the velocity and antenna temperature of the ith velocity channel, and the sum over the n ch velocity channels lying in the specified velocity interval (Narayanan, Walker & Buckley 1998) and s T A , the rms noise in a single channel, assumed to be constant for all the channels in the velocity window, is given by This equation is used to derive the error bars for the centroid velocity measurements presented below.
6 N G C 1 3 3 3 -I R A S 2

Previous observations
NGC 1333 is a reflection nebula associated with the L1450 dark cloud, which lies at a distance of ,220 pc (Černis 1990) in the Perseus molecular cloud complex (Sargent 1979). To the south of the reflection nebula lies a ,450-M ( molecular core, with a central density and temperature of 10 4 cm 23 and ,18 K respectively (Lada et al. 1974;Warin et al. 1996). This is a site of highly active low-and intermediate-mass star formation, evidenced by the large concentration of T Tauri and Herbig Ae/Be stars Lada, Alves & Lada 1996), Herbig -Haro objects (Bally, Devine & Reipurth 1996), and bipolar jets and outflows (Hodapp & Ladd 1995;Liseau, Sandell & Knee 1988) in and around the core. Loren (1976) originally suggested that the star formation in NGC 1333 is being triggered by the collision of two dense molecular clouds. More recently, Warin et al. (1996) have argued that the morphology of the region supports a sequential star formation scenario, where outflows produced by one generation of stars compress the surrounding gas, which triggers collapse and further star formation. Jennings et al. (1987) listed nine compact IRAS sources associated with the NGC 1333 cloud, five of which drive bipolar molecular outflows (Liseau et al. 1988). Of these, IRAS 2 and IRAS 4 have no associated on-source optical or near-infrared emission. A full submillimetre continuum survey of NGC 1333-IRAS 2 was carried out by Sandell et al. (1994) using the JCMT. A strong, compact peak was found at all wavelengths, with low-level emission at 800 mm extending north-west and south-east of the peak, in a 'flattened bar' morphology. Sandell et al. derived an envelope mass of ,0.8 M ( for the compact submillimetre source. The total far-infrared luminosity of IRAS 2 was estimated to be 17 L ( by Jennings et al. for a distance of 220 pc. An upper limit to its total luminosity of 26 L ( was found by Ward- , who quoted the ratio of bolometric to submillimetre luminosity to be L BOL / L SUBMM # 130, below the Class 0 threshold of 200 (André et al. 1993).
Two bipolar outflows have been associated with IRAS 2, observed in COðJ ¼ 3 ! 2Þ emission , several transitions of CS Langer, Castets & Lefloch 1996), and 2.12-mm shock-excited molecular hydrogen emission (Hodapp & Ladd 1995). The outflows are oriented nearly perpendicular to each other on the plane of the sky. However, since they typically affect only the extreme line wings, they do not alter the results of our modelling.   (Lovas 1992). The observing runs during which each transition was observed are indicated in the last column (cf. Table 2). Also indicated are transitions which were obtained from the JCMT data archive. The COðJ ¼ 3 ! 2Þ and CSðJ ¼ 7 ! 6Þ transitions were observed simultaneously, by selecting a front-end frequency which placed one transition in each sideband. (hereafter IRAS 2). Both the HCO 1 ðJ ¼ 4 ! 3Þ and CSðJ ¼ 5 ! 4Þ spectra show asymmetric double-peaked line profiles, with the blue peak stronger. The absorption dip in the HCO 1 ðJ ¼ 4 ! 3Þ line profile is particularly strong, and coincides in velocity with the peak of the rarer isotopomer lines. Also visible in the main line HCO 1 and CS spectra are broad, low-level wings, which are probably tracing the outflowing gas associated with this source. The peak of the C 18 OðJ ¼ 2 ! 1Þ line is shifted towards the blue relative to the H 13 CO 1 ðJ ¼ 4 ! 3Þ peak. Since H 13 CO 1 ðJ ¼ 4 ! 3Þ is expected to be the optically thinner of the two lines, we interpret this as showing that the C 18 O line is partially optically thick in this object. Several properties of the spectra agree with the predicted qualitative signatures of infall in protostellar envelopes (e.g. Zhou et al. 1993;Myers et al. 1995). For optically thin lines, radiative transfer models predict a symmetric single-peaked line profile, centred on the systemic velocity. Such a profile is shown by the H 13 CO 1 ðJ ¼ 4 ! 3Þ line, suggesting a systemic velocity of 7.7 km s 21 . The asymmetric double-peaked line profiles skewed towards blue velocities seen in both the HCO 1 ðJ ¼ 4 ! 3Þ and CSðJ ¼ 5 ! 4Þ spectra are predicted by radiative transfer models of infalling envelopes for optically thick lines. The velocity of the self-absorption minimum may coincide with the systemic velocity, or be redshifted with respect to it, depending on the optical depth of the line and the velocity and density structure of the envelope.

Molecule Transition
The C 18 OðJ ¼ 2 ! 1Þ line profile, which is slightly asymmetric, and blueshifted with respect to the systemic velocity, is also consistent with infall model predictions for a line with small but non-negligible optical depth. The only disagreements with the qualitative expectations of a pure infall model are the strong, highvelocity wings seen in several of the spectra, which we ascribe to the outflow (infall models do predict high-velocity wings, but at a lower level than observed here), and the very small redshift of the peak of the CSðJ ¼ 7 ! 6Þ line profile, with respect to the inferred systemic velocity. Mardones et al. (1997) obtained on-source spectra of IRAS 2 in CSðJ ¼ 2 ! 1Þ, N 2 H 1 ðJ F1F2 ¼ 1 01 ! 1 02 Þ and H 2 COðJ K21K1 ¼ 2 12 ! 1 11 Þ; using the IRAM 30-m telescope. The CS line is singlepeaked and approximately symmetrical, peaking at the systemic velocity (7.7 km s 21 ), and the N 2 H 1 line has an infall type asymmetrical profile. Within the infall scenario, this would suggest that the CS line is more optically thin than the N 2 H 1 line, contrary to expectations. The H 2 CO line profile observed by Mardones et al. has a complicated, multiple-peaked structure, with a very broad redshifted wing, suggesting that this line is sensitive to the outflow. The line core is symmetrical and skewed to the red, which contrasts with the blue-skewed profiles shown in Fig. 19. The complete set of spectral line observations at the position of this object do not, therefore, paint a wholly consistent picture. However, they are in reasonable agreement with the qualitative expectations for an infalling envelope, despite the fact that outflow emission is probably contributing significantly to the line centres in many of the observed transitions.
The spatial dependence of the observed line profiles may be able to shed light on this question. Grid maps of IRAS 2 in HCO 1 ðJ ¼ 4 ! 3Þ and CSðJ ¼ 5 ! 4Þ emission are presented in Figs 20 and 21 respectively. In both maps the emission is centrally peaked, and the strength of the self-absorption feature decreases away from the peak. The off-source spectra in the maps show an interesting variety of profiles. In the HCO 1 map there is a clear trend towards broader linewidths from south to north, which is also discernible in the CS map. The CS map spectra show broader and relatively stronger wing emission than the HCO 1 spectra, and may therefore be more sensitive to outflowing gas. In a detailed submillimetre spectral line study of the NGC 1333-IRAS 4 protobinary source, Blake et al. (1995) came to a similar conclusion, finding that the CS abundance is considerably enhanced in the outflow, and that HCO 1 appears to be less depleted than neutral species in the dense protostellar envelope.
North of the map centre, a dramatic reversal of the line asymmetry in the HCO 1 map is seen, and the CS spectra show weaker blue-asymmetries. One possible explanation for this is that the emission from the northern (redshifted) outflow lobe is enhancing the redshifted emission, perhaps due to an interaction with ambient gas at this position.
Centroid velocity maps calculated from the HCO 1 and CS grid  maps are shown in Fig. 22. The interpretation of centroid velocity maps requires some care, since the centroid velocity may be quite sensitive to the velocity range over which it is calculated. To illustrate this, we show HCO 1 ðJ ¼ 4 ! 3Þ and CSðJ ¼ 5 ! 4Þ centroid velocity plots calculated for two different velocity windows, one covering the entire line profile, including the high-velocity wings ð3:7-11:7 km s 21 Þ, and one covering mainly the line centre ð6:2-9:2 km s 21 Þ. All of the plots show a clear velocity gradient along an axis with a position angle between 08 and 208. The velocity gradient is in the same sense as the north-south bipolar outflow, suggesting a possible physical connection with it. The east -west outflow may be responsible for the small east -west velocity gradient apparent in the centroid velocity plots calculated over the wider velocity window.
The 'wide window' centroid velocity plots clearly show a much stronger velocity gradient than the 'narrow window' plots, particularly in the CSðJ ¼ 5 ! 4Þ maps. This is consistent with the idea that this gradient is caused by the north-south bipolar outflow, and also supports the suggestion that CS is a more sensitive outflow tracer than HCO 1 . The agreement in the direction of the centroid velocity gradient for both the narrow and wide velocity windows suggests that the outflow even affects the line profile shapes in the line centres. This adds a serious complication to the interpretation of the line centre profiles in terms of infall models, since it is then extremely difficult to separate reliably the contribution of the envelope and outflow to the line profile.
We now examine a possible alternative explanation for the velocity gradient seen in the centroid velocity maps, and the  reversal of line asymmetry seen in the HCO 1 ðJ ¼ 4 ! 3Þ maprotation of the protostellar envelope. Rotation will tend to produce a centroid velocity gradient across the envelope, and radiative transfer models which include both rotation and gas infall predict reversals of the infall asymmetry at certain off-centre positions (e.g. Zhou 1995), although the degree of reversal in the HCO 1 spectrum north of the map centre in Fig. 20 is larger than typically predicted by such models. Since the line wings of our CS and HCO 1 observations are clearly affected by the outflow, we concentrate on centroid velocities calculated over the line centre only. Fig. 23 shows linear fits to the combined CSðJ ¼ 5 ! 4Þ and HCO 1 ðJ ¼ 4 ! 3Þ centroid velocities calculated over the velocity window 6:7-8:7 km s 21 , along north-south and east -west axes through the source. The fit to the north-south centroid velocity gradient is 7:7ð^0:9Þ Â 10 23 km s 21 arcsec 21 . For a distance of 220 pc, this is equivalent to 7:2^0:8 km s 21 pc 21 . The corresponding numbers for the east -west velocity gradient are 28:0ð^0:8Þ Â 10 25 km s 21 arcsec 21 and 7:5ð^0:7Þ Â 10 22 km s 21 pc 21 . The interpretation of centroid velocities is complicated by optical depth effects, since both the transitions are optically thick in the line centre. Furthermore, infall motions in the envelope will tend to skew the profiles of optically thick lines towards bluer velocities.
One indication of this is the fact that the centroid velocities along the east -west axis are blueshifted by . 0.1 km s 21 , on average, with respect to the systemic velocity. The measured centroid velocity gradient for an optically thick line is therefore an unreliable measure of the actual velocity gradient in the object, and a proper treatment requires full 3D radiative transfer modelling, or the use of optically thin lines to eliminate optical depth effects. However, for the purposes of the present discussion we assume that our centroid velocity gradient measurements give a reasonable indication of the actual line-of-sight velocity gradient across the object due to rotation.
Using this assumption, a lower limit for the binding mass can be estimated by hypothesizing that the rotation is centrifugally balanced by gravity, and that the rotation axis lies in the plane of the sky. This gives a lower limit, because the actual rotation rate may be greater if the rotation axis is inclined out of the plane of the sky, and magnetic fields and pressure gradients may also contribute to the support of the envelope. Furthermore, if part of the envelope is infalling, as some of the observations suggest, then at some level the gravitational forces must dominate over the centrifugal support.
Assuming spherical symmetry, the mass M interior to a radius r is then given by where r is the distance from the centre of rotation, v(r) is the lineof-sight velocity at r (relative to the systemic velocity), V(r) is the angular velocity of the rotation at r, G is the gravitational constant, and i is the inclination angle of the rotation axis from the plane of the sky. Using equation (31) with i ¼ 0, we find a lower limit of 0.4 M ( for the enclosed mass inside a 30-arcsec (6600-au) radius. Sandell et al. (1994) estimated a total dust and gas mass in the IRAS 2 protostellar envelope of 0.79 M ( , and to this must be added the mass of the central protostar. This analysis therefore does not contradict the hypothesis that the line-centre centroid velocity gradient is tracing rotation about an east -west axis projected on to the plane of the sky. If this were correct, the most likely model would be that the east -west bipolar outflow is driven by the Class 0 source, and the protostellar envelope is rotating about this outflow axis. The older north-south outflow would then have to be driven by a separate, more evolved, and as yet unobserved object. This seems unlikely, and we consider it more likely that the velocity gradient is tracing the north-south outflow rather than rotation.

S E R P E N S S M M 1 -4 7.1 Previous observations
The Serpens molecular cloud is a low-to intermediate-mass starforming region which has been the subject of intensive observational study (e.g. Eiroa 1991;Eiroa & Casali 1992;White, Casali & Eiroa 1995, hereafter WCE;Hogerheijde et al. 1998;Testi & Sargent 1998;Davis et al. 2000;Testi et al. 2000;Williams & Myers 2000). The distance to the region has been the source of some debate, with estimates ranging from 250 pc (Chavarria et al. 1988) to 750 pc (Zhang et al. 1988), based on spectroscopic parallaxes of stars believed to be physically related to the cloud. More recent studies (de Lara, Chavarria & Lopez-Molina 1991;Straižys, Č ernis & Bartašiūtė 1996) have identified mis-classifications of the spectral types of some of the stars used in the previous distance determinations, and suggest an actual distance of 300^30 pc. Eiroa & Casali (1992) carried out a comprehensive near-infrared survey of the region and positively identified 51 young stellar objects associated with the cloud. They estimated a lower limit of .450 pc 23 for the average stellar density over a region of .0.5-pc diameter, and a star formation efficiency in the range 8-28 per cent, suggesting that a bound cluster may be forming. The peak of the star formation activity is centred close to the Serpens reflection nebula (SRN), which is a diffuse optical and near-infrared nebula (Gómez de Castro, Eiroa & Lenzen 1988;Eiroa & Casali 1992) powered by the luminous pre-main-sequence star SVS2 (Strom, Vrba & Strom 1976).
Maps in ammonia (Ungerechts & Güsten 1984), C 17 O (WCE), CS (McMullin et al. 1994b) and far-infrared emission (Zhang et al. 1988; show the gas and dust in the central region to be concentrated in two cores, aligned approximately north-west -south-east. Estimates of the gas density and temperature in the cores lie in the range 10 4 -10 5 cm 23 and 20-30 K respectively, with higher temperatures and densities towards some of the embedded sources. The ,70-M ( north-western core contains the well-studied far-infrared source Serpens FIRS1 (Harvey, Wilking & Joy 1984), which is associated with a highly collimated, high-velocity (,300 km s 21 ) radio continuum jet (Curiel et al. 1993(Curiel et al. , 1996, and a linear near-infrared feature, approximately aligned with the north-western lobe of the radio jet (Eiroa & Casali 1989). This source was labelled SMM1 (Casali, Eiroa & Duncan 1993, hereafter CED), and is the strongest of the submillimetre continuum sources in the cloud.
The more massive south-eastern core includes the reflection nebula, and contains a small cluster of compact submillimetre continuum sources (CED), two of which (SMM2 and SMM4) have no near-infrared counterparts. Barsony (1997) used a maximumcorrelation algorithm to enhance the spatial resolution of the IRAS maps of this region, and derived upper limits for the far-infrared luminosities of the submillimetre sources. SMM1-4, and S68N (McMullin et al. 1994b) were identified as Class 0 sources. Hurt, Barsony & Wootten (1996) have observed each of these sources in the high-density molecular tracer H 2 CO. In the H 2 COðJ K21K1 ¼ 3 03 ! 2 02 Þ transition, four out of the five sources showed asymmetrical line profiles suggestive of infall. The remaining source, SMM1, showed a symmetrical double-peaked profile. Gregersen et al. (1997) included Serpens SMM1-4 in their HCO 1 survey of Class 0 sources, and found strong infall signatures towards SMM2 and SMM4.
The COðJ ¼ 2 ! 1Þ outflow map of WCE shows possible bipolar outflow associations with SMM1 (south-east -north-west) and SMM4 (approximately north-south), and perhaps SMM2 (south-east -north-west). Establishing the presence or absence of outflows driven by the other Serpens objects is hindered by the considerable confusion from overlapping outflow lobes. Fig. 24 shows the HCO 1 ðJ ¼ 4 ! 3Þ and H 13 CO 1 ðJ ¼ 4 ! 3Þ spectra observed towards the Serpens Class 0 sources SMM 1 -4. Of the four sources, SMM2 and SMM4 show line profiles most suggestive of infall. SMM4 in particular shows the classic signature of a double-peaked self-absorbed main-line profile skewed heavily towards the blue, with an absorption minimum redshifted relative to the peak of the optically thin isotopic line. The SMM1 main-line profile shows a deep absorption minimum, close to the peak velocity of the isotopic line, and a very mild blue asymmetry. Strong high-velocity wings are also visible, which are probably tracing the energetic bipolar jet and outflow driven by this source. The HCO 1 ðJ ¼ 4 ! 3Þ spectrum observed towards SMM3 does not show any sign of self-absorption, but a strong blueshifted wing, probably tracing an outflow, is clearly seen. In the following we concentrate primarily on the sources SMM2 and SMM4, as the best protostellar infall candidates in the Serpens cloud. Fig. 25(a) shows the on-source spectra observed towards Serpens SMM2. Blue asymmetric line profiles are seen in most of the transitions, although CSðJ ¼ 7 ! 6Þ emission was not detected above the ,0.2-K noise of the observation. There is a remarkable similarity between line profiles of the main-line HCO 1 transitions and their isotopic counterparts, although the signal-to-noise ratio in the H 13 CO 1 ðJ ¼ 4 ! 3Þ line profile may be too low to attach much weight to this. Fig. 26(a) shows a grid map of Serpens SMM2 in the HCO 1 ðJ ¼ 4 ! 3Þ transition. This map was obtained from the JCMT data archive, and the observing runs during which these data were taken are described in WCE. The HCO 1 ðJ ¼ 4 ! 3Þ emission is seen to peak north-west of the continuum peak. However, the bright line profiles observed towards the north-western corner of the map include a contribution from the outer envelope of SMM4, which lies at an offset of (225,150), and shows extended bright emission in HCO 1 ðJ ¼ 4 ! 3Þ -see Fig. 27. This may in part explain the similar velocity structure of the isotopic and main-line profiles discussed above, if the emission arises from overlapping clumps with different velocities along the line of sight. There is no reliable way to separate the overlapping contributions to the line profiles, and we therefore do not discuss this source any further.

New data
The spectra observed towards the position of the Serpens SMM4 submillimetre continuum peak are shown in Fig. 25(b). The optically thin lines -H 13 CO 1 ðJ ¼ 3 ! 2Þ and H 13 CO 1 ðJ ¼ 4 ! 3Þ -peak at a velocity of 7.95 km s 21 . The C 18 OðJ ¼ 3 ! 2Þ line peaks at 7.7 km s 21 , suggesting that there is a non-negligible optical depth in this line, as we have found previously. There is good qualitative agreement between the line profiles shown and typical radiative  transfer model predictions for an infalling envelope (see the discussion above of the NGC 1333-IRAS 2 on-source spectra).
Both the main-line HCO 1 transitions show the asymmetric double-peaked line profile skewed towards blue velocities. We also note that the COðJ ¼ 3 ! 2Þ line also shows this type of profile. However, CO line profiles are not as robust infall indicators as HCO1, since they suffer greater confusion from foreground material and outflowing gas. The isotopic HCO 1 line profiles shown are also in qualitative agreement with model predictions, being single-peaked, and peaking between the velocity of the absorption minimum and the blueshifted peak of the main-line profile.
Further qualitative agreement between model predictions and observations is shown by three isotopic CO lines. Models predict that as the optical depth of the observed transition increases, the line profiles progress from single-peaked and symmetric, to singlepeaked and skewed to the blue, to double-peaked with the blue peak brighter than the red (Myers et al. 1995). This is just the progression observed in the CO line profiles as the isotopic abundance increases. The CS line profiles represent the intermediate optical depth stage, being single-peaked and skewed to the blue. Taken as a whole, therefore, the set of line profiles shown in Fig. 25(b) are remarkably consistent with the qualitative expectations for an infalling envelope. Fig. 26(b) shows a grid map of the SMM4 envelope in the HCO 1 ðJ ¼ 4 ! 3Þ transition. The strongest emission is found at the origin and at the grid position immediately north of the origin. Blue-skewed line profiles are evident in most of the spectra shown, although double-peaked line profiles are found only in the central and western map columns. Comparison of the spectra in the southeast and north-west corners of the map reveals a shift in the peak velocity of about 1 km s 21 towards the blue along the south-eastnorth-west direction -see Fig. 27(a). A centroid velocity map calculated from this grid map is shown in Fig. 27(b), from which a south-east -north-west centroid velocity gradient of ,30 km s 21 pc 21 is measured. Also apparent in this plot is a 'bluebulge' encroaching on to the redshifted side of the velocity gradient from the blueshifted side. As discussed above, this is a predicted signature of infall in the presence of rotation.
A larger area HCO 1 ðJ ¼ 4 ! 3Þ raster map of SMM4 was also made, and is shown in the left-hand panel of Fig. 28 as an integrated intensity map. The HCO 1 ðJ ¼ 4 ! 3Þ emission is seen to be strongly enhanced in the region of SMM4 and is well resolved in our beam (the transition from dark to light contours in the figure marks the half-maximum level). The emission is approximately centred around the offset position (0, 15) on the map. Towards this point the integrated intensity flattens off into a plateau, and there is a shallow minimum at the central point itself, surrounded by horseshoe-shaped peak. The lower contours reveal an extension of the peak towards the north-east, and a flattened linear feature oriented north-west -south-east.
The right-hand panel of Fig. 28 shows the centroid velocity map calculated from the raster map. The south-east -north-west velocity gradient is again very apparent. Comparison of integrated intensity and centroid velocity maps shows that the region of brightest integrated HCO 1 ðJ ¼ 4 ! 3Þ emission coincides with the blueshifted part of the velocity gradient, which makes the symmetric appearance of the plot about the line rest velocity all the more remarkable. The axis of the steepest velocity gradient ðP:A: . 2408^58Þ is close to (but not exactly coincident with) the axis of the elongated feature in the integrated intensity map ðP:A: . 2508^38Þ.
The steepest velocity gradient is 24:1ð^0:15Þ Â 10 22 km s 21 arcsec 21 . which implies a gradient of 28^1 km s 21 pc 21 , for a distance of 300 pc. The centroid velocities along this axis clearly shift towards bluer velocities at positions near to the central peak. This 'blue-bulge' signature, parallel to the rotation axis of a rotating, infalling protostellar envelope, was first predicted by the radiative transfer analysis of Walker et al. (1994).
As before, we estimate the lower limit on the mass assuming the observed velocity gradient is due to rotation and the rotation axis lies in the plane of the sky. We consider the mass enclosed within the limits of the steep velocity gradient between position offsets of 220 and 110 arcsec from the HCO 1 peak, and obtain a lower limit of ,1.9 M ( (for a distance of 300 pc).  estimated a gas mass of 3 M ( surrounding SMM4, derived from dust continuum emission, which, along with the mass of any embedded protostars, is sufficient to provide the centripetal force required to bind the rotation, as long as the rotation axis is close to the plane of the sky. Nevertheless, if our analysis is correct, centrifugal effects are likely to have a very significant influence on the equilibrium and dynamics of the envelope gas. The elongation seen in the integrated HCO 1 ðJ ¼ 4 ! 3Þ map, which is aligned approximately perpendicular to the inferred rotation axis, may well in fact be caused by centrifugal flattening of the outer envelope (e.g. Pudritz & Norman 1986).
The alternative possible explanation for the observed velocity gradient in the HCO 1 ðJ ¼ 4 ! 3Þ maps could be outflowing gas. Gregersen et al. (1997) detected a velocity gradient in the same sense and along the same axis in their HCO 1 ðJ ¼ 3 ! 2Þ map of SMM4, but interpreted it in terms of outflow rather than rotation. SMM4 has previously been associated with a north-south bipolar outflow (WCE), but not aligned with the gradient we observe. Examination of the channel maps in our own data likewise reveals no evidence for an outflow in this direction. Thus we believe our rotation hypothesis to be the more likely explanation. We approached the modelling by first searching for a fit to the HCO 1 ðJ ¼ 4 ! 3Þ profile, choosing realistic physical parameters. We aimed to achieve consistency in our model fits with previously determined properties of the objects, such as the envelope mass and ambient gas temperature. Once a suitable fit was found, the same model was used to find the best fit to the H 13 CO 1 ðJ ¼ 4 ! 3Þ profile, but in this case the only free parameter was the ½H 12 CO 1 =½H 13 CO 1 ratio, which we required to lie in the range 50-100. Depending on Figure 28. HCO 1 ðJ ¼ 4 ! 3Þ integrated intensity (left) and centroid velocity (right) maps of Serpens SMM4, over the velocity interval 4-12 km s 21 . The peak contour of the integrated intensity map is 18.6 K km s 21 , and the contour spacing is 0.6 K km s 21 (.1s ). In the centroid velocity map, the minimum and maximum contours levels are at 6.8 and 8.8 km s 21 respectively, and the contour spacing is 0.1 km s 21 . The rest velocity of the source (.7.95 km s 21 ) lies in the transition interval between the dark and light contours. The cross indicates the effective FWHM beam size. the quality of this fit, the model was either rejected, or tested further against the CS line profiles.
We attached less weight to the CS profiles in the model fitting, since CS tends to be more affected by outflow emission than HCO 1 (e.g. Blake et al. 1995). Langer et al. (1996) mapped the entire NGC 1333 core region in the CSðJ ¼ 5 ! 4Þ, ðJ ¼ 3 ! 2Þ and ðJ ¼ 2 ! 1Þ transitions using the IRAM telescope, and found that the CS emission extended over a very wide area. The CS profiles may therefore be significantly affected by emission and absorption from the gas in the foreground core.
The parameters used in the model fit are listed in Table 5. The outer radius of the cloud was set at 10 000 au, and an innermost shell radius of 15 au was used. A distance to the cloud of 220 pc was assumed (Č ernis 1990). The total envelope mass within the outer radius of the cloud inferred from the molecular hydrogen density profile is 3.7 M ( . The values of r inf used in the infall velocity and density profiles were 6000 and 900 au respectively, and a eff ¼ 0:29 km s 21 was used in both cases. No acceptable fit could be found which used identical values of r inf for both profiles. In other words, we could not find a solution where the density and velocity profiles taken together are consistent with the singular isothermal sphere model. The fit to the HCO 1 ðJ ¼ 4 ! 3Þ line profile is very good between 21 and 11 km s 21 (relative to the systemic velocity of 7.75 km s 21 ). At larger velocities in both the red and blue wings of the line, the model significantly underestimates the emission. The most likely origin for the excess high-velocity emission is from one or both of the bipolar outflows associated with this source -the clear break in the slope of the observed line profile at .^1.3 km s 21 also suggests a separate origin for the highvelocity emission. We have not tried to estimate how much the outflow contributes to the line profile at lower velocities, so this is another source of uncertainty in the model.
The H 13 CO 1 ðJ ¼ 4 ! 3Þ model fit agrees well with the peak velocity and intensity in the observed line, although the model profile has slightly overestimated the linewidth. The width of the H 13 CO 1 ðJ ¼ 4 ! 3Þ line in the model is determined mostly by the systematic and turbulent velocities near the centre of the cloud, since the warm dense gas in this region has the strongest intrinsic emission, and for an optically thin tracer this emission is unattenuated by the outer envelope on its way to the observer. The asymmetry in the main-line profile is sensitive mainly to the systematic and turbulent velocities at larger radii. If the infall velocity at large radii is chosen to produce agreement with the asymmetry in the HCO 1 ðJ ¼ 4 ! 3Þ profile, then the inferred velocity at small radii always produces too much broadening in the H 13 CO 1 ðJ ¼ 4 ! 3Þ line, as a result of the v / r 21=2 velocity law.
The largest discrepancy between the observations and the model is in the CSðJ ¼ 5 ! 4Þ line, where the model predicts a substantially higher intensity in the line core than is observed. Better fits to the overall intensity of the CS line (but not the line shape) could be found by reducing the CS relative abundance, although this then underestimated the intensity of the CSðJ ¼ 7 ! 6Þ line. This might be explained if the CS relative abundance  (FWHM) 0:60 km s 21 X HCO + 6:0 Â 10 29 X CS 6:0 Â 10 29

12
C/ 13 C 9 0 Figure 30. Model fits using the same model parameters as in Fig. 29, but with the infall velocity set to zero inside a radius of 800 au, and a factor of 4 reduction in the CS relative abundance. may vary with radius, but we did not pursue this. The value of the CS relative abundance used in the model fits was simply optimized to give the best agreement with the CSðJ ¼ 7 ! 6Þ line. Fig. 30 shows an alternative fit to the HCO 1 ðJ ¼ 4 ! 3Þ and H 13 CO 1 ðJ ¼ 4 ! 3Þ line profiles obtained by setting the infall velocity to zero inside a radius of 800 au, keeping all other parameters the same. The main-line profile is affected remarkably little by this dramatic alteration to the inner velocity field (although the fit is slightly worse), which underlines the point that each observed transition constrains only a limited domain of parameter space in the model envelope. The CS line profiles (in which the CS relative abundance has been lowered by a factor of 4 relative to the original model fit) are still not well fitted by this model. In both Figs 29 and 30 the predicted H 13 CO 1 ðJ ¼ 4 ! 3Þ and CSðJ ¼ 7 ! 6Þ line profiles are almost identical to each other, and we could not produce simultaneous fits to the velocity widths of both observed profiles. Since CS emission shows a greater tendency for contamination by the outflow, we are inclined to attribute the larger velocity width on the CSðJ ¼ 7 ! 6Þ line to the outflow, and therefore we tend to prefer the former fit described above.
In Fig. 31 we compare the off-centre HCO 1 ðJ ¼ 4 ! 3Þ line profiles with the model predictions. Because of the assumption of spherical symmetry, the predicted line profiles depend only on the impact parameter of the line of sight relative to the origin, and not on its direction. The decrease in line intensity with increasing impact parameter is reasonably well matched by the model. However, the line shapes at the off-centre positions show poorer agreement. With the exception of the spectrum at the (0,115) position, the off-source line profiles are all skewed towards the blue. The model profiles are much more symmetric, as a result of the fall-off in infall velocity at large radius, which is again a consequence of the v / r 21=2 infall velocity law. It is very difficult to accommodate the strong high-velocity emission seen in several of the off-source positions -most prominently at (0,115) and (115,215) -within a plausible infall model, particularly given the highly asymmetric distribution of the high-velocity emission in both position and velocity, which cannot be fitted by a spherically symmetric model.
We explored the possibility above that the centroid velocity gradient measured across the NGC 1333-IRAS 2 envelope using HCO 1 ðJ ¼ 4 ! 3Þ and CSðJ ¼ 5 ! 4Þ spectral line maps could be explained by rotation, and concluded that although rotation could not be ruled out, outflow was the more likely explanation. On the basis of the results in Fig. 31, we believe that the north-south bipolar outflow, and not rotation, is responsible for the centroid velocity gradient.
The mass of the central protostars can be approximately estimated from the infall velocity profile used in the model fit, applying the free-fall equation M * ¼ rv 2 / 2G to the inner part of the profile. Using this relation we obtain an estimate of . 0.3 M ( . However, this must be regarded as somewhat uncertain. The density and velocity profiles adopted in this model are inconsistent with the SIS model, despite individually having the SIS model forms. The infall radius used in the velocity profile is ,7 times greater than the infall radius used in the density profile, indicating that the infall is more well developed than a SIS model envelope with a similar density profile. Hydrodynamical (e.g. Foster & Chevalier 1993) and magnetohydrodynamical (e.g. Ciolek & Mouschovias 1994) simulations of the evolution of dense cloud cores during the pre-stellar phase predict that by the time the density profile approximates closely to a singular distribution, substantial infall velocities have already developed. Notwithstanding the uncertainties in our model fit, our results support this picture of a non-static initial condition for collapse.

Serpens SMM4
Fig. 32 shows the model fits to the CS and HCO 1 on-source spectra observed towards Serpens SMM4. The approach to the modelling was similar to that used for NGC 1333-IRAS 2, starting with the HCO 1 ðJ ¼ 4 ! 3Þ line, and then making adjustments to produce the optimal agreement with the remaining lines. To produce a fit to the HCO 1 ðJ ¼ 4 ! 3Þ line required the infall radius in the SIS model form of the velocity profile to be set well beyond the outer radius of the cloud, i.e., the entire protostellar envelope must be infalling.
If the infall radius is set inside the outer radius of the cloud, then absorption from the static gas beyond this radius tends to produce an absorption dip close to the systemic velocity, in disagreement with the observed profile. As the SIS model velocity profile is inappropriate in this case, we used a simple r 21/2 power law, which is the asymptotic form of the inner part of the SIS model profile, to model the velocity profile instead.
The parameters used in the model are given in Table 6. The outer radius of the model envelope was taken as 10 000 au, and we used an innermost shell radius of 62 au. The assumed distance to the source was 300 pc. The total envelope mass implied by the model is ,5.1 M ( , roughly consistent with the estimate of ,3 M ( derived from the IRAS fluxes by .
As with NGC 1333-IRAS 2, the infall velocities in the model fit   (FWHM) 0.80 km s 21 X HCO + 3:5 Â 10 29 X CS 3:5 Â 10 29 12 C/ 13 C 9 0 are much greater than would be predicted by the SIS model given the density profile. From the infall velocity profile we estimate the mass of the central protostellar system to be M * ¼ rv 2 / 2G ¼ 2:1 M ( : The temperature profile is normalized at a considerably higher level than in the IRAS 2 fit, which would suggest that SMM4 is considerably more luminous. A high accretion luminosity for this object is also suggested by the large central protostellar mass and mass accretion rate (.6 Â 10 25 M ( yr 21 Þ.  estimated the luminosity of SMM4 as 9 L ( , after deconvolving the HIRES IRAS images of the Serpens region to obtain the far-infrared fluxes. Jennings et al. (1987) estimated the luminosity of IRAS 2 as 17 L ( , also using IRAS fluxes, so there are clearly large uncertainties on these estimates. The parameters used in our model fits are inconsistent with these luminosities, and for our model to be correct, either the SMM4 luminosity calculated from the IRAS data, or our assumed protostellar radius, would have to be a significant underestimate. The predicted H 13 CO 1 linewidths are considerably broader than the observations, as we found with NGC 1333-IRAS 2, although Figure 33. Comparison of the off-centre HCO 1 ðJ ¼ 4 !3Þ line profiles observed towards Serpens SMM4 (solid lines) with the model predictions (dashed lines). The offsets (in arcsec) from the IRAS 2 position are given in the title of each plot. In the model spectra we have included a solid-body rotation of 28 km s 21 pc 21 , about an axis passing through the centre of the cloud at a position angle of +408. The centre of the cloud was assumed to lie at the offset (0,+3.5).
the peak intensities are in very good agreement. Again, we suggest this may be indicating deviations of the infall velocity profile below the r 21/2 free-fall profile in the centre of the cloud. We found above that the SMM4 envelope shows strong evidence for rotation ðV sin i . 28 km s 21 pc 21 Þ, and suggested that the flattening seen in the spatial distribution of integrated HCO 1 ðJ ¼ 4 ! 3Þ emission may be caused by centrifugal support of the outer envelope. Thus centrifugal braking may retard the infall velocity near the centre of the cloud.
The calculated CSðJ ¼ 5 ! 4Þ and CSðJ ¼ 7 ! 6Þ line profiles match the observed line intensities reasonably well, but the peak velocities of the calculated profiles are less redshifted than the observations. The enhanced blueshifted wings seen in both of the observed CS spectra strongly suggest that outflow contributes to the emission in these lines. As discussed above, CS is believed to be enhanced relative to HCO 1 in bipolar outflows (Blake et al. 1995). Fig. 33 shows the fits to the off-centre HCO 1 ðJ ¼ 4 ! 3Þ line profiles towards SMM4 from our best-fitting model. We included a solid-body rotation of 28 km s 21 pc 21 in the model (assuming a rotation axis in the plane of the sky), consistent with the above centroid velocity gradient across SMM4. As the HCO 1 ðJ ¼ 4 ! 3Þ spectra at the (0,0) and (0,17) positions have approximately equal intensities, we assumed that the centre of the cloud lies at the offset (0,13.5). Although the model profiles somewhat underestimate the observed lines, we are encouraged by the overall agreement of the trends in the observed lineshapes across the envelope. It is likely that a better fit to the line intensities requires a shallower density profile than predicted by the SIS model.
In contrast to IRAS 2, we find little evidence for outflow emission in the off-centre line profiles, and essentially all of the high-velocity emission can be accounted for by the infall. We consider these results to provide convincing evidence that both infall and rotation are taking place in the Serpens SMM4 envelope. The magnitude of the rotational velocity is large enough for centrifugal effects to have a significant influence on the gas dynamics.

SUMMARY
In this paper we have discussed the physical constraints which can be used to limit the number and range of parameters needed to define the spherically symmetric models in the Stenholm radiative transfer code. We have studied the dependence of the predicted line profiles on a number of model parameters. Increasing the infall velocities in the envelope tends to produce stronger asymmetry between the blue-and redshifted peaks of double-peaked line profiles, without strongly affecting the separation between the peaks. The value of the turbulent velocity dispersion is the most significant factor in determining the velocity separation between the peaks. We carried out a preliminary investigation into the effect of adopting different power-law dependences of the turbulent velocity width on the gas density, and found the predicted line profiles to be extremely sensitive to the value of the power-law exponent.
We have studied the effect on the predicted line profiles of superimposing a solid-body rotation on to the infall velocity field. The effect of the solid-body rotation on the on-source doublepeaked line profiles was to 'smooth out' the double-peaked structure. The centroid velocities and the line shapes of the offcentre profiles were found to be significantly altered by the solidbody rotation. On the redshifted side of the rotational velocity gradient, a reversal of the asymmetry in the line profiles was seen. We measured the centroid velocity gradient predicted by our model for a solid-body rotating envelope, which was found to slightly underestimate the actual rotational velocity gradient used in the model.
We selected the objects NGC 1333-IRAS 2 and Serpens SMM4 for detailed modelling using the Stenholm radiative transfer code, confining our search to the SIS model forms for the density and velocity profiles in the modelling, although we did not require these profiles to use the same SIS model parameters in the models. Our best model fits to both the NGC 1333-IRAS 2 and Serpens SMM4 data used velocity profiles which were more strongly developed than the SIS model would predict, given the fit to the density profile. This is consistent with theoretical predictions that at the instant a central protostar is formed, the gas in the envelope has already accelerated to substantial velocities (e.g. Foster & Chevalier 1993;JHC), which contrasts with the static initial state assumed by the SIS model, and is one of the principal distinguishing features between the SIS model and other collapse models.
Our fits to the H 13 CO 1 line profile observations of both objects slightly overestimated the profile widths. We could produce better fits to these lines by lowering the infall velocity in the centre of the cloud, although this affected the fits to the main-line profile slightly. The infall velocity may therefore be retarded from the r 21/2 free-fall relation towards the centre of the cloud, possibly by centrifugal braking. In the case of IRAS 2, the discrepancies between the model predictions and the observations for the offsource profiles could not possibly be accommodated within a plausible model of infall and rotation, and we attribute them to the outflow.
The model fits to the off-centre HCO 1 ðJ ¼ 4 ! 3Þ line profiles of Serpens SMM4 included the solid-body rotation inferred from the HCO 1 ðJ ¼ 4 ! 3Þ centroid velocity gradient in the calculation. The predicted and observed profiles showed good agreement in the overall line shape, although the model slightly underestimated the line intensities. We thus conclude that Serpens SMM4 is an infalling and rotating protostar.