4 Cloud parameters from the nonlocal model

4.1 Line fitting by SimLine

We performed non-local radiative transfer simulations using the line radiative transfer code SimLine introduced in detail in Appendix B. SimLine is a FORTRAN code to compute the profiles of molecular rotational lines in spherically symmetric clouds with arbitrary density, temperature and velocity distribution. It consists of two parts: the self-consistent solution of the balance equations for all level populations and energy densities at all radial points and the computation of the emergent line profiles observed by a telescope with finite beam width and arbitrary offset. The optical depths in the lines may vary from minus a few, corresponding to weak masking, to several thousand.

Already in the spherically symmetric description of a core we face a large number of parameters. For all quantities (hydrogen density, kinetic temperature, velocity dispersion, molecular abundances) a radial function has to be found. Regarding the limited amount of information available from the three transitions, this leaves many options open. We decided to assume simple power-law radial functions and a central region with constant parameters for all quantities in the core simulations. This reduces the number of parameters to two (central value and radial exponent) for each gas property, plus the outer and inner radius.

The parameter fit procedure used the multidimensional downhill simplex algorithm from Press et al. (1992). Although it is not the most efficient way in terms of convergence speed it turned out to be very robust in all situations considered. Because a downhill simplex code does not necessarily find the global minimum of the $\chi ^2$ function we performed for each core several runs with randomly chosen initial simplex covering a large part of the physically reasonable parameter space. For all clouds we made at least 30 runs to get a rough idea of the topology of the $\chi ^2$function. For cores like W33 this turned out to be sufficient since only one large minimum showed up which was found in almost half of the runs. The other extreme is S106 where we needed almost 1000 runs to be sure that we found the global minimum. Here, the $\chi ^2$ function was quite complex with numerous local minima. Future improvements of the fit procedure should include more sophisticated algorithms like simulated-annealing approaches.

The noise in the line profiles produces some graininess of the $\chi ^2$function when directly fitting the measured profiles. This results in a very slow convergence of the algorithm close to the minimum. A considerable acceleration can be obtained by not fitting the measured noisy data but a smooth approximation to them. The line profiles were represented by a superposition of a Gaussian and a Lorentzian profile. This allows a good characterisation of all measured profiles including the reproduction of asymmetric profiles, self-absorption dips and line wings.

Taking the three measured transitions and their spatial extent as the quantities to be reproduced by a $\chi ^2$ fit we find that we are able to derive at most six parameters to a reasonable accuracy. Fits with seven or eight free parameters, although still slightly improving the numerical $\chi ^2$ value, do not produce any significant changes above the noise limit. This is comparable to the results by Young et al. (1998) fitting the full position velocity map of a particular core in a single transition. There is however a number of additional parameters where the model can derive certain limits. Hence, we have to decide first which parameters should be fitted directly, which parameters may be constrained by preventing reasonable fits when outside a certain range, and which parameters can be guessed independently from a physical line of reasoning. This is discussed in detail for all quantities in the following subsections.

Unfortunately, it is impossible to give an comprehensive error estimate for the parameters derived from the simulations. This would need a description of the six-dimensional surface in the parameter space within which none of the observational error bars is exceeded. According to the complex topology of the $\chi ^2$ function it is not possible to give an easy description for the six-dimensional valley around the global minimum or its boundaries.

As a simple alternative we performed only one-dimensional variations to get a rough estimate of the maximum error in the parameters that we must expect. After the convergence of the $\chi ^2$ fit, each of the fit parameters was varied up and down until one of the computed lines deviated by the assumed maximum observational error of 25% from the measured lines. The central values of the different functions were varied independently. When changing the inner radius, the central values of the gas parameters were adjusted to keep the functions in the power-law region unchanged. When varying exponents the corresponding central values were corrected in such a way that the parameters at the density of $2\times10^{6}$ cm^-3 remained constant. By the selection of this density as the fix point when changing the slope, we scan about the maximum possible range of the exponents. This provides a conservative estimate of the maximum error.

Figure 1: Line profiles observed at the central position of S255 in four CS transitions together with best-fit models. The dashed line represents the microturbulent approximation, the dotted line stands for clumping in velocity space with cell sizes of 0.01 pc and the solid line includes the additional effect of density clumping where the combined abundance is reduced by a factor 100.

4.2 The turbulence description

In spite of the large number of free parameters in the cloud models it is impossible to fit the line profiles with a smooth density and velocity distribution. The typical self-absorbed profiles known for all microturbulent codes (White 1977) appear in this case. Thus we have to take into account the effects of internal clumping and turbulence leading to a more realistic picture and preventing strong self-absorption.

