text:
# Application of time-dependent density-functional theory to electron-ion coupling in ethylene

## I Introduction

The time-dependent density functional theory (TDDFT) offers a compromise between economy and accuracy for the calculation of electronic excitations in atoms, molecules and condensed matter. The electronic excitation energies are typically reproduced by a few tenths of an eV for states with large oscillator strengths, and their strengths are typically reproduced within 25% accuracy . We are now interested in seeing how much further the theory can be applied at this level of accuracy and, in particular, how well it works in describing the coupling of electronic excitations to the nuclear degrees of freedom.

In a previous publication, we examined this question taking the benzene molecule for a case study. There we found that the TDDFT worked rather well at describing the nominally forbidden transitions, reproducing oscillator strength over 3 orders of magnitude with errors at the 30% level. Besides giving transition strength to the symmetry-forbidden transitions, the electron-nucleus coupling induces a vibrational splitting of each electronic line. While the detailed description requires the Franck-Condon factors of the individual vibrational states, the overall widths can be calculated with much simpler theory. For benzene, we found that the overall widths associated with the vibrational splittings were reproduced within 25% accuracy.

The ethylene is an even simpler molecule whose absorption spectrum has been subject to much theoretical study by the methods of quantum chemistry. The observed photoabsorption spectrum in ethylene shows a strong, structured peak in the region of 7-8.5 eV excitation that may be largely attributed to the $`\pi `$-$`\pi ^{}`$ transition. The studies of ref.  and  present detailed calculations of the vibronic couplings, but they do not make a direct comparison to experiment. We also note the calculations by Ben-Nun and Martinez  in which the subsequent nuclear motion is treated quantum mechanically. The goal of the present article is less ambitious: we aim to see how well the density functional theory (DFT) works for electron-nucleus coupling at short times, before the system responds dynamically. We will focus on the $`\pi `$-$`\pi ^{}`$ transition, which is of general interest for the optical properties of conjugated carbon systems.

The basic properties of the transition, besides its mean energy, are its width (about 1 eV FWHM) and its oscillator strength, measured to be $`f=0.42`$ . A feature of particular interest is the long progression of strength on the low-energy side of the peak, extending down to about 6 eV. We shall also apply the TDDFT to this feature and directly compare the theoretical and the experimental strength in the region of 6-7 eV.

## II Calculational Method

We use the DFT to calculate the ground potential energy surface (PES) as well as the excited PES. In the present work, we use a simple density functional based on the local density approximation (LDA) .

Denoting the nuclear coordinates by $`Q`$, the ground PES $`E_{gs}(Q)`$ is obtained in the static theory by minimizing the energy functional of the orbital variables. The optimized $`\varphi _i`$ satisfy the Kohn-Sham equations

$$\frac{^2}{2m}\varphi _i+\frac{\delta }{\delta n}\varphi _i=ϵ_i\varphi _i.$$
(1)
Here $``$ is the energy functional including Coulomb interactions but not including the electrons’ orbital kinetic energy.

The excited states are calculated by solving the linearized equations of the TDDFT. These equations are very close to the RPA equations, differing only by the presence of an exchange-correlation term in the interaction. For each ground state orbital $`\varphi _i`$, there are two excited wave function, $`\varphi _i^+and\varphi _i^{}`$. The equations they satisfy are

$$\frac{^2}{2m}\varphi _i^\pm +\frac{\delta }{\delta n}\varphi _i^\pm ϵ_i\varphi _i^\pm +\frac{\delta ^2}{\delta n^2}\delta n\varphi _i=(ϵ_i\pm \omega )\varphi _i^\pm ,$$

with the normalization condition, $`\varphi _i^+|\varphi _i^+\varphi _i^{}|\varphi _i^{}=1`$. The physical quantities extracted from the solutions are the eigenvalues $`\omega `$, representing vertical excitation energies, and the transition densities $`\delta n=_i\varphi _i(\varphi _i^++\varphi _i^{})`$ which represent the matrix elements of the density operator between the ground and the excited states. The excited PES is given by the sum of the ground PES and the excitation energy $`\omega `$,

$$E_{ex}(Q)=E_{gs}(Q)+\omega (Q).$$
(2)
To describe the oscillator strength distribution, we treat nuclear motion classically except for the zero-point motion on the ground PES. This is called the reflection approximation . Also, we neglect couplings between different vibrational degrees of freedom and treat the fluctuations in various coordinates $`Q`$ independently. The reflection formula for the strength function is

$$\frac{df}{dE}=f(Q)\left(\frac{dE_{ex}}{dQ}\right)^1|\mathrm{\Psi }(Q)|^2,$$
(3)
where the coordinate $`Q`$ in this expression is obtained by inverting the equation $`E=E_{ex}(Q)`$. In the formula, $`\mathrm{\Psi }(Q)`$ denotes the vibrational nuclear wave function in the ground PES, normalized as $`|\mathrm{\Psi }|^2𝑑Q=1`$. Also $`f(Q)`$ denotes the oscillator strength at a fixed $`Q`$ between the ground and the excited state. In ref. , we used an even simpler approximation, namely taking the vibrations to be harmaonic and treating other factors in eq. (3) as constant. Then the width is Gaussian with an energy variance $`\sigma _E`$ given by

$$\sigma _E=\frac{dE_{ex}}{dQ}|_0\sigma _Q,$$
(4)
where $`\sigma _Q`$ is the variance of $`Q`$ in the ground state vibrational wave function. For a Gaussian profile, the variance is related to the width $`\mathrm{\Gamma }`$ by $`\mathrm{\Gamma }=2.3\sigma _E`$.