SimLine treats turbulence and clumping in a local statistical approximation following Martin et al. (1984) (see Appendix B.2). The cloud material is subdivided into small coherent units (clumps) with only thermal internal velocity dispersion. The relative motion of many units then provides the full velocity profile. Martin et al. (1984) showed that the effective optical depth of such an ensemble can be considerably reduced compared to the microturbulent approximation. The relative reduction of the total optical depth depends on the optical depth, i.e. of size the coherent units.

This description does not need any assumptions on the nature of the turbulence creating the internal cloud structure. The reduction occurs in the same way whether clumps are units of the same velocity in a medium of constant density, representing the behaviour of incompressible turbulence, or whether they are density enhancements in a thin inter-clump medium. The different nature of these two scenarios has to be taken into account, however, when computing the excitation. Both the average gas density, providing a measure for the column density and thus the line intensity, and the local density, providing the collisional coupling to the gas, enter the balance equations. In case of coherent units in velocity space both densities agree.

For density clumps we use the additional simplification to treat the cloud locally as a two-component medium neglecting the contribution of the inter-clump medium to the radiative transfer. Then the collisional excitation is provided by the density in the clumps and the column density is provided by the average density. The ratio between the two quantities reflects the filling factor of the volume occupied by dense clumps. From the mathematical point of view this is equivalent to the treatment of the cloud as a homogenous medium where the clumping occurs only in velocity space and the abundance of the radiating molecules is reduced. Hence, it is impossible to separate the influence of the filling factor from that of the molecular abundance so that the line fit provides only a combined quantity which we will denote combined abundance in the following. Beside the modification of the molecular abundance the statistical turbulence description introduces the size of the coherent units as an additional parameter.

Figure 1 demonstrates the influence of the turbulence. It shows the central line profiles of the best fits to the S255 observations using different assumptions on the turbulent nature. We selected S255 here because we can exploit the advantage of additional data for the CS 10-9 transition measured by Plume et al. (1997) providing more constraints on the cloud model. In all models the core is optically thin in CS 7-6 and 10-9, so that the reduction of the optical depth by the turbulent clumping produces only minor changes in these lines. In CS 2-1 the differences are, however, most obvious. In the microturbulent description and for the incompressible turbulence we find self-absorbed line profiles. A reasonable fit to the observations is only possible using the turbulence model including clumps with enhanced density. The reduction of the effective optical depth of the cloud by increasing the optical depth of the coherent units produces narrower lines providing a better fit to the observed line profiles.

We find that the parameter fits do not provide accurate values for the clump size and the combined abundances reflecting the volume filling factor. They only constrain an interval of possible values. We obtained good fits to the observations for the full range of clump sizes between about 0.005 and 0.05 pc. This size scale is confirmed by independent determinations of clump sizes by high-resolution observations in some of our cores. For NGC 2024 Mezger et al. (1992) determined a radius of dense condensations around 0.015 pc and the interferometric studies by Wiesemeyer et al. (1997) showed values between 0.005 and 0.01 pc. In W3 Tieftrunk et al. (1998) observed compact clumps with a size of 0.02 pc and in OMC-2 Chini et al. (1997) found dust condensations with radii between 0.01 and 0.05 pc.

For a better comparison to turbulence theory we prefer to specify the cell size in terms of the turbulent correlation length $l_{\rm corr}$ which should be on the order of 0.1 pc (Miesch et al. 1994; Goodman et al. 1998). The size of the units which are coherent with respect to the line radiative transfer is smaller by the ratio of the thermal line width to the total velocity dispersion. Thus the size range found corresponds to correlation lengths between 0.04 and 0.4 pc. In the following computations we use the intermediate value of 0.1 pc as correlation length for all clouds.

Regarding the combined abundances we find two classes of objects. The majority of cores, including the example of S255, is well fitted by values between 10^-11 and 10^-10, whereas a second class, consisting of W3, Serpens, and S106 needs values between 10^-10 and 10^-9 for a good fit. Because we do not know the molecular abundances there is no way to translate these values directly into clump filling factors. Assuming the CS abundances of $1.3\times10^{-9}$ to $1.3\times10^{-8}$ obtained by Hatchell et al. (1998) for several star-forming cores, the first class corresponds to filling factors around 0.01, whereas the CS abundance of $4\times10^{-10}$from Plume et al. (1997) corresponds to a filling factor of 0.1. The filling factors in the second group are ten times higher accordingly. In all following computations we have used a combined abundance factor of $3\times 10^{-11}$for the first and $3\times 10^{-10}$ for the second group. In the translation to cloud masses we will use the intermediate CS abundance of $4\times10^{-9}$.

As an additional parameter quantifying turbulence in a one-dimensional cloud model we have to take a radial variation of the turbulent velocity distribution into account to explain the observed size-line width and size-density relations (Larson 1969). Exponents of this radial dependence between about 0.1 and 0.7 are typically discussed (see e.g. Goodman et al. 1998). We left the width and the exponent of the turbulent velocity distribution as free parameters and obtained exponents between 0.15 and 0.65.  

4.3 The density structure

Any physically reasonable cloud model should include a spatial dependence of the gas parameters. Collapse simulations might provide reasonable model assumptions for the density structure. Bodenheimer & Sweigart (1968) and Shu (1977) have shown that an isothermal sphere evolves into a power law density profile $n_{\rm H_2} \propto r^{-\alpha\rm _n}$ with $\alpha\rm _n=2.0$. Homologous collapse simulations provide an exponent $\alpha\rm _n=3.0$ (Dickel & Auer 1994) and the free-fall collapse discussed by Welch et al. (1987) results in a density exponent $\alpha\rm _n=1.5$. The inside-out collapse model (Shu 1997) combined two regions of different exponents and recent more sophisticated collapse simulations (see e.g. Basu & Mouschovias 1995) show more complex density structures with an average exponent $\alpha\rm _n$ between 1.5 and 1.7. Dust observations of the density profile of protostellar cores show evidence both for cores with a typical r^-2 profile and for cores with flat density structure and a sharp outer edge (André et al. 1999). Thus we expect exponents between about 1.5 and 2.0 in our power-law density model which is bound by an outer cut-off and a central constant region. Although this simple model may not reflect the whole complexity of the density profile we can hardly derive any more information from the limited observations available. We left the central density and the density exponent as free parameters to be fitted. The resulting exponents span the relatively wide range between 1.1 and 2.2.

In the fit of the radii confining the power-law density profile we face two problems. We cannot distinguish changes of the model parameters on the smallest size scales where even the highest observed transition is thermalised because of the high density. Thus we can only set an upper limit $R_{\rm cent}^{\rm max}$to the radius of the central region where a transition from the power law behaviour in the envelope to constant parameters might occur. Only in four clouds - Serpens, Mon R2, $\rho$ Oph A, and the b-component of W49A - the strength of the CS 7-6 line sets an upper limit to the density so that we can derive the inner radius directly from the observations. In all parameter fits, we left the central radius as a free parameter. After the fit we increased the radius until the maximum deviation in one of the lines reached 5%. This provides a reasonable upper limit to the inner radius in all cases where the line profiles are independent of the density structure below this limit and gives only small modifications in the four cases where the inner radius was already well constrained by the fit.

The outer radius of the cloud is also quite uncertain. We can easily provide a value for the extent of gas at densities above about $3\times10^5$ cm^-3 based on the spatial extent of the CS 2-1 emission and the line profiles. But it is not possible to derive a reliable value for the extent of low density gas. We can only give a lower limit to the outer radius and thus the mass of the massive cores considered. As the density exponents are typically shallower than -2, a majority of mass could be present beyond this radial limit at low densities invisible in CS. We excluded the outer radius from the parameter fit using a sufficiently large value for all clouds and performed later a separate run reducing the radius to find the minimum outer radius $R_{\rm cloud}^{\rm min}$ in a way equivalent to the maximisation of the inner radius described above.  

4.4 The temperature structure

The temperature distribution of massive cores is still a matter of debate (cf. Garay & Lizano 1999). During early phases of cloud collapse the temperature should remain constant as long as the core remains optically thin. Deviations are to be expected, however, as protostellar sources are formed in most of our cores, leading to an internal heating of the cloud. Moreover in thin outer regions external heating can be important. Based on several observational results Scoville & Kwan (1976) set up a spherical cloud model with a warm inner region resulting in an temperature profile $T_{\rm kin} \propto r^{\alpha\rm _T}$ with $\alpha\rm _T\approx -0.4$.

Hence, we should also derive the core temperature and the temperature exponent from the radiative transfer model. In a first run we have investigated the influence of a temperature gradient in S255 when fitting only the central line profiles. We compared the best fit models to the S255 observations using either a constant cloud temperature, a temperature decaying with the exponent -0.4, or the temperature exponent as a free parameter. The latter case provided a best fitting exponent of -0.12. All three fits showed an excellent agreement in terms of the central line profiles falling almost exactly on the solid curves in Fig. 1. Thus, it is impossible to favour a certain exponent from the least squares fit of the line profiles only.