Our computational method is quite different from usual quantum chemistry methods in that we employ a three-dimensional Cartesian grid representation for the orbital wave functions. Only valence electrons are explicitly included in the functional. More tightly bound electrons are taken into account by pseudopotentials constructed with a commonly used prescription . Further details on the computation may be found in our previous publications, e.g. . The important numerical parameters in the computations are the mesh spacing $`\mathrm{\Delta }x`$ and the shape and size of the volume on which the grid is defined. In our previous work on conjugated carbon systems, we used a grid spacing of $`\mathrm{\Delta }x=0.3`$ Å. With that mesh, orbital energies are converged to better than 0.1 eV. One can use the same coarse mesh to compute the coupling to the vibrational degrees of freedom, provided the ground PES is independently known. In ref.  we only considered small amplitudes of nuclear motion, obtaining the necessary information from empirical data on the vibrations. Here we want to consider larger amplitudes of motion, beyond the harmonic regime. This requires a more accurate calculation of the ground PES, obtainable with a finer mesh.

## III Results

### A Ground PES

We first consider the ground PES. To obtain converged results for equilibrium nuclear positions and vibrational frequencies, we found it necessary to use mesh spacing of $`\mathrm{\Delta }x=0.2`$ Å. The orbital wave functions are defined on the grid points that lie within a sphere. We found that a sphere of radius $`R=5`$ Å is adequate to obtain converged ground orbital energies. Some properties related to the ground PES are shown in Table I. The C-C equilibrium distance is reproduced to 1% accuracy. The curvature of the potential energy surface for the two important modes, the C-C stretch and the torsion, control the respective vibrational frequencies. As may be seen from the Table, the empirical vibrational frequencies  are reproduced to about 5% accuracy. We also quote the results of the ab initio calculations of ref.  and , which show more or less the same accuracy. We conclude that the quality of the DFT results are encouraging to proceed with the calculation of anharmonic and large amplitude fluctuation effects.

### B Vertical Excitations

Next we turn to the TDDFT excitation energies at the equilibrium geometry (call the Franck-Condon point, FCP). The TDDFT predictions with the LDA energy functional are shown in Table II. The $`\pi `$-$`\pi ^{}`$ transition appears at 7.5 eV excitation energy with an oscillator strength of $`f=0.30`$. This energy is about 0.2 eV lower that the center of gravity of the observed peak. Also, comparing to the measured strength of 0.42 in this region, we infer that additional electronic transitions are likely to be present. Indeed, the TDDFT produces other states in the same region of excitation. The three additional states have a Rydberg character with diffuse orbitals. Only the lowest of these, which is $`s`$-like and has no angular nodes at large distance, is optically active. It has a transition strength of $`f=0.08`$. Adding this to the $`\pi `$-$`\pi ^{}`$ strength, the total strength is only 0.04 units abelow the observed value. We note that our oscillator strength are quite close to those obtained by quantum chemistry methods, eg. ref. . In principle, the LDA cannot give reliable excitation energies for the Rydberg states, due to an incorrect asymtotic potential. However, the present results are within 0.3 eV of the more refined theoretical values in the last column. It might be possible to improve the energies by using a GGA density functional as in ref. , but for the present study we will use the LDA and confine our attention to the $`\pi `$-$`\pi ^{}`$ transition.

### C Coupling to the nuclear coordinates

The nuclear coordinates that have the largest effect on the $`\pi `$-$`\pi ^{}`$ transition are the C-C stretch coordinate and the torsion angle coordinate. We shall consider both of these in the sections below. The CH2 wagging coordinates also play a role, first in mixing the $`\pi `$-$`\pi ^{}`$ excitation with Rydberg excitations , and then at large torsion angles providing a direct path from the excited PES to the ground PES . Since we ignore the Rydberg states and do not treat the the multidimensional character of the nuclear coordinate space at all, we shall neglect the wagging and other modes that couple indirectly.

#### 1 C-C stretch motion

Historically, the first candidate for the source of the width of the $`\pi `$-$`\pi ^{}`$ transition was the coupling to the C-C stretch coordinate. Calculating the excited PES along the stretch coordinate according to eq. (2). we find the results shown in Fig. 1. The excited PES has a minimum at $`R_{CC}=1.5`$ Å with an energy 7.0 eV. The slope of the PES at the FCP is given by $`dE_{ex}/dQ=8`$ eV/Å. This is all we need to estimate the width using eq. (4). Taking $`\sigma _Q=0.04`$ Å from empirical vibrational frequency, we find a variance of $`\sigma _E0.3`$ eV for the absorption peak. The corresponding width is $`\mathrm{\Gamma }=0.7`$ eV, roughly 2/3 of the observed peak width. We expect that the widths from fluctuation of different coordinates add in quadrature. Thus, the width due to fluctuation in the C-C distance must be augmented by a similar width from other coordinates.

The observed peak has an apparent asymmetry, falling off more steeply on the low-energy side. In principle the asymmetry can be calculated using the full expression eq. (3). We first solve the Schroedinger equation on the ground PES to get the wave function $`\mathrm{\Psi }(Q)`$, and then apply eq. (3) directly. The result is shown in Fig. 2. Indeed, a significant asymmetry is predicted, reflecting the anharmonicity in the ground PES and the curvature of the excited PES. However the width is still too small. We mention that the calculation of ref.  also predict an asymmetric peak. Their width includes all the vibrational modes, and comes out somewhat larger than the observed width.

#### 2 Coupling to torsional motion

Experimentally, the strong absorption peak at 7-8 eV has a progression of strength at low energies, seen down to nearly 6 eV . The C-C stretch PES doesn’t drop to that low an energy, leaving the torsional mode as the most likely source of the strength. Combining the inertia $`I_\theta `$ (with $`\mathrm{}^2/I_\theta 0.0048`$ eV/r2) of the torsional motion with empirical vibrational frequency $`\omega _{torsion}`$, the zero-point fluctuation variance in the torsional angle is given by $`\theta ^2^{1/2}=\sqrt{\mathrm{}/2I_\theta \omega _{torsion}}8^{}`$. This is rather soft, and suggests that coupling to this coordinate could be significant. Unfortunately, the $`\pi `$-$`\pi ^{}`$ excitation mixes with the Rydberg states at finite torsional angles, splitting the original $`\pi `$-$`\pi ^{}`$ transition into several components. Since the Rydberg excitation energies are not accurate in the LDA, we choose to ignore the splitting of the $`\pi `$-$`\pi ^{}`$ state, taking the diabatic state at the mean excitation energy associated with an excitation along the $`x`$ coordinate. The diabatic excited PES is shown in Fig. 3. Using it in eq. (3), we obtain the strength distribution shown in Fig. 4. Qualitatively the measured fall-off is reproduced, but in detail the predictions are not as accurate as we had found for the forbidden transitions in benzene. In the region of 6.5-7.0 eV, the theoretical strength is too high, by about a factor of three. The agreement is much better at lower energies, down to 6.1 eV, the lowest energy measured in ref. . However, the theory should be more reliable the closer to the FCP, so the low energy agreement is probably a fortuitous cancellation of opposite-sign errors.

A likely cause of the discrepancy just below 7 eV is the error in the TDDFT excitation energy. In Fig. 2, we saw that the strength around the main peak in the LDA calculation is somewhat lower than the observed one. This is within a typical error of the TDDFT for excitation energy. A overall shift of 0.2 eV in the excited state PES would improve the description of strength function both in the peak region and in the region as far down as 6.5 eV. Shifting the excited state PES in that way, on the other hand, would destroy the present agreement below 6.5 eV. However, this large-angle region of PES is only accessible by a deep tunneling of the nuclear motion. It might be that the naive reflection treatment becomes inaccurate. For example, it might be that the tunneling in the excited state vibrational wave function becomes important. It may also be that thermally promoted tunneling on the ground PES becomes significant in the deep tunneling region, even though the thermal excitation probabilities are small. However, we leave it for future work to investigate these possibilities.

## IV Dynamics on the excited PES

We now discuss briefly about nuclear dynamics on our excited PES. An interesting measurement of a dynamic time scale on the excited PES was reported in ref. , made by a two-photon pump-probe technique. In the experiment the two photon energies were 6.2 eV and 4.64 eV, and the average time delay between them could be varied. The ionization rate was measured as a function of the delay time. The observed ionization yield peaked when the pulse associated with the lower energy photon was delayed by about 50 fs from the pulse of the more energetic photon. Clearly, the ethylene must first be excited by the 6.2 eV photon before it can absorb the lower energy photon. The experiment was interpreted in terms of lifetime of the state excited by the 6.2 eV photon, with an extracted value of $`30\pm 15`$ fs. The ionization threshold of ethylene is 10.5 eV, so the total energy available is only 0.3 eV above the threshold.

While the complete dynamics on the excited state PES are beyond the TDDFT, the PES we have already constructed can be used to estimate the accelerations in the short-time domain after excitation. According to the calculation in the previous section, the absorption strength at 6.2 eV comes from the torsional zero-point fluctuation the hydrogens rotated to a angle of about $`35^{}`$. Let us assume that the ionization potential of the molecule is independent of the angle. Then, as long as the excitation energy is kept in the electronic degrees of freedom, the ionization can take place. However, the gradient in the PES will cause an acceleration in that coordinate, converting energy from potential to kinetic. The kinetic energy associated with nuclear motion cannot be used in the ionization process, so these accelerations give the first quenching of the ionization rate. To make a concrete estimate, we use the slope of the TDDFT PES to determine the acceleration. At $`35^{}`$, the PES is still rapidly falling, with a slope given by $`dE_{ex}/d\theta 2`$ eV/rad2. Taking this as the force in Newton’s equation, the system acquires a angular velocity given by

$$\frac{d\theta }{dt}=\frac{1}{I_\theta }\frac{dE_{ex}}{d\theta }t$$

with a corresponding kinetic energy

$$K=\frac{1}{2}I_\theta \left(\frac{d\theta }{dt}\right)^2.$$

Evaluating these expression numerically, we find that it only takes 5 fs for the kinetic energy to grow to 0.5 eV. This would put the potential energy below the ionization threshold, and indicates that the quenching should happen very rapidly.