In Table 6 we see the resulting values for the other core parameters. The average temperature given here is computed as the mass-weighted average up to $R_{\rm cloud}^{\rm min}$. As main difference we find a kind of compensation between temperature and density exponent. The sum of both exponents is kept approximately constant to fulfil the constraints given by the line ratios. Temperature and density gradient act in a similar manner, leading to higher excitation in regions which are either denser or warmer. The variation of the density gradient, however, changes the extent of the emission in the CS 2-1 transition (the change is much smaller in the higher transitions). Consequently, we can constrain the temperature exponent when taking the observed size of the source into account. The beam convolved FWHM of the CS 2-1 emission in S255 falls between 1.5' and 2.1' clearly excluding models with steep temperature gradients.

 

 

Table 6: Resulting cloud parameters for S255 assuming different exponents for the radial temperature dependence and fitting only the central line profiles.

$\alpha\rm _T$	$T_{\rm c}$	$\langle T \rangle$	$\alpha\rm _n$	$\langle n_{\rm H_2} \rangle$	$R_{\rm cloud}^{\min}$	$M_{\rm cloud}$	$\mbox{{\it FWHM}}^{\rm a}$
	[K]	[K]		[cm-3]	[pc]	[$M_{\odot}$]	[']
0.0	47	47	-1.5	3300	1.1	1200	1.8
-0.12	73	52	-1.3	2500	1.3	1500	2.2
-0.4	220	52	-1.0	1000	2.5	7000	3.6

$^{{\rm a}}$ CS 2-1 transition.

In the fit procedure applied to all cores we have thus included the fit of the spatial extent. In this way we can constrain the temperature exponent but we are not able to derive exact values. In the parameter fit we always start with an isothermal cloud and change the temperature gradient in steps of 0.1 until the model provides a simultaneous good fit to the central line profiles and the spatial FWHM of the emission. It turns out that the observations of all clouds except W33, W3 and Mon R2 can be fitted by an isothermal model. This indicates that a large part of the gas mass in the clouds is characterised by a uniform kinetic temperature. This isothermal behaviour, however, does not support the assumption of constant excitation temperatures like in escape probability model. The excitation temperatures show local variations which are steeper than the density gradient. In the derivation of the kinetic temperature structure we have to keep in mind, however, that we are not very sensitive to the exact value of the kinetic temperature gradient. Observations of higher transitions like CS 14-13 would be needed to obtain reliable values including full error estimates.

Figure 2: Observed line profiles and best fit model to the observations of NGC 2024 using a CS to C34S abundance ratio of 13.

4.5 Other parameters

Beside the temperature gradient we can expect the formation of compact H II regions in the centre of a star forming core. Using the parameters of the central H II region derived by Dickel & Auer (1994) for W49 (R=0.2 pc, $n_{\rm e}=2.6\times10^4$ cm^-3, $T_{\rm e}=10^4$ K) we have compared the resulting CS lines when either including or neglecting the H II region in the radiative transfer computations (see Appendix B.3). We find that the influence of the H II region is negligible in this example for the four CS transitions considered here.

In general H II regions have two dominant effects on the line profiles. The integrated brightness changes the molecular excitation throughout the cloud and is visible as continuum emission. Moreover, the molecular material in front of an H II region, appears partially in absorption. A large H II region with electron densities being a factor 10 or more higher than in the example above, would thus result in distinct changes in the lines. Depending on the configuration and velocity structure, the CS lines may appear in absorption or with P Cygni profiles. Moreover, a strong continuum emission would be observed at all frequencies considered. However, bright H II regions with high electron densities are unlikely to be that extended (Wood & Churchwell 1989). For ultracompact H II regions the change of the line profiles by absorption is negligible due to the small angular size of the H II region so that the molecular excitation is the remaining effect. From the lack of a bright continuum underlying the lines we can, however, exclude configurations with a bright H II region here. Weaker compact or ultracompact H II regions - although possible - were not included in the fit computations as they would only influence material in their close environment which cannot be resolved in this study.

Regarding the chemical evolution of massive cores one should also expect a variation of the molecular abundances of CS and C³⁴S. However, our present knowledge is still insufficient to guess reliable values here (cf. Bergin & Langer 1997). As discussed in Sect. 4.2, it is also not possible to fit the abundance independently from the clump filling factor, so that we adopted here a constant CS abundance of $4\times10^{-9}$ relative to H₂ ignoring any radial variation of the abundances.

An additional test is only possible for NGC 2024 where C³⁴S is sufficiently bright so that we could include it in the fits. Figure 2 shows the best fit model to all six available line profiles. Equivalent to the results from the escape probability model we get the best match for a relative molecular abundance X(CS)/X(C³⁴S $) \approx 13$. All six lines are simultaneously fitted.

None of the cores except W49A show clearly asymmetric profiles as a signature of systematic internal velocities. The weak wings seen in a few of the other lines at low resolution are insufficient to derive any collapse or outflow model. Thus we have fitted all cores with a static model.

  

Table 7: Direct fit parameters and their likely ranges obtained in the $\chi ^2$ minimisation.

$$\begin{tabular}{lr@{\,}lr@{\,}lllr@{\,}lrll} \hline \multicol... ...&[31:40] & 0.0 & 1.2\,[0.8:1.8] & 0.48 \\ \hline \end{tabular}$$

4.6 Comparing the massive cores

Table 7 lists the resulting best fit parameters for all cores with their error intervals. As discussed in Sect. 4.3 the intervals for the radii are in most cases only limited at one end. The central density refers to the value at the radius $R_{\rm cloud}^{\rm min}$.

Comparing the different cores we find that the inner and outer radii are mainly determined by the selectional bias from the observability as one massive core. The most distant cores are only detectable with the KOSMA telescope if they are relatively large whereas at small distances only small cores are unresolved.

The minimum central density fitting the CS lines varies between 10⁶ cm^-3 for W49A^(b) and $8\times10^7$ cm^-3 for NGC 2024. The high value derived for NGC 2024 is due to the availability of the CS 10-9 observations tracing higher densities. With the three CS lines measured for most cores, only the density range below about 10⁷ can be reliably traced, so that we have to take the central density for all clouds except W49A^(b), Serpens, Mon R2, and $\rho$ Oph A as lower limits. The density exponent covers the range between -1.1 and -2.2 where the majority of clouds shows values around -1.6 corresponding to large-scale collapse models (Sect. 4.3).

Although we don't have a sample where we can expect to set up statistically significant correlations we can interpret some general relations. It turns out that the parameters are not completely independent of each other. Whereas the majority of clouds shows an average temperature between 20 and 50K, the three clouds which needed a larger combined abundance factor corresponding to a higher clump filling factor in the turbulence description (W3, S106, Serpens) also tend to require relatively high temperatures. There is a clear correlation between the cloud temperature and the turbulent line width indicating that heating and turbulent driving might have a related cause. The cores with a significant temperature exponent also show a relatively steep density exponent whereas a steep density exponent itself does not necessarily require a temperature exponent. These internal relations should be explained from the physical nature of the clouds.

 

 

Table 8: Resulting cloud parameters for S255 using either three or four transitions in the fit.

$n_{\rm l}$	$R_{\rm cent}^{\rm max}$	$R_{\rm cloud}^{\rm min}$	$n_{\rm cent}$	$\alpha\rm _n$	$T_{\rm cent}$	$\langle n_{\rm H_2} \rangle$	$M_{\rm cloud}$
	[pc]	[pc]	[cm-3]		[K]	[cm-3]	[$M_{\odot}$]
4	0.06	1.1	$1.5\times10^7$	-1.5	47	3600	1000
3	0.11	1.1	$5.3\times10^6$	-1.4	49	3700	1000

There are some peculiarities concerning four massive cores. The fit of the second component of W49A needs a large inner region with constant parameters or a very shallow decay of the density profile. Here, the main information that we get from the line profiles is the central density whereas we can hardly constrain the density exponent in the power-law range. Regarding the error bars of the parameters for Serpens we see that the the density structure of this relatively nearby cloud is well constrained by its observed size whereas the temperature structure is relatively poorly determined. In $\rho$ Oph it was not possible to fit the line profiles simultaneously with the size of the smallest core resolved in the observations. This can be explained by the clumpy structure of the whole region where the excitation of the core at the central position cannot be treated separately from the other clumps. Thus we have used for the core model a size which is six times larger than the smallest resolved clump and contains most of the strong emission observed. To fit the profiles and spatial extent of Serpens and S106 we have to assume a relatively large combined abundance (see Sect. 4.2) and a low hydrogen density in the clumps. This peculiarity could be removed when assuming some foreground CS that increases the apparent size of the massive core in CS 2-1 but does not contribute to the excitation in the core.