However, the fact that vibronic structure is seen in the strength function tells one that the motion along the torsional coordinate is not very dissipative. The system can accelerate to very large angle, bounce, and return to a state with no kinetic energy in the nuclear motion. Since we cannot say what the true irreversible rates are, the considerations here can only give a lower bound on the effects of a time delay. We note that the calculation of ref.  treats the large angle motion as well and obtains results consistent with the observed time correlation.\[A

## V Conclusions

In this article, we aimed to elucidate how well the TDDFT works in describing electron-nucleus coupling. For this purpose, we analyzed the $`\pi `$-$`\pi ^{}`$ transition of ethylene molecule in the 6-8 eV region. The DFT and TDDFT were employed to construct the ground and excited PES in the LDA. A simplified treatment in the reflection approximation is made for the Frank-Condon factor. The TDDFT describes reasonably well the vertical excitation energy and the magnitude of oscillator strength. The roles of the C-C stretch motion and the torsional motion on the spectral line-shape were investigated. We found that the C-C stretch coordinate contributes substantially to the width of the transition, but the torsional coordinate is responsible for the low energy tail. There the transition strength is reproduced by about a factor of three over 4 orders of magnitude. We thus conclude that the TDDFT survives the first tests as a useful approximation to treat electron-nucleus couplings in excited states. As next steps, a more sophisticated treatment of the quantum nuclear motion would be required, as well as an improved treatment of the exchange-correlation potential beyond the LDA.

## Acknowledgment

We thank T.J. Martinez for helpful discussions, and P.A. Snyder for correspondance leading to a corrected Table 2. This work was supported by the Department of Energy under Grant DE-FG03-00-ER41132, and facilitated by the Computational Materials Science Network.


format: text

text:
# Constraints on Star Formation from the Close Packing of Protostars in Clusters

## 1 Introduction

The mm-wave continuum sources (hereafter MCS) in Ophiuchus (Motte, André & Neri 1998; Johnstone et al. 2000), Serpens (Testi & Sargent 1998; Testi et al. 2000), and Orion (Johnstone et al. 2001) should reveal the process of their own formation. They appear to be the last phase in the gas before it collapses into stars. According to Motte et al., the typical density of an MCS in the Oph B2 core is $`6\times 10^7`$ H2 cm-3, the size is $`3000`$ AU, the mass is $`0.3`$ M, and the projected separations are $`3`$ arcmin $`0.1`$ pc $`20,000`$ AU. This separation corresponds to a spatial density on the order of $`10^3`$ pc-3, considering that the overall core size is about the same 0.1 pc. This density is comparable to that in embedded star clusters (McCaughrean & Stauffer 1994; Carpenter et al. 1997).

The relative velocities between the MCS are probably $`2^{1/2}`$ times the three-dimensional virial speeds in their cloud cores, which for Oph B2 is $`\left(6GM/5R\right)^{1/2}2`$ km s-1 for a core mass of $`M=50`$ M and a core radius $`R=0.055`$ pc (Motte et al.). Sources outside of dense cores may have smaller velocity dispersions (Testi et al. 2000; Belloche, André, & Motte 2001). If the dispersion is small and the MCS are far apart, then they will not interact before they collapse into stars and the final stellar mass function will be about the same as the MCS mass function. However, if the dispersion and density of MCS in a cloud core are high enough, then some may collide before they collapse, and this can affect the mass function (e.g., Bonnell, Bate, & Zinnecker 1998).

The relative collision rate should have been higher in the past if the MCS passed through a lower density state. The collision rate scales approximately with the square of the MCS radius, $`R^2`$, and the collapse rate scales with $`R^{3/2}`$, from the square root of density with mass conservation. The ratio of the collision rate to the collapse rate therefore scales with $`R^{7/2}`$, which is a sensitive function of MCS size. This implies that even if the MCS are not collision-dominated now, they did not have to be much larger in the past before their internal dynamical evolution was severely affected by interactions. If such interactions were destructive for the MCS or detrimental to their stellar-like mass function (Motte et al. 1998), then they could not have formed by contraction from the cloud core through a sequence of nearly constant mass, but had to form from smaller masses by coagulation or accumulation through a sequence of near-constant density, or they had to form elsewhere.

Here we determine the relative collision rates for the MCS in core regions B and C of Ophiuchus. We also determine the relative rates for idealized MCS with radial density profiles from Whitworth & Ward-Thompson (2001) and a mass distribution from the initial stellar mass function (IMF).

The ratio of the collision rate to the collapse rate increases with the mass of an MCS. This means that massive clusters with the same or higher densities will have more severe crowding problems than the cores in Ophiuchus, which is making only low-mass stars. This implies that the contrast between the gas density in an MCS and the gas density in the ambient cloud core must increase with cluster mass or density. This effect is modelled in Section 3.

An important simplification in this problem is that collisional destruction of MCS should be about the same for interactions with other MCS as it is for interactions with protostars and stars. This is because the grazing collision distance between one MCS of mass $`M_i`$ and another of mass $`M_j`$ having the same density is equal to the tidal distance for the destruction of the first MCS by a point source of mass $`M_j`$. Tidal destruction occurs when an MCS enters a region where the average density is the same as its own density or higher. Thus our result is unchanged if only a fraction of the final stars are present as MCS at any one time. Most likely, there is a continuous conversion of existing MCS into stars while new MCS are forming out of the ambient cloud core gas. The ratio of the numbers of MCS to final stars will be about the same as the ratio of the dynamical times, which is the inverse square root of the density contrast.

We also consider whether the turnover mass in the IMF might depend on the cluster core density, taking higher values for denser clusters because of limitations from crowding. The effect is probably small, but it may be observable under extreme conditions. For example, coalescence in the protostellar phase of very dense clusters could shift a significant fraction of low-mass MCS to higher mass. Considering how close the Salpeter IMF is to an $`M^2dM`$ form (this form has equal total mass in equal logarithmic intervals of mass), this shift means that an IMF which extends as a Salpeter power law to low mass but falls below the Salpeter power law at high mass in a low density environment might shift to one that turns over at low mass and has a Salpeter power-law slope at high mass in a high density environment. Observations of the IMF support these changes with stellar density, although other explanations are possible.

## 2 The Relative Collision Rate

The collision rate of an MCS is determined from the space density of its neighbors, the relative velocity dispersion, and the gravitationally-enhanced collision cross section. The collapse rate is determined from the MCS density. The ratio of these two rates indicates the relative probability of collisions versus isolated collapse. If the ratio exceeds unity then collisions should be important. Here we use the observed sizes and masses for MCS tabulated by Motte, Andre & Neri (1998) and we determine the ratio of these rates to see how important collisions are. We then consider the same ratios again for the MCS with hypothetically lower densities. The results are then generalized to the case in which the number of MCS is given by the initial mass function and each MCS has an internal density profile from Whitworth & Ward-Thompson (2001).

### 2.1 Observed mm-wave continuum sources in Ophiuchus

A field of $`N`$ MCS with radii $`R_j`$, mass $`M_j`$, and uniform velocity dispersion $`v`$ inside a cloud core with radius $`R_{cloud}`$ and mass $`M_{cloud}`$ has grazing (hard-sphere) collisions with a particular MCS of radius $`R_i`$, mass $`M_i`$, and velocity $`v`$ at the rate

$$\omega _i=\frac{\pi }{4\pi R_{cloud}^3/3}\underset{j=1}{\overset{j=N}{}}\left(R_i+R_j\right)^22^{1/2}v\left[1+\frac{2G\left(M_i+M_j\right)}{2v^2\left(R_i+R_j\right)}\right].$$
(1)
The virial speed in the cloud core gives the three-dimensional velocity dispersion,

$$v=\left(3GM/5R\right)^{1/2},$$
(2)
which is an approximation valid for a uniform cloud; the extra factors of $`2^{1/2}`$ near $`v`$ and $`2`$ near $`v^2`$ in equation (1) are from relative motion assuming the same speed for all MCS. The collapse rate for a uniform MCS of density $`\rho _i=3M_i/\left(4\pi R_i^3\right)`$ is

$$\omega _{\mathrm{collapse},\mathrm{i}}=\left(32G\rho _i/3\pi \right)^{1/2}.$$
(3)
The ratios of the collision rates to the collapse rates for MCS in the B and C cores of Ophiuchus are shown in Figure 1, plotted as solid symbols versus the MCS mass. The ratios are in the range from 0.001 to 0.1, increasing slightly with mass for Oph B. Such low values indicate that collisions are slower than internal dynamical processes. If the MCS are gravitationally unstable to collapse into stars, then they will do this independently, without interactions.

The same ratios are plotted again in Figure 1 as open symbols for a hypothetical previous state in which each MCS was larger than it is today by a factor of 2.5. The MCS masses are assumed to be unchanged. A small difference in MCS size corresponds to a large difference in relative collision rate, so with larger sizes, some MCS have ratios of rates larger than 1. These MCS could not have reached their present state by purely dynamical collapse without strongly interacting with previous neighbors. Internal evolution that is slower than the collapse rate makes previous interactions even more likely, as would the presence of a larger number of MCS in the previous cloud core.

The dynamical time scale in the present-day MCS is only $`4000`$ years at the average MCS density of $`6\times 10^7`$ H2 cm-3. The dynamical time in each whole core is $`10`$ times longer because of the $`100\times `$ lower average core density. If the MCS and the cores around them both evolved dynamically on their own time scales, then the MCS would have had twice their present radii only several $`\times 10^4`$ years ago, when the cores as a whole were not much different than they are today. This is such a short time ago compared to the age of the whole Ophiuchus star forming region that it seems unlikely the MCS recently suffered some collision without producing tidal features such as bridges or tails. Most likely the MCS have not been in a lower density state since they were part of the ambient cloud core gas.

### 2.2 Protostars with Smooth Density Profiles

The MCS in Ophiuchus are denser than the average cloud cores by a factor of $`100`$, suggesting a mismatch in thermal pressure. The boundaries of the MCS must therefore have either a gradual density gradient, taking the MCS density down to the cloud core density, or a shock front at which the core thermal pressure balances an inter-MCS ram pressure from relative motions or turbulent flows. Boundaries like this are required whether or not the MCS are self-gravitating.

Dense mm-wave sources are often observed to have density gradients like the one fitted by Whitworth & Ward-Thompson (2001), which is

$$\rho (r)=\frac{\rho _{\mathrm{flat}}}{\left[1+\left(r/R_{\mathrm{flat}}\right)^2\right]^2},$$
(4)
for radius $`R_{\mathrm{flat}}`$ and density $`\rho _{\mathrm{flat}}`$ in the central region where the profile flattens. The density contrast is represented by

$$𝒞=\rho _{\mathrm{flat}}/\rho _\mathrm{s},$$
(5)
where $`\rho _\mathrm{s}`$ is the surface density of each MCS. The parameter $`𝒞`$ depends on the definition of the surface. We assume that the MCS surface occurs where its density equals the average density in the cloud core (e.g., $`𝒞=50`$ is the ratio of the average MCS density to the average cloud core density in Oph B2; $`𝒞=187`$ is the average for Oph C – Motte et al. 1998).

Equations (4) and (5) give the radius of the MCS, $`R_s`$, as a function of $`𝒞`$ and $`R_{\mathrm{flat}}`$,

$$R_\mathrm{s}=aR_{\mathrm{flat}}\mathrm{for}a=\left(𝒞^{1/2}1\right)^{1/2}.$$
(6)
For mass

$$M=_0^{R_\mathrm{s}}4\pi r^2\rho (r)𝑑r,$$
(7)
$$R_{\mathrm{flat}}=\left(\frac{M}{2\pi \rho _{\mathrm{flat}}}\right)^{1/3}\left(\mathrm{arctan}(a)\frac{a}{1+a^2}\right)^{1/3}.$$
(8)
The average density of an MCS is

$$\rho _{\mathrm{av}}=\frac{3M}{4\pi R_\mathrm{s}^3}=\frac{3\rho _{\mathrm{flat}}}{2a^3}\left(\mathrm{arctan}(a)\frac{a}{1+a^2}\right),$$
(9)
giving an average collapse rate

$$\omega _{\mathrm{collapse},\mathrm{a}}=\left(32G\rho _{\mathrm{av}}/3\pi \right)^{1/2}.$$
(10)
The collapse rate in the flat part of the density profile is

$$\omega _{\mathrm{collapse},\mathrm{f}}=\left(32G\rho _{\mathrm{flat}}/3\pi \right)^{1/2}.$$
(11)
A cluster of MCS with smooth density profiles has a higher relative collision rate than what we have just calculated for the MCS in Ophiuchus because smoothly tapered MCS are bigger and their average densities are comparable to the cloud core density, giving them a longer collapse time. For a general case, the collision rate for an MCS depends on the space density of other MCS, which is a function of mass through the IMF.

A convenient form for the IMF at low to moderate mass is

$$n(M)dM=n_0M^{2.3}\left[1\mathrm{exp}\left(M^2/M_\mathrm{t}^2\right)\right]dM,$$
(12)
which has a Salpeter slope of $`2.3`$ toward higher mass and a turn-over toward lower masses at $`M_\mathrm{t}`$ (see Elmegreen 1997). The turnover mass may equal the thermal Jeans mass unless crowding limitations are important; we discuss these limitations in Section 4. We introduce a high-mass modification to the IMF in Section 3 (see Fig. 3).

To get $`n_0`$, we set the lower mass limit from observations of Ophiuchus, we equate the total MCS mass in the IMF to the observed total, and we constrain the number of MCS with a mass greater than some maximum mass to be unity. In this way, the IMF has the correct total mass and it also falls off gradually at high mass to give only one most massive star. The minimum mass, $`M_{\mathrm{min}}`$ is set equal to 0.1 M for both Oph B2 and Oph C (see Fig. 1). The total MCS masses in both Ophiuchus B and C are equal to the fraction $`ϵ=0.13`$ times the cloud core masses.

The collision rate for MCS $`i`$ is now obtained by integrating over the IMF,

$$\omega _{\mathrm{collision}}\left(M_i\right)=\frac{\pi 2^{1/2}v}{4\pi R_{cloud}^3/3}_{M_{min}}^{\mathrm{}}n(M_j)(R_i+R_j)^2[1+\frac{2G\left(M_i+M_j\right)}{2v^2(R_i+R_j)]}.]dM_j$$
(13)
The results of this model are shown in Figure 2. The ratios of the collision rates to the average collapse rates for MCS are plotted as a function of MCS mass. Solid lines are for an idealized model of Oph B and dashed lines are for Oph C. Two different values of the turn-over mass in the IMF, $`M_t`$, are assumed. The lower set of curves is for the MCS cores only, which come from the flat part of the Whitworth & Ward-Thompson profile only. That is, it uses $`\rho _{\mathrm{flat}}`$ for the collapse rate, and $`R_{\mathrm{flat}}(M)`$ for the collision cross section.

Figure 2 reproduces the ratios of the collision rates to the collapse rates found in Figure 1 when only the flat parts of the density profiles are considered. These ratios are small, so the flat parts should not collide with each other before they collapse into stars. However, if the MCS have envelopes where the density drops down to the average value in the cloud core, which corresponds to an extra factor of $`a3`$ in radius, then the envelopes will interact significantly before they get involved with the collapse. Thus MCS in dense cluster environments should not have significant envelopes.

## 3 Globular Clusters

Other models were considered for more massive clusters. When the cluster mass is large, the largest MCS mass is also large, forcing us to modify the IMF so that, for example, 1000 M stars do not form. We follow the discussion in Elmegreen (2000) and use for this case an IMF of the form:

$$n(M)dM=n_0M^{2.3}\left[1\mathrm{exp}\left(M^2/M_\mathrm{t}^2\right)\right]\mathrm{exp}\left(M/100\mathrm{M}_{}\right)dM,$$
(14)
which has a cutoff at high mass. The IMFs in equations (12) and (14) are shown in Figure 3 (multiplied by $`M`$ to convert to a $`\mathrm{log}M`$ coordinate). The dotted line is from equation (12) for $`M_t=0.1`$; it blends with the solid line at low mass, which is from equation (14) using the same $`M_t`$. The dashed line is from equation (14) using $`M_t=1`$ M and the same total cluster mass as the solid line. The two straight lines show slopes of $`1.3`$ and $`1.6`$ for comparison; these bracket the IMF observations in the indicated mass range. There are no observations of an IMF for masses much larger than 100 M, so the fall off is purely conjectural at this time.

The result for massive clusters are shown in Figures 4 and 5. Figure 4 shows the ratio of rates versus MCS mass for two clouds with $`M_{\mathrm{cloud}}=6000`$ and $`6\times 10^5`$ M, and with equal star formation efficiencies ($`ϵ=0.3`$) and densities ($`\rho _{\mathrm{cloud}}=7.8\times 10^{20}`$ gm cm-3), corresponding to radii of $`R_{\mathrm{cloud}}=`$ 1.08 pc and 5 pc, respectively. The density contrast, $`𝒞=160`$, was also taken to be the same in each case, and comparable to that in Ophiuchus. Two values of the turnover mass, $`M_t`$ are considered, while $`M_{\mathrm{min}}=0.1`$ M for all cases.

Figure 4 shows several things. First, the maximum mass increases with the total stellar mass, $`ϵM_{\mathrm{cloud}}`$, as expected for a statistical ensemble (i.e., the solid lines go to larger protostellar mass than the dashed lines). Second, the ratio of collision to collapse rates increases with MCS mass because of the increasing collision cross section. Third, the ratio is much lower for the MCS cores than for envelopes, as in Figure 2; it is also lower for higher turnover mass, $`M_t`$. The reason for the $`M_t`$ dependence is that higher $`M_t`$ corresponds to relatively fewer low mass stars, and this decreases the number density for collisions in the field.

Figure 4 suggests a problem with star formation in massive clusters. If the densities of the MCS-equivalent objects are the same as in the Ophiuchus cores, then the envelopes of the massive stars will strongly interact with most of the lower mass MCS, probably accreting them. Even the cores of the MCS in the highest mass clouds (solid lines in the figure) will interact (by tidal forces) strongly with other MCS and with the protostars already formed. This implies that massive clusters cannot have MCS with the same densities as those in Ophiuchus (i.e., $`10^7`$$`10^8`$ m(H2) cm-3).

Figure 5 shows the same clouds and clusters as in Figure 4 but now with compression factors ($`𝒞=59`$, 3460) that scale with the first (dashed lines) and second (solid lines) powers of the Mach number from virialized motions, assuming a sound speed of 0.3 km s-1. The assumption here is that MCS might be produced in self-gravitating shocks that have average compression factors proportional to some power of Mach number. This power is likely to be close to 1 when the shock moves across magnetic field lines, and 2 when the shock moves parallel to the field lines. The shocks presumably arise from supersonic turbulence in the cloud cores where MCS form, so the turbulent speed is taken to be the cloud virial speed, which is the same as $`v`$ in equation (13).

Figure 5 also shows cases with very large turnover masses, $`M_t=10`$ M. When $`M_t1`$ M, only the Mach-squared scaling is strong enough to make the MCS so small that their cores avoid severe collisions. In all cases, the envelopes collide. McKee & Tan (2002) discuss the last stages of massive star formation and note that it requires very high initial densities and pressures. This is consistent with the high value of $`𝒞`$ found here from close-packing constraints in the massive clusters where massive stars tend to form. One implication of this result is that the formation of stars from MCS in massive clusters should be faster than in low mass clusters.

## 4 IMF Dependence on Cluster Mass and Density

Figures 4 and 5 suggest that the MCS mass where collisions are important ($`\omega _{\mathrm{collision}}/\omega _{\mathrm{collapse}}1`$) becomes insensitive to the IMF turnover mass, $`M_t`$, when the ratio between the maximum stellar mass and $`M_t`$ is large. Once this ratio exceeds several hundred, most collisions are between MCS with masses above $`M_t`$ and then the value of $`M_t`$ does not matter.

Higher $`M_t`$ means fewer low-mass collision partners and more freedom for the massive MCS to collapse without interference from neighbors. For a given modest MCS density contrast, $`𝒞`$, denser and more massive stellar clusters should have higher IMF turnover masses because of crowding limitations. However, sufficiently large $`𝒞`$ gives all cores a low ratio of collision to collapse rates, and then the IMF does not depend on crowding.

An increase in thermal temperature in a cloud core should decrease the density contrast in shocks and decrease $`𝒞`$. This increases the relative importance of collisions and can increase $`M_t`$ by removing the low mass MCS. The dependence of $`M_t`$ on $`𝒞`$ for a constant ratio of rates is approximately $`M_t𝒞^2`$ from Figure 5. This comes from the mutual dependencies of $`\omega _{\mathrm{collision}}/\omega _{\mathrm{collapse}}`$ on $`𝒞^{3.5/3}`$ (see Sect. 5) and $`M_t^{0.6}`$ for the “core only” case with $`M_t`$ in the range from 1 to 10. From these mutual dependencies, $`\omega _{\mathrm{collision}}/\omega _{\mathrm{collapse}}`$ is constant if $`M_t^{0.6}𝒞^{3.5/3}`$, which gives $`M_t𝒞^{1.94}`$.

Such a $`𝒞`$ dependence is comparable to the dependence of the thermal Jeans mass ($`M_J`$) on compression ratio, which would predict $`M_J𝒞^2`$ for $`𝒞`$ proportional to the square of the Mach number. This comes from the dependence of the Jeans mass on $`T^2`$ for constant pressure, and the dependence of the Mach number on $`T^{1/2}`$ for constant virial speed in the cloud.

The similar dependencies of $`M_t`$ and $`M_J`$ on $`𝒞`$ imply that the thermal Jeans mass and the crowding limitations both determine the turnover mass in the IMF. If two regions have different temperatures with all else the same, then $`M_t`$ should be larger in the warmer region in proportion to the thermal Jeans mass and the crowding conditions will stay about the same. This means that warm cloud cores can produce a top-heavy IMF because of a combination of increased thermal Jeans mass and increased crowding.

No clear dependencies between cluster density and IMF have been found (Massey & Hunter 1998), although some IMF observations suggest they could be present. For example, the clustered embedded stars in Orion have a shallower IMF than the unclustered embedded stars (Ali & DePoy 1995), and the LMC clusters in regions of high stellar density tend to have shallower slopes than those in regions of low stellar density (J.K. Hill et al. 1994; R.S. Hill et al. 1995). Similarly, the extreme field stars in the LMC have a steeper IMF slope (Massey et al. 1995) than the dense and massive cluster 30 Dor, which has a shallow slope like the Salpeter value (Massey & Hunter 1998). Shallow slopes correspond to proportionally more massive stars, and this can be the result of a previous stage in the evolution of the cluster where the low mass MCS coalesced into higher mass MCS. The problem with these interpretations is that other processes, such as mass segregation, can change the slope of the IMF in the same way as MCS coalescence.

Variations in $`M_t`$ have also been difficult to establish observationally. The old globular clusters in the Milky Way halo tend to have normal $`M_t`$, namely about half a solar mass (Paresce & De Marchi, 2000), whereas a young globular cluster in M82 has $`M_t`$ higher than normal by a factor of $`5`$ (Smith & Gallagher 2001). The cluster masses and densities are about the same.

This lack of a clear correlation between the IMF and the density or mass of a cluster (over and above possible changes in the thermal Jean mass) implies that crowding effects are not particularly important in either the protostellar phase or in an earlier phase, such as that producing the MCS in Ophiuchus. This can only be the case if the MCS that eventually form stars build up from smaller objects at nearly constant density, and if this density scales with at least the square of the Mach number in the cloud core.

Crowding and coalescence in the MCS phase will be very important in the cores of globular clusters if the temperature is large or the compression factor scales directly with the Mach number, as for a strongly magnetic gas. In this case, supermassive stars can form by the coalescence of MCS. Such stars may either form massive black holes directly, or they may continue to coalesce as stellar remnants until several together make a massive black hole. This is the scenario for black hole formation in dense clusters proposed by Ebisuzaki et al. (2001), but modified slightly here to include the MCS phase, when collisions are more important than in the main-sequence or post-main sequence phases.

## 5 Scaling of the Results

The curves in Figures 4 and 5 have a similar dependence on MCS mass that comes mostly from the mass dependence for the square of the MCS radius, which is approximately as $`M^{2/3}`$. Lower $`M_t`$ clusters have slightly higher ratios of rates because of the correspondingly larger MCS spatial densities. Lower $`M_t`$ clusters also produce slightly steeper slopes in the figures. The steeper slope at low $`M_t`$ arises because the terms $`(R_i+R_j)`$ and $`(M_i+M_j)`$ are dominated by the radius and mass of the test object, $`R_i`$ and $`M_i`$, when most of the field MCS are at much smaller radii and masses. This occurs at low $`M_t`$, considering that the integral in equation (13) is dominated by the field MCS mass at the peak of the IMF, which is close to $`M_t`$. When $`M_t`$ is high, most of the field interactors are large and massive, and then the dependencies of $`(R_i+R_j)`$ and $`(M_i+M_j)`$ on $`R_i`$ and $`M_i`$ alone diminish. The test mass, $`M_i`$, is the quantity plotted on the abscissa in the figures.

The vertical shifts between curves may be understood in terms of the dependencies of the ratio of rates on the compression factor, $`𝒞`$, and the cloud mass, $`M_{cloud}`$. In Figure 4, the cloud mass in two cases differs by a factor of 100 while the cloud density and the MCS densities and radii are the same. The cloud mass enters into the velocity dispersion $`v`$, the cloud radius, $`R_{cloud}`$, and the total number of MCS, which appear in equation (13). For the assumed constant cloud density, the only remaining dependence is on $`v`$, which scales as $`M_{cloud}^{1/3}`$. Thus the dashed and solid lines for both core and envelope models are displaced vertically in the figure by $`100^{1/3}=4.64`$, which is 0.66 in the log, as observed.

In Figure 5, the compression factor in two cases differs by a factor of 58.6 for the same cloud mass and radius. For the envelopes, the MCS radii scale with $`𝒞^{1/4}`$ times the core radii, and the core radii scale with $`𝒞^{1/3}`$ for a given MCS mass. Thus the MCS envelope radii scale with $`𝒞^{1/12}`$ at constant mass when $`𝒞`$ is large. The ratio of collision to collapse rates scales with radius as $`R^{3.5}`$, which is as $`𝒞^{3.5/12}`$. Considering the ratio 58.6 for $`𝒞`$, this means that the dashed and solid curves for the envelope case in Figure 5 should be vertically displaced by $`58.6^{3.5/12}=0.30`$, which is $`0.5`$ in the log, as observed approximately. For the cores, the MCS radii scale with $`𝒞^{1/3}`$ for a given MCS mass and so the ratio of rates scales with $`𝒞^{3.5/3}`$. This gives a vertical displacement between the dashed and solid lines for the cores that is equal to a factor of $`58.6^{3.5/3}=0.0085`$, which is $`2.1`$ in the log, as observed.

The displacement between the envelope and core cases results mostly from the $`𝒞^{1/4}`$ ratio between envelope and core radii, given the $`R^{3.5}`$ dependence for the ratio of rates at constant MCS mass. Thus in Figure 5, the vertical displacement for the $`𝒞=59`$ case is $`59^{3.5/4}=35.4`$, which is $`1.55`$ in the log, and the vertical displacement for the $`𝒞=3460`$ case is $`3460^{3.5/4}=1249`$, which is $`3.1`$ in the log. Both are approximately correct compared to the dislacements in the Figure.

## 6 Conclusions

The mm-wave continuum sources in Ophiuchus are not interacting strongly at the present time and should collapse to stars independently, preserving their masses (Motte et al. 1998). Previous versions of these MCS should have interacted more strongly if their densities were lower at fixed masses, and MCS in more massive clusters should interact more strongly too. If such interactions are destructive, or if they modify the mass function, then such a previous stage could not have occurred. This might imply that the MCS formed elsewhere and fell into the cloud core after they had their high densities. However, such a model would not explain the hierarchical positions of these MCS, as seen in both Ophiuchus and Serpens (Testi, et al. 2000); MCS ballistic motions over one or more core dynamical times would mix them up.

A high ratio of rates might also imply that collisions were important in the past, particularly among the most massive MCS. However, pure collisional agglomeration leads to a shallow mass function, $`M^{1.5}dM`$ (Field & Saslaw 1965), not the observed $`M^{2.3}dM`$ function at high MCS mass (Motte et al. 1998; Testi & Sargent 1998; Johnstone et al. 2001). Accretion of ambient gas by moving protostars may give a better mass function (Bonnell et al. 1998), but such models are not collision-dominated like the high $`\omega _{\mathrm{collision}}/\omega _{\mathrm{collapse}}`$ cases here.

The most likely scenario for the formation of the MCS in Ophiuchus and elsewhere is that they grow inside high-density turbulent shocks and become unstable when they reach a sufficiently large mass. Such a process would not have a previous stage where low density MCS were strongly interacting. It would also account for the hierarchical positioning of MCS in cloud cores (because turbulence makes hierarchical structure), the inferred Mach number dependence of the MCS density with increasing cloud core mass, the lack of smoothly varying envelopes, and the required rapid formation time of the MCS, considering their high densities. Models of turbulence-induced star formation are in Elmegreen (1993), Padoan (1995), Heitsch, Mac Low, & Klessen (2001), Klessen, Heitsch, & Mac Low (2000), Klessen (2001), Ossenkopf, Klessen, & Heitsch (2001), Padoan, et al. (2001), and Wuchterl & Klessen (2001). Other formation scenarios might be preferred outside of dense clusters.

For turbulence-induced star formation, the number of stars that finally forms in a cluster should be larger than the number of MCS that are present at any one time because the MCS phase is shorter than the lifetime of the whole core. The ratio of the number of MCS to the number of final star counts should be the ratio of corresponding evolution times, which is the inverse square root of the ratio of densities, $`𝒞^{1/2}`$. This ratio equals $`0.1`$ for the Ophiuchus cores and could be as small as $`0.01`$ for super star clusters ($`M10^5`$ M). The MCS-equivalent phase in massive clusters should also be associated with higher densities, possibly in proportional to the square of the Mach number. Thus objects with densities of $`10^7`$ cm-3 in Ophiuchus would appear in a similar phase of evolution having densities of $`10^910^{10}`$ cm-3 in super star clusters. Of course, observations with a sensitivity to a particular density will always see the observable value, so comparisons between low and high mass clusters should look for the same pre-stellar phases in terms of structure and kinematics when comparing densities.

The IMF should be affected by crowding in the sense that warmer, more crowded environments should have a higher low-mass limit to the power law part and proportionally more massive stars to conserve total mass. The temperature dependencies of the thermal Jeans mass and the low mass limit from crowding are about the same.