text:
# Evaporation rate of hot Jupiters and formation of Chthonian planets

## Abstract

Among the hundred of known extrasolar planets, about 15% are closer than 0.1 AU from their parent stars. But there are extremely few detections of planets orbiting in less than 3 days. At this limit the planet HD 209458b has been found to have an extended upper atmosphere of escaping hydrogen. This suggests that the so-called hot Jupiters which are close to their parent stars could evaporate.

Here we estimate the evaporation rate of hydrogen from extrasolar planets in the star vicinity. With high exospheric temperatures, and owing to the tidal forces, planets evaporate through a geometrical blow-off. This may explain the absence of Jupiter mass planets below a critical distance from the stars. Below this critical distance, we infer the existence of a new class of planets made of the residual central core of former hot Jupiters, which we propose to call the “Chthonian” planets.

Institut d’Astrophysique de Paris, CNRS, 98bis boulevard Arago, F-75014 Paris, France

Following the recent discovery of a significant escape of atomic hydrogen from the planet orbiting HD 209458 (Vidal-Madjar et al. 2003), we propose to evaluate the escape flux from the upper atmosphere of hot Jupiters under the influence of the strong tidal forces from their parent stars. This escape flux had previously been estimated (Guillot et al. 1996), but the conclusion that the mass loss is not significant was based on two hypotheses which need to be revisited.

First, the black body radiative equilibrium used to calculate the temperature of the upper atmosphere (Schneider et al. 1998) is inappropriate because it does not apply to the low density upper atmosphere. Observations in the Solar System show that the temperature of a planetary upper atmosphere (thermosphere, exosphere) is much higher than the effective temperature of the bottom atmosphere. For example, whereas the temperatures in the Earth and Jupiter are around 200 K and 120 K respectively at the level of the tropopause, they reach 1000 K in the thermosphere of these two planets. Although these observed high temperatures in the planets of our own Solar System remain unexplained, there are some clues that a combination of the extreme and far ultraviolet fluxes with the Solar wind is responsible for the heating.

The second important hypothesis is tidal force, which modifies the gravity of hot Jupiters. The common hypothesis is to neglect that effect as for isolated planets far from their star. However, tidal forces have a significant influence on the density distribution in the upper atmosphere of hot Jupiters.

At a given temperature the escape flux can be calculated using the Jeans’ escape estimate. Jeans’ escape refers to the escape of atoms and molecules whose velocities are in the tail of the Boltzmann distribution and which have enough energy to escape the planet gravity. The flux is calculated at the top of the thermosphere, that is the level above which the atoms and molecules have no collisions. The exobase is usually defined as the location above which the mean free path is larger than the scale height of the atmosphere. Here we extend this idea by defining the exobase as the location above which the mean free path is larger than the distance to the Roche lobe of the planet. Thus, the exobase is the location above which the atoms and molecules can definitively escape the planet.

We computed the escape flux of atomic and molecular hydrogen as a function of the upper atmosphere temperature. From this, we can derive the corresponding life time needed to evaporate the total mass of the planet. For $`T\begin{array}{c}>\\ \end{array}7000`$ K, the H$`\mathrm{i}`$ escape flux is larger than the minimum flux of $`10^{10}`$ g s-1 determined from Ly-$`\alpha `$ observations (Vidal-Madjar et al. 2003).

Although the mechanism responsible for the heating of the upper atmosphere of the planets in the Solar System is not fully identified, we can hint a simple estimate of a plausible temperature from the comparison of the heating and cooling mechanisms. A lower limit of the heating can be estimated from the energy flux of the stellar extreme ultraviolet and the Ly-$`\alpha `$ photons. The cooling is due to a combination of the heat conduction toward the cooler bottom atmosphere, the collisional ionization and excitation of the H$`\mathrm{i}`$ electronic levels, and the evaporation at the top of the upper atmosphere. Thus, a lower limit to the temperature can be estimated by the energy balance. With the parameters of HD 209458b, we evaluate the temperature in the upper atmosphere to be at least $`10000`$ K. With this temperature the present evaporation rate of hydrogen from HD 209458b must be $`1`$ to $`5\times 10^{11}`$ g s-1. This is in agreement with the lower limit found from the HST Ly-$`\alpha `$ observations. Assuming an age of $`5\times 10^9`$ years, HD 209458b may have lost about 1% to 7% of its mass.

Finally, we also estimated the life time of a given planet as a function of its mass and orbital distance. We found that planets with orbital periods shorter than 2 to 3 days have very short life time. This may explain that only few planets have been detected with periods shorter than 3 days. Moreover, low-mass planets have also short life time. This means that the nature of these low-mass hot Jupiters should evolve with time. These planets should loose a large fraction of their hydrogen. This process may lead to planets with an hydrogen-poor atmosphere, or even with no more atmosphere at all. The emergence of the inside core of former and evaporated hot Jupiter may give rise to a new class of planets which we proposed to call “Chthonian” planets in reference to the Greek god Khthôn (“Chthonian” is used to name the Greek deities who come from hot infernal underground).

References:

Guillot, T., Burrows, A., Hubbard, W. B., et al. 1996, ApJL, 459, L35

Schneider, J., Rauer, H., et al. 1998, ASP Conf. Ser., 134, 241

Vidal-Madjar, A., et al. 2003, Nature, 422, 143 (and this proceedings)

We thank J. C. McConnell, C. D. Parkinson, G. E. Ballester & M. Mayor for their help.


format: text

text:
# Bulk antiferromagnetism in 𝐍𝐚_0.82𝐂𝐨𝐎_𝟐 single crystals

## Abstract

Susceptibility, specific heat, and muon spin rotation measurements on high-quality single crystals of $`\mathrm{Na}_{0.82}\mathrm{CoO}_2`$ have revealed bulk antiferromagnetism with Néel temperature $`\mathrm{T}_\mathrm{N}=19.8\pm 0.1`$ K and an ordered moment perpendicular to the $`\mathrm{CoO}_2`$ layers. The magnetic order encompasses nearly 100% of the crystal volume. The susceptibility exhibits a broad peak around 30 K, characteristic of two-dimensional antiferromagnetic fluctuations. The in-plane resistivity is metallic at high temperatures and exhibits a minimum at $`\mathrm{T}_\mathrm{N}`$.

Long-standing interest in the interplay between spin, charge, and orbital degrees of freedom in cobalt oxides has increased following the recent discovery of superconductivity in the hydrated cobaltate $`\mathrm{Na}_{0.35}\mathrm{CoO}_21.3\mathrm{H}_2\mathrm{O}`$ (Ref. takada03, ). Its transition temperature, $`\mathrm{T}_\mathrm{c}=4.6\mathrm{K}`$, is much lower than those of the layered copper oxides, but higher than that of $`\mathrm{Sr}_2\mathrm{RuO}_4`$, the only other known layered transition metal oxide which exhibits superconductivity. The crystal structure of this compound contains distorted, edge-sharing $`\mathrm{CoO}_6`$ octahedra such that the Co sites form a triangular lattice lynn03 ; jorgensen03 . The electronic structure is quasi-two-dimensional singh00 ; koshibae03 ; honerkamp03 ; renner03 . Unconventional superconductivity has been proposed based on NMR experiments waki03 ; fujimoto03  and model calculations singh00 ; baskaran03 ; ogata03 ; motrunich03 , but many important questions about the symmetry of the order parameter and the character of the quasi-particles at the Fermi surface have yet to be answered.

The triangular $`\mathrm{CoO}_2`$ layers of $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ are isostructural to those of $`\mathrm{Na}_{0.35}\mathrm{CoO}_21.3\mathrm{H}_2\mathrm{O}`$ (Refs. fouassier73, ; jansen74, ; balsys96, ). It may therefore serve as a reference system with respect to the effects of doping and magnetic correlations. $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ is metallic tanaka94  and exhibits an unusual Hall effect wang03 , as well as a very large thermopower terasaki97 . The thermopower increases with increasing x motohashi01 , which in a local-moment picture corresponds to a progressive dilution of magnetic $`\mathrm{Co}^{4+}`$ ($`S=1/2`$) with nonmagnetic $`\mathrm{Co}^{3+}`$ ($`S=0`$). Following observations of a pronounced field dependence, a magnetic mechanism has been invoked to explain the large thermopower wang03 . It is therefore important to determine the magnetic ground state and exchange parameters in this material.

We have taken a first step in this direction by synthesizing high-quality single crystals of $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ with x$`0.8`$ and studying their resistivity, specific heat, and magnetic susceptibility as a function of temperature. An antiferromagnetic state with Néel temperature $`\mathrm{T}_\mathrm{N}=19.8\pm 0.1`$ K and an ordered moment perpendicular to the $`\mathrm{CoO}_2`$ layers was observed. Muon spin rotation ($`\mu `$SR) measurements confirm that the static magnetic order encompasses nearly 100% of the volume. The in-plane resistivity shows a minimum around the magnetic phase transition. The susceptibility exhibits a broad peak around 30 K, characteristic of two-dimensional antiferromagnetic fluctuations. These findings may be indicative of a quasi-two-dimensional commensurate spin density wave.

Cylindrical single crystals of diameter 6 mm and length 80 mm were grown in an optical floating-zone furnace. The initial polycrystalline material was prepared using a mixture of $`\mathrm{Na}_2\mathrm{CO}_3`$ and $`\mathrm{Co}_3\mathrm{O}_4`$ with a Na:Co ratio of 0.8:1. The powders were calcined at 750 C for 12 hours and then reacted at 850 C for a day with intermediate grindings. The mixture was pressed to form a cylinder and premelted prior to growth. The molten zone was passed through the feed rod under oxygen flow at a rate of 2 mm/h. Details will be published elsewhere chen . Pieces cut from two of the resulting crystals were analyzed by inductively coupled plasma atomic emission spectroscopy and atomic absorption spectroscopy. The Na:Co ratios were found to be $`0.82\pm 0.02`$ and $`0.83\pm 0.02`$, respectively. The virtually identical compositions attest to good reproducibility of the crystal preparation. The Néel temperatures of both samples were also identical within experimental error. Single-crystal and powder x-ray diffraction confirm a nearly pure phase with space group P63/mmc, in good agreement with Ref. balsys96,  and with Raman scattering experiments raman . The lattice parameters were determined by single-crystal x-ray diffraction along $`(h,h,0)`$ and $`(0,0,l)`$ at room temperature, resulting in values of $`a=2.843(1)\mathrm{\AA }`$ and $`c=10.687(5)\mathrm{\AA }`$. The susceptibility, specific heat, and resistivity measurements were carried out in Quantum Design MPMS and PPMS systems.

Fig. 1 shows the magnetic susceptibility measured in a 5T magnetic field oriented either along the $`c`$-axis or perpendicular to it. Upon cooling, the anisotropy increases, and a broad peak develops around $`31\pm 3`$ K. This is indicative of short-range, Ising-type, quasi-two-dimensional antiferromagnetic correlations. A sharp decrease of the $`c`$-axis susceptibility around 20 K heralds an antiferromagnetic transition with ordered moments along $`c`$. The inset in Fig. 1 shows that the transition is sharp and only weakly affected by the field. Curie tails, presumably due to paramagnetic impurities, are suppressed in a 5T field.

A broad hump in the susceptibility is evident at around $`285\pm 5\mathrm{K}`$, for all sample orientations and at all applied magnetic fields. This temperature corresponds roughly to the onset of the low-temperature splitting of a phonon mode at $`\nu 550\mathrm{c}\mathrm{m}^1`$, as observed in IR conductivity measurements (not shown). This suggests that this feature in the susceptibility arises from a structural phase transition intrinsic to $`\mathrm{Na}_{0.8}\mathrm{CoO}_2`$, rather than from an impurity phase such as CoO (an antiferromagnet with $`\mathrm{T}_\mathrm{N}=292\mathrm{K}`$). This interpretation is consistent with x-ray powder diffraction measurements, which exhibit no traces of CoO. The origin of the structural transition remains to be determined.

The specific heat, $`C_P`$, also exhibits a sharp anomaly at $`19.8\pm 0.1`$ K, which indicates the onset of long-range antiferromagnetic ordering (Fig. 2). The entropy contained in this anomaly amounts to 0.08(1) J/molK and corresponds to about 10% of the entropy of the Co4+ spin-1/2 system. There is a slight excess heat capacity above $`\mathrm{T}_\mathrm{N}`$, extending up to $``$ 35 K. These observations are consistent with the short-range magnetic fluctuations above $`\mathrm{T}_\mathrm{N}`$ visible in the magnetic susceptibility measurements.

The specific heat at low temperatures contains a contribution linear in temperature which is revealed in a $`C_P/T=\gamma +\beta T^2`$ plot (inset in Fig. 2). It is slightly sample-dependent and amounts to 8.4(3) and 10.4(4) mJ/molK2 for the two measured crystals. We ascribe it to the Sommerfeld term from the conduction electrons. If the cubic term $`\beta `$ is attributed to the contribution of acoustic phonons alone (ignoring possible contributions from magnons), it corresponds to a Debye temperature of $``$ 420 K.

Muon-spin-rotation ($`\mu `$SR) experiments have been performed using the GPS setup at the $`\pi `$m3 beamline at the Paul-Scherrer-Institute (PSI) in Villigen, Switzerland, which provides 100% spin-polarized muons. The $`\mu `$SR technique is especially suited for the investigation of magnetic materials with small magnetic moments. In particular, it allows one to study the homogeneity of the magnetic state on a microscopic scale, and also to access its volume fraction. The accessible time scale is 10-6-10-10 s. The muons are implanted into the bulk of the material, where they thermalize rapidly without any noticeable loss in their initial polarization. They stop at well-defined (but in general unknown) interstitial lattice sites. In the present compound, they most likely form a muoxyl bond with an oxygen atom in the CoO2 planes or occupy empty Na sites in between the CoO layers. The whole ensemble of muons is randomly distributed throughout a layer with a thickness of about 200 $`\mu `$m, and thus probes a representative part of the sample volume. The muons decay (with a mean lifetime of $`\tau _\mu =`$ 2.2 $`\mu `$s) into two neutrinos and a positron which is preferentially emitted along the direction of the muon spin at the instant of decay. The time evolution of the spin polarization, $`P(t)`$, of the muon ensemble can thus be obtained via the time-resolved detection of the spatial asymmetry of the positron emission rate.

Fig. 3a shows $`\mu `$SR spectra obtained at T=7.5 K in a zero-field configuration. The normalization of the total asymmetry of the polarization function is based on the measurement in a weak transverse field of 100 Oe at 50 K, where the sample is in the paramagnetic state. Spectra for two different orientations of the muon spin-polarization with respect to the crystallographic $`c`$-axis are shown. It is evident that the spectra contain several oscillating components, which means that the muons experience several well-defined local magnetic fields. Even without any further analysis, one can immediately conclude that a large fraction of the sample volume must exhibit static magnetic order. We analyzed the spectra using a relaxation function of the form

$$P(t)=P(0)\underset{𝑖}{}A_i\mathrm{cos}(2\pi \nu _\mu ^it)\mathrm{exp}(\lambda ^it).$$
(1)
An excellent description of the data can be obtained with a sum of three oscillatory terms with frequencies and relaxation rates (at low temperature) $`\nu _\mu =1.25`$, 2.6 and 3.2 MHz and $`\lambda =0.8,`$ 1.2 and 2.8 $`\mu `$s$`^1,`$ plus one non-oscillatory term with a very small relaxation rate of $`\lambda =0.02`$ $`\mu `$s-1 that accounts for the component of the muon spin parallel to the internal fields (lines in Fig. 3a). The frequencies of the oscillating components do not depend on the direction of the muon spin, while their amplitudes, $`A_i`$, exhibit a strong variation. The amplitudes for the different muon-spin orientations can be used to determine the direction of the local magnetic fields, as well as the corresponding magnetic volume fractions. The spectra suggest that the muons experience three different well-defined local magnetic fields which point in different directions. The highest field is parallel to the $`c`$-axis (within the error bar) and amounts to about 17% of the entire signal. The 2.6 MHz signal is at an angle of about 25 with respect to the $`c`$-axis and comprises nearly 33% of the total amplitude. Finally, the 1.2 MHz signal is at an angle of 43 and represents about 50% of the entire signal. Based on this analysis, we thus find that the volume fraction of the magnetic state is in excess of 95%: that is, the entire sample is magnetically ordered within the error of our measurement.

Evidence for antiferromagnetic ordering in NaxCoO2 was also obtained in a previous $`\mu `$SR study sugiyama03 . While the Néel temperatures we find are in good agreement with the data of Sugiyama et al. sugiyama03 , there are also substantial differences. First, they obtained much smaller volume fractions of the magnetically ordered phase (about 20% in their polycrystalline Na0.75CoO2 sample and 50% in their Na0.9CoO2 single crystal). This suggests that their samples might have been less homogeneous, possibly due to an inhomogeneous distribution of Na. Another important difference may exist concerning the intercalation of water, which is known to occur rather rapidly under moist conditions, especially for powder samples. Secondly, the magnitude and number of the deduced muon spin precession frequencies are also rather different. In their polycrystalline samples, they also observe three muon frequencies, at about 3.3, 2.6 and 2.1 MHz. However, their relative amplitudes are different from those observed in our crystals, so that, for instance, the largest contribution comes from the highest-frequency mode. In their single-crystal sample, Sugiyama et al. observe only one frequency, though there is some evidence for at least one additional frequency sugiyama03 .

This leaves us with the important question of whether the different local magnetic fields are a consequence of different crystallographic muon stopping sites, or whether they should be interpreted in terms of different magnetic environments. The latter may be caused by commensurate magnetic order with a large unit cell, but could also arise from a macroscopically inhomogeneous magnetic state. As noted by Sugiyama et al., the local fields experienced by the muons in different crystallographic muon sites should be mostly of dipolar origin. Potential muon stopping sites are close to the oxygen ions and near Na(1) or Na(2) vacancies. We performed a calculation of the dipolar fields for the case of Co moments directed along the $`c`$-axis and exhibiting A-type antiferromagnetic order. We find that the magnetic fields near the Na(1) and the Na(2) vacancies are vanishingly small. We therefore suggest that the muons are located primarily near the oxygen ions, forming a muoxyl bond, as in other oxide compounds including the cuprate high-$`\mathrm{T}_\mathrm{c}`$ superconductors muonsite . Assuming an ordered moment of 0.3 $`\mu _B`$, we obtain a local magnetic field of about 2 MHz at the oxygen site.

The three different muon-spin precession frequencies observed in the experiment could be due to the presence of a structural minority phase with a monoclinic distortion, as discussed in Ref. fouassier73, . However, if a minority phase were present in our sample, its Néel temperature would have to be very similar to that of the majority phase. This conclusion is based on the temperature dependence of the local magnetic fields, which is shown in Fig 3b. It is evident that all three local fields exhibit a similar T-dependence, i.e. they decrease towards the magnetic transition at 19.8 K evident in the magnetization and the specific heat data.

It is also possible that the three frequencies reflect a commensurate spin order in which the same local field is repeated in every third crystallographic unit cell. However, our data do not support an interpretation in terms of an incommensurate magnetic order sugiyama03 , which should give rise to much larger depolarization rates and to a different form of the relaxation function. We have tried fitting with Bessel functions, but could not obtain a reasonable fit. Since $`\mu `$SR data provide only local information, it is not straightforward to deduce further information concerning the magnetic ordering pattern. Neutron scattering experiments will be required to settle this important issue.

Fig. 4 shows the in-plane resistivity measured on a rectangular crystal (4 $`\times `$ 3 $`\times `$ 0.6 mm3) with a four-point ac-technique ($`\nu `$ = 119 Hz). The resistivity is frequency-independent between 19 and 1000 Hz and exhibits metallic character for $`\mathrm{T}>\mathrm{T}_\mathrm{N}`$, with a residual-resistivity ratio $`\rho `$ (300K)/$`\rho `$ (4.2K) $``$ 20. The temperature derivative of the resistivity changes sign at $`\mathrm{T}_\mathrm{N}`$ (insets in Fig. 4). This is the behavior expected if a spin density wave opens a small gap at the Fermi surface.

In conclusion, the thermodynamic, transport, and $`\mu `$SR data reported here are consistent with an intrinsic, commensurate antiferromagnetic spin density wave transition at 19.8 K, but incompatible with incommensurate magnetic order sugiyama03 . Our data may prove valuable in elucidating the origin of the large thermopower of $`\mathrm{Na}_\mathrm{x}\mathrm{CoO}_2`$ and of the unconventional superconducting state in the closely related compound $`\mathrm{Na}_{0.35}\mathrm{CoO}_21.3\mathrm{H}_2\mathrm{O}`$.

We acknowledge the technical support of A. Amato at the PSI, and E. Brücher and G. Siegle at the MPI.


text:
# Testing Bimodal Planet Formation

## 1. Test of Bimodality

Fischer, Valenti & Marcy (2003) have demonstrated that the frequency of extra-solar giant planets is a strong function of metallicity $`Z`$. Santos, Israelian & Mayor (2004) have confirmed this result and have further shown that there is a sharp break in frequency at $`Z=0.02`$, which can be represented algebraically by,

$$f=0.025+16(Z0.02)\mathrm{\Theta }(Z0.02),$$
(1)
where $`\mathrm{\Theta }`$ is a step function. See their figure 7. We suggest that the two terms of this equation correspond to two different modes of giant-planet formation, the first being metallicity-independent and the second being strongly dependent on metallicity. For example, in its simplest form, the disk-instability model of Boss (1995) would not be expected to depend on metallicity. By contrast, the more conventional picture of gas accretion onto rock-ice cores might well be very sensitive to metallicity.

If the two mechanisms generated planets with substantially different distributions in eccentricity, period, or mass, then these should be revealed in the observed distributions of these properties as functions of metallicity. That is, the planets with $`Z<0.02`$ should entirely reflect one distribution, while those with $`Z0.02`$ should predominantly reflect the other.

To test this possibility, we make three plots. In Figures 1, 2, and 3, we show the eccentricities $`e`$, periods $`P`$, and mass functions $`M\mathrm{sin}i`$ of the 109 planets cataloged by Santos et al. (2004) versus the metallicities of their 97 host stars. We exclude HD 47536 as its planet’s properties are not well established. We take planet properties from The Geneva Extrasolar Planet Search Programmes website111http://obswww.unige.ch/$``$udry/planet/planet.html.

No strong pattern emerges from any of these three plots. The low-$`Z`$ stars show a modest deficit of short-period planets, but a Kolmogorov-Smirnov (KS) test shows a significance of only $`p=0.17`$. Similarly, the low-$`Z`$ stars show a modest deficit of high-mass planets, but this is even less significant ($`p=0.21`$). However, because only 23 of the 109 planets have metallicities $`Z<0.02`$, these plots and tests would only be sensitive to the most glaring differences between the two putative populations.

Since there are marginal signals in each of Figures 2 and 3, we have also plotted (but do not show) $`P`$ vs. $`M\mathrm{sin}i`$, with different symbols for the planets of stars above and below $`Z=0.02`$. However, we find no enhanced signal in this plot, nor indeed in similar plots of $`e`$ vs. $`P`$ and $`e`$ vs. $`M\mathrm{sin}i`$.

We can also check whether planetary multiplicity is bimodal in character. Of the 75 planet-bearing stars with $`Z0.02`$, 9 have more than one planet. If multiplicity were independent of metallicity, we would expect 2.64 of the 22 planet-bearing stars with $`Z<0.02`$ to have multiple planets. In fact, there is only one such metal-poor multiple system. Again, however, no strong significance can be attached to this difference because, based on Poisson statistics, such a shortfall would be expected 26% of the time.

Now that a clear break point has been established in planet frequency at $`Z=0.02`$, it is critical to establish a larger sample of planets of sub-solar metallicity stars. The relatively low frequency of these planets will make this search difficult, partly because a large number of such stars must be searched and partly because, as a consequence of this, these stars will tend to be farther away and so fainter. Nevertheless, since a comparison of the planet properties of host stars above and below $`Z=0.02`$ is the simplest test for bimodal planet formation, the additional effort required to enhance the low-metallicity sample would be well justified.

This work was supported by grant AST 02-01266 from the NSF.


text:
# Controlling the coexistence of structural phases and the optical properties of gallium nanoparticles with optical excitation

## Abstract

We have observed reversible structural transformations, induced by optical excitation at 1.55 $`\mu m`$, between the $`\beta `$, $`\gamma `$ and liquid phases of gallium in self-assembled gallium nanoparticles, with a narrow size distribution around 50 nm, on the tip of an optical fiber. Only a few tens of nanowatts of optical excitation per particle are required to control the transformations, which take the form of a dynamic phase coexistence and are accompanied by substantial changes in the optical properties of the nanoparticle film. The time needed to achieve phase equilibrium is in the microsecond range, and increases critically at the transition temperature.

We report that the coexistence of different crystalline and disordered phases in nanoparticles can be controlled by optical excitation in a *controllable*, *continuous* and *reversible* fashion. This has been observed in nanoparticles of polymorphic elemental gallium. Shifts in the balance between phases are accompanied substantial changes in the optical properties of the nanoparticle film. The light-induced transformations display phenomenological features typical of 1st-order Ehrenfest phase transitions (such as a reflectivity hysteresis) but simultaneously display characteristics of 2nd-order transitions (such as critical dependencies on temperature of the ‘susceptibility’ and relaxation time of the stimulated response). Our study is motivated by a desire to understand the exciting physics of phase equilibria in nanoparticles Smirnov (1994); Berry et al. (1984); Berry and Smirnov (2000) and in particular in metallic nanoparticles, which have the potential to play a key role in future highly integrated photonic devices as the active elements of waveguiding Krenn et al. (1999) and switching Zheludev (2002) structures.

We studied light-induced structural transformations in gallium nanoparticles by monitoring the optical reflectivity of nanoparticle films. The particles, typically 50 nm in diameter with a relatively narrow size distribution ($`\pm 14`$ $`nm`$), were prepared on the tips of silica optical fibers, using the recently developed light-assisted self-assembly technique MacDonald et al. (2002); Fedotov et al. (2003a). This process produced a nanoparticle film on the fiber’s core ($`9`$ $`\mu m`$ in diameter) comprising $``$$`2.10^4`$ nanoparticles. Phase transitions in the nanoparticles ware stimulated by a $`1.4`$ $`mW`$, $`1.55`$ $`\mu m`$ diode laser launched into the fiber. Its output was modulated with $`50\%`$ duty cycle at $`6.2`$ $`kHz`$ providing an average excitation power of about 22 $`nW`$ per nanoparticle. Another $`0.4`$ $`mW`$ cw diode laser operating at $`1.31`$ $`\mu m`$ and phase-sensitive detection apparatus were used to monitor the reflectivity of the film. The detection system had an overall bandwidth of $`300`$ $`kHz`$.

The nanoparticle film’s reflectivity shows a very wide ($`>`$$`100`$ $`K`$) hysteresis, with two distinct steps, at $`T_0^{}233`$ $`K`$ and $`T_0253`$ $`K`$, in the rising temperature part of the cycle (Fig. 1a). The presence of these steps indicates that there are two structural phase transitions and therefore that either the nanoparticles undergo a transition from one solid phase to another and then from that phase to the liquid, or that two solid phases coexist at low temperatures and undergo the transition to liquid at different temperatures. Analysis of gallium’s phase diagram Bosio (1978) and results from energy-dispersive x-ray diffraction studies of gallium nanoparticles Di Cicco (1998) suggest that in the present case the solid phase with the highest melting point is $`\beta `$-gallium, and that the other solid phase is $`\gamma `$-gallium.

The presence of pump excitation at $`1.55`$ $`\mu m`$ changes the reflectivity of the nanoparticle film by several percent in a reversible and reproducible fashion. The magnitude $`\mathrm{\Delta }`$ of the induced reflectivity change (Fig. 1b) has a non-zero value at all temperatures, with pronounced peaks that coincide with the steepest parts of the corresponding reflectivity versus temperature plot. With increasing temperature, and the relative phase shift (Fig. 1c) between the pump modulation and the change in probe reflectivity $`\varphi `$ $`0`$ across the whole range from $`90`$ $`K`$ to $``$$`253`$ $`K`$ $`(T_0)`$, indicating that the reflectivity responds quickly to optical stimulation, with relaxation times smaller than $`10`$ $`\mu s`$. At $`T_0`$ and $`T_0^{}`$, $`\mathrm{\Delta }`$ increases by a factor of $``$$`8`$, and above $`T_0`$ the phase shift $`\varphi `$ becomes negative. Above $``$$`255`$ $`K`$, $`\varphi 180^{}`$, meaning that optical excitation actually reduces the reflectivity of the nanoparticles in the molten phase. The relaxation time $`\tau `$ of the nonlinearity increases steeply in the proximity of the peak with scaling dependence $`\tau (T_0T)^{1.7}`$, while the magnitude of the response in the same temperature range follows $`\mathrm{\Delta }(T_0T)^{1.5}`$ (Fig. 2). Remarkably, the phase shift between the pump modulation and the probe reflectivity change remains almost zero at $`T_0^{}`$, indicating that the response and relaxation times remain short ($`<10`$ $`\mu s`$). With decreasing temperature, the pump and probe remain in anti-phase until $``$$`180`$ $`K`$, then $`\varphi `$ returns to zero by $``$$`150`$ $`K`$. A small, smooth peak is seen in $`\mathrm{\Delta }`$ at the re-crystallization temperature $`T_0^{}`$.

The phenomenological picture seen here is unusual and is quite different from what is normally observed in the bulk where first order phase transitions are characterized by a discontinuous change in the state of the body, effected by means of a sudden rearrangement of the crystalline lattice. Many features of our experimental results may be explained by consideration of the fact that under the confined conditions of the nanoparticle geometry, there is a continuous dynamic coexistence of structural forms in the nanoparticles. Such behavior has been studied theoretically for small clusters, and it was predicted that there would be a temperature interval wherein two different phases coexist Berry et al. (1984); Smirnov (1994). The surface of a particle, where atoms have fewer nearest neighbors than internal atoms, becomes very important and acts as a boundary at which transformation processes can start. We believe that this concept is largely applicable to the present case and it brings us to a model in which we consider an ellipsoidal nanoparticle with a core in one structural phase covered by a homogeneous ‘shell’ of another phase. In the simplest case of a phase transition to the melt, the light-induced behavior of such a particle would be analogous to the temperature-driven ‘surface melting’ effect that has been seen in lead nanoparticles Garrigos et al. (1989) and found to be thermodynamically reversible within a narrow temperature range Peters et al. (1997). In the presence of light, the equilibrium position of the phase boundary, or to put it another way, the thickness of the surface layer, will be determined by both temperature and light intensity. At a certain temperature $`T_a`$ below the phase transition point $`T_0^{^{}}`$ the influence of light on the surface layer’s thickness becomes apparent as optical excitation changes the reflectivity of the film (see Fig. 1b). With increasing temperature or level of optical excitation, the surface layer’s thickness increases until the transformation of the core to the ‘surface’ phase is completed. When, at $`T_b`$, the core of the particle is fully consumed by the new phase the nanoparticle becomes stable against a return to the old phase because this would require the creation of a nucleation center. However, if the temperature or level of optical excitation is reduced *before* the transformation to the new phase is complete, i.e. while a nucleus of the old core phase is still present, the transformation is reversed and the skin layer shrinks to an appropriate equilibrium position. Thus, reversibility is provided in the temperature range between $`T_a`$ and $`T_b`$. This whole process is then replicated around the next phase transition temperature $`T_0`$ (see Fig. 3a). From the standpoint of phase transition theory the observed effect shows some features of second-order transitions: the thickness of the surface layer is characterized by the correlation length of near surface fluctuations $`\xi (T_0T)^\nu `$. The susceptibility is then a scaling function of the correlation length: $`\chi =(\delta \eta /\delta h)_T(T_0T)^{\gamma _1}`$, where $`\gamma _1=1.5`$ in our experiment. When the optical excitation is terminated, the skin layer retreats and the reflectivity is restored to its original level. Like the susceptibility, the relaxation time has scale-invariant form, i.e. $`\tau (T_0T)^{\gamma _3}`$, where $`\gamma _3=1.7`$ in our experiment. The longer relaxation times near the phase transition temperatures result from the increased thickness of the surface layer and reduced velocity of the interface between the phases. In terms of Landau theory this situation may be described by taking the free energy of a particle to be a polynomial function of the order parameter $`\psi r^{1/3}`$, where $`r`$ is the radius of the solid core within a particle of radius $`R`$. As is usual in the vicinity of a second order phase transition, the value of the order parameter will characteristically be proportional to $`(T_0T)^{1/2}`$, giving $`r(T_0T)^{1.5}`$. The reflectivity of the nanoparticle film will then, to a first approximation, be proportional to the volume of the molten phase: $`V_m=4\pi R^3[(Rr)/R((Rr)/R)^2+((Rr)/R)^3/3]t^{1.5}`$, where $`t=(TT_0)/(T_0)`$, thus producing correct estimates for the experimental values of the critical indices. Intriguingly, the relaxation time for the first transition at $`T_0^{^{}}`$ is much shorter than for the second one at $`T_0`$. This may be an indication that the first transition occurs between two solid phases and would be consistent with the thermodynamic argument that the enthalpic and entropic contributions to the Gibbs free energy appear to enter differently into the kinetics of solid-solid and solid-liquid transitions because of the difference in the configurational and vibrational space available in the solid and liquid phases Tsao et al. .

The strength of the phase coexistence concept is supported by calculations of the optical properties of gallium nanoparticle films on a dielectric substrate performed using a recently developed effective-medium model for densely packaged nanoshells Fedotov et al. (2003b). For the purposes of our calculations, the dielectric constants of $`\beta `$ and $`\gamma `$ gallium, which are much closer to those of a free-electron metal than those of the $`\alpha `$ phase, were estimated by using the damping constant in Drude’s free-electron model as a fitting parameter to produce the nanoparticle film reflectivity levels shown in Fig. 2a. These calculations (also detailed in Ref. Fedotov et al. (2003b))confirmed that the presence on each nanoparticle of a shell just a few nanometres thick in a phase different from the core can produce a change in reflectivity sufficient to explain our experimental data.

A thermally activated transition due to laser-induced heating can explain certain characteristics of the effect. For instance, by assuming a local light-induced temperature increase of $`0.42`$ $`K`$, one can derive a good facsimile of the experimental peaks in $`\mathrm{\Delta }`$ at $`T_0^{^{}}`$ and $`T_0`$ from the reflectivity data in Fig. 1a. However, there are serious discrepancies between the results of this thermal model and the experimental results, primarily at temperatures more than a few degrees below the peaks, where the observed effect is larger than predicted. The thermal model also fails to explain important details of $`\varphi `$’s dependence on temperature. This suggests that another, temperature-independent, non-thermal excitation mechanism is also contributing to the effect. This mechanism may be especially important for gallium which is know to have covalent bonds within some of its crystalline structures Bernasconi et al. (1995). The covalent bonds’ absorption line encompasses the pump wavelength Gong et al. (1991), so optical excitation may result in bonding-antibonding transitions, which destabilize the structure Siegal et al. (1995). A ‘defect’ or ‘inclusion’ of a new phase is thus created, changing the optical properties of the ‘host’ phase at temperatures far below its transition point. Furthermore, the latent heat of transition from the defect-containing structure to a new phase is reduced in proportion to the number of defects. The latent heat controls the balance between phases Trittibach et al. (1994), so this defect creation eventually shifts the phase equilibrium and promotes the formation of a thicker layer of the new phase without any increase in temperature. This process could provide an additional, non-thermal contribution the peak at $`T_0`$ and $`T_0^{^{}}`$. Such ‘optical melting’ mechanisms, proceeding through light-induced destabilization of the crystalline bond structure, have been analysed previously for gallium MacDonald et al. (2001) and selenium Poborchii et al. (1999).

In conclusion, we have found that light can stimulate reversible structural transitions in gallium nanoparticles that belong to a novel class of surface-driven excitation-induced phase transitions and have some striking similarities with a transition recently observed in superfluid helium droplets Kimball et al. (2002) and a transformation believed to occur in confined superconductors, wherein the existence of a ‘boundary phase’, characterized by spontaneous fluctuational nucleation of short-lived vortices is critical Sobnack and Kusmartsev (2001). If there are substantial differences between the optical properties of the phases involved, as there are in gallium, this type of transition provides a new means of obtaining a large optical nonlinearity, which may be used to control light with light in nanostructured materials such as ‘plasmon waveguides’ made from arrays of metallic nanoparticles.

The authors acknowledge the financial support of the EPSRC (UK).


text:
# The Formation of Self-Gravitating Cores in Turbulent Magnetized Clouds

## 1 Introduction

In the standard theory of isolated star formation, the presence of magnetic fields in molecular clouds plays a vital role. Based on the Jeans argument, observed high-density molecular clouds should collapse within a few freefall times to form stars. Instead, the molecular clouds were thought to be quasi-stable. Turbulence in molecular clouds and magnetic field support were proposed to explain these apparently stable clouds (e.g. McKee, 1999; Williams, Blitz, & McKee, 2000; Vázquez-Semadeni et al., 2000). Molecular clouds were suggested to be supported magnetostatically (e.g. Mouschovias & Spitzer, 1976; Fiege & Pudritz, 1999) or dynamically by magnetohydrodynamic (MHD) waves, especially Alfvén waves (e.g. Dewar, 1970; Shu, Adams, & Lizano, 1987). Alfvén waves are transverse waves and produce perturbation perpendicular to the mean magnetic field. McKee & Zweibel (1995) pointed out that Alfvén waves can lead to an isotropic turbulence pressure to counteract gravitational collapse, provided that the waves are neither damped nor driven. However, this mechanism requires a negative radial gradient in wave sources in the cloud in order to support the cloud from collapsing (Shu, Adams, & Lizano, 1987).

Observationally, the question of whether magnetic fields are sufficiently strong to support molecular clouds alone from collapsing remains unresolved. Crutcher (1999) concludes that static magnetic fields are insufficient to support the observed clouds and that MHD waves are equally important in cloud energetics. McKee (1999) however points out that ambipolar diffusion might already have altered the mass-to-flux ratio observed by Crutcher. Observations by Bourke et al. (2001) yield similar conclusions to Crutcher (1999), although they caution that the measured magnetic field strength may be significantly lower than the true values due to the low volume filling fraction of dense cores. Nakano (1998) argues that observed cloud cores cannot be magnetically supported. A magnetically supported core should not have a column density much higher than its surroundings, and the core cannot maintain a large velocity dispersion inside it for its whole lifetime. However, observed cloud cores generally have superthermal velocity dispersions, and column densities much higher than their surroundings, suggesting that most cloud cores are not magnetostatically supported.

Significant progress in the understanding of the roles of turbulence and MHD waves in supporting gravitationally unstable cloud cores has been made recently using 3D numerical simulations (e.g. Vázquez-Semadeni, Passot, & Pouquet, 1996; Mac Low et al., 1998; Padoan & Nordlund, 1999a; Ostriker, Gammie, & Stone, 1999, for a review, see Mac Low & Klessen 2004). Klessen, Heitsch, & Mac Low (2000, hereafter Paper I), using models computed with both ZEUS-3D and a smoothed particle hydrodynamics (SPH) code, conclude that hydrodynamical turbulence can prevent global collapse but not local collapse under typical molecular cloud conditions. Heitsch, Mac Low, & Klessen (2001, hereafter Paper II) include magnetic fields and conclude that local collapse cannot be prevented much longer than a global free-fall time by magnetized turbulence, in the absence of mean-field support. However, they find that magnetic fields delay local collapse by decreasing local density enhancements via magnetic pressure. Stars begin to form quickly when local density enhancements collapse. This result favors the dynamical picture of molecular clouds being a transient feature in the interstellar medium (Ballesteros-Paredes et al., 1999; Elmegreen, 2000; Hartmann, Ballesteros-Paredes & Bergin, 2001) instead of existing for many free-fall times.

Another interesting result from the 2563 resolution simulations of Papers I and II is that several cores begin to evolve into flattened (or disk-like) objects, which leads to the speculation that they were beginning to resolve the formation of protostellar cores. Unfortunately, the resolution and number of disk-like cores were not high enough to perform a statistically meaningful study. Therefore, these flattened cores were not discussed in Paper II.

Simulations suggest that supersonic turbulent flows, driven at large scale by supernovae and density waves, produce shock-compressed gas sheets or filaments that fragment into dense cores, and drive the observed supersonic turbulence in the clouds (Mac Low & Klessen, 2004). Differential rotation of galactic disks (Fleck, 1981) provides another plausible driving mechanism to maintain interstellar turbulence, especially at galactic radii with little massive star formation. The coupling from large-scale shear down to turbulent scales could be through the magnetorotational instability (e.g. Balbus & Hawley, 1998; Sellwood & Balbus, 1999; Kim, Ostriker, & Stone, 2003). Ionizing radiation (McKee, 1989; Vázquez-Semadeni, Passot, & Pouquet, 1995; Bertoldi & McKee, 1997; Kritsuk & Norman, 2002), and stellar winds from massive stars may also be important, but supernovae (e.g. Norman & Ferrara, 1996; Avillez, 2000) probably provide most of the energy required to maintain interstellar turbulence as well as the turbulence within molecular clouds (Mac Low & Klessen, 2004).

Recent numerical studies indicate that supersonic turbulence may play an important role in shaping the initial mass function (IMF) of the stellar population that forms in molecular clouds (Padoan & Nordlund, 2002). Observations of the IMF appear consistent with a general form having a power law in the high-mass wing and then flattening and turning over at the low mass end (e.g. Kroupa, 2002, and references therein). This is an important constraint on star formation theories. Interestingly, and perhaps significantly, the probability distribution function (PDF) of gas density in isothermal supersonic turbulence is well approximated by a log normal distribution. Padoan, Nordlund, & Jones (1997) and Padoan & Nordlund (2002) suggest that the IMF and density PDF in turbulent clouds are intimately related, and predict that the mass distribution of dense cores is controlled by the density, the Mach number, and velocity dispersion of the supersonic turbulence. Previous simulations of molecular cloud collapse (e.g. Padoan et al., 2001; Klessen, 2001; Bate, Bonnell, & Bromm, 2003; Gammie et al., 2003) indicate that the mass distribution of cores (or clumps) roughly resembles a log normal distribution, but with the high-mass wing following a power law approximately consistent with the Salpeter law. However, these simulations either lack magnetic fields or do not have a sufficient number of gravitationally bound cores for accurate statistical analysis, because of low numerical resolution. Higher resolution MHD turbulence simulations are required to verify the importance of super-Alfvénic turbulent fragmentation to the stellar IMF.

In this paper, we report on results of our latest three-dimensional (3D), turbulent simulation of a magnetized molecular cloud, increasing the spatial resolution in each dimension by a factor of two relative to the simulations in Paper II, using the ZEUS-MP code. The excellent parallel scalability of ZEUS-MP allows us to study MHD turbulent molecular clouds on a 5123 grid, sufficient to study the global physical properties of cores and their formation process. In §2, we describe the techniques, assumptions, and parameters used in our simulations. In §3, we discuss the statistical properties of the super-Alfvénic turbulent cloud before gravitation is turned on. The global physical properties of self-gravitating cores and their time evolution are reported in §4. We also examine the turbulent fragmentation theory by comparing its predictions to our simulation results. In §5, we perform a resolution study and investigate the physical properties of some disk-like cores observed at later stages of our simulation. We discuss our results and present conclusions in §6.

## 2 Simulations

Our simulations are performed on 256 SGI/Cray Origin 2000 processors at the National Center for Supercomputing Applications using the ZEUS-MP code. ZEUS-MP is a multi-physics, massively-parallel, message-passing code for astrophysical fluid dynamic simulations in 3D, developed by the Laboratory for Computational Astrophysics (lca.ncsa.uiuc.edu) at NCSA. ZEUS-MP is a distributed memory version of the shared-memory code ZEUS-3D, which uses block domain decomposition to achieve scalable parallelism (Norman, 2000). The code includes ideal hydrodynamics, ideal MHD, and self-gravity. Self-gravity is implemented using the FFTW library (Frigo & Johnson, 1998).

Our scale-free MHD turbulent simulations are basically the same as model Eh1w in Paper II, but performed at resolutions up to 5123. All parameters are given in normalized units, where physical constants are scaled to unity. The total mass in the box $`M`$ = 1 and the side length $`L`$ = 2. The average density $`\rho _o`$ = 0.125. The initial conditions of our simulation are identical to model Eh1w in Paper II. In this model, the sound speed $`c_s`$ = 0.1. The corresponding Jeans length $`\lambda _\mathrm{J}`$ = 0.5, and Jeans mass

$$\mathrm{M}_\mathrm{J}\rho _o\lambda _\mathrm{J}^3=\left(\frac{\pi }{G}\right)^{3/2}\rho _o^{1/2}c_s^3=0.0156.$$
(1)
Therefore, the mass in the box $`M`$ = 64 MJ and the length of the box $`L`$ = 4 $`\lambda _\mathrm{J}`$. The magnetic field is initially uniform and along the z-direction with strength $`B_o`$ = 0.188. The ratio $`M/M_{\mathrm{cr}}`$ = 8.3, which means our entire cloud is magnetically supercritical (see e.g. Mouschovias & Spitzer, 1976). The ratio of thermal to magnetic pressure is 0.9 and the rms Mach number for this parameter set is $`M_{\mathrm{rms}}`$ = 10. Isothermality is assumed throughout the simulation, which is justified as the cooling time is much shorter than the dynamical time inside high density molecular clouds. No ambipolar diffusion is included in the calculation.

Periodic boundary conditions are applied in all three dimensions. In effect we simulate a small portion of a larger cloud. We do this for simplicity, to avoid addressing the structure of the non-isothermal boundaries of molecular clouds. We do not expect that our results will be contradicted by unified models that include larger-scale flows in the interstellar medium. In our analysis, we examine column-densities through one realization of our periodic cube, which should give statistically valid answers to the questions we ask for molecular clouds of finite size.

The procedure of setting up the turbulent flow is explained in Mac Low (1999). We briefly review it here. We assume that each component of the velocity perturbation is a Gaussian random field with flat power spectrum in the range 1 $`|k|`$ 2. Random phases and amplitudes are generated in that spherical shell in Fourier space, and then transformed back into real space to generate each component of the driving velocity perturbation. When all three components are obtained, the amplitude of the velocity is scaled to a desired initial root mean square velocity, which is unity in our case. This velocity field is then normalized at each time step to ensure a constant kinetic energy input rate and added in as a perturbation to the evolved velocity field. The dynamical behavior of isothermal self-gravitating gas is scale free. Equations (2) to (4) can be used to scale back to astrophysical units with a mass scale of the thermal Jeans mass, length scale given by the thermal Jeans length, and a time scale given by the free-fall timescale (see Paper I),

$$L=0.89\mathrm{pc}\left(\frac{c_s}{0.2\mathrm{kms}^1}\right)\left(\frac{n}{10^4\mathrm{cm}^3}\right)^{1/2}$$
(2)
$$M=413\mathrm{M}_{\mathrm{}}\left(\frac{c_s}{0.2\mathrm{kms}^1}\right)^3\left(\frac{n}{10^4\mathrm{cm}^3}\right)^{1/2}$$
(3)
$$\tau _{\mathrm{ff}}=\left(\frac{3\pi }{32G\rho }\right)^{1/2}=0.34\mathrm{Myr}\left(\frac{n}{10^4\mathrm{cm}^3}\right)^{1/2}$$
(4)
Therefore, if we choose a high-density scaling such as the BN region in Orion, in which $`c_s`$ = 0.2 km s-1 and the number density $`n`$ = 105 cm-3, the region in the simulation box is 0.28 pc in size. The total mass in this region is 130.6 M and one system time unit, which is 0.65 free-fall time, is 0.07 Myr.

A total of 110,000 CPU hours were required for the 5123 zone simulation to run to 4.325 time units, or $``$ 2.82 $`\tau _{\mathrm{ff}}`$. In our case, $`\tau _{\mathrm{ff}}`$ = 1.53 system time units. At the beginning of the simulation, no gravitational forces are applied. Turbulence is allowed to develop from the relatively smooth initial conditions until $`t`$ = 2.0 time units. Gravitation is then turned on. Density fluctuations generated by the supersonic turbulence in converging and interacting shock fronts that locally exceed the Jeans limit begin to contract. In the 2563 simulations reported in Paper II, there were only a dozen gravitationally collapsing cores, which were under-resolved and insufficient for statistical analysis. In our 5123 simulation, there are 15 cores that satisfy our definition of gravitationally bound at $`t`$ = 2.0 when gravity is turned on. By the end of the simulation at $`t`$ = 4.325, we have 83 gravitationally bound cores. The simulation was terminated then because too many cells violate the Jeans condition (Truelove et al., 1997) which dictates that the local Jeans length must be resolved by several zones. Also, the densities of some cores became so high that the isothermal assumption would be invalid. Henceforth, we adopt the same time convention as in Paper II, where the time $`t`$ = 0 is defined as the time that gravity is turned on. According to this definition, the simulation is terminated at $`t`$ = 2.325 system time, or 1.52 $`\tau _{\mathrm{ff}}`$. We will use free-fall time as the time unit for the rest of the paper.

## 3 Statistical Properties

### 3.1 Power Spectra

Kolmogorov (1941) heuristically derived a simple scaling for the 1D energy spectrum in incompressible turbulence, $`E(k)k^{5/3}`$, where $`k`$ is the magnitude of the wavenumber. Astrophysical fluids are magnetized and highly compressible. Dynamically important magnetic fields and strong shocks should interfere with eddy motions, which in turn will affect the final velocity power spectrum. However, observations and some simulations (e.g. Paper I; Lazarian, Pogosyan, & Esquivel, 2002; Cho, Lazarian, & Vishniac, 2003) suggest that MHD compressible turbulence also closely follows the Kolmogorov spectrum. The theory of Goldreich & Sridhar (1995) points out that it is the perpendicular motion of sub-Alfvénic incompressible turbulence that behaves as Kolmogorov, $`E(k)k^{5/3}`$. Some other MHD turbulence simulations suggest that the scaling law for MHD compressible turbulence is not exactly Kolmogorov. Boldyrev, Nordlund, & Padoan (2002) reported that the power spectrum of the solenoidal component of the velocity field in a super-Alfvénic compressible turbulence is $`E(k)k^{1.74}`$, steeper than the Kolmogorov spectrum. In the latest simulations by Vestuto, Ostriker, & Stone (2003) on 2563 and 5123 grids, the slope of the power spectra varies from Kolmogorov to even larger than the $`E(k)k^2`$ Burgers’ model of shock-dominated turbulence (Burgers, 1974) for models ranging from sub-Alfvénic to super-Alfvénic. The measurement of the slopes in their simulations might be affected by the short (or nonexistent!) inertial range of their power spectra, however.

Figure 1a shows the power spectra of density, velocity, and the energies of our turbulent cube at $`t`$ = 0, after the turbulence has fully developed, but before gravity has been turned on. The long dashed line is the driving spectrum. The power spectrum of velocity averaged over all directions is basically a Kolmogorov power spectrum (the thin straight line in Figure 1a). The inertial range of the velocity power spectrum spans more than an order of magnitude of the wavenumber (2 $`|k|`$ 30). The power spectrum of the magnetic field shows a similar result. The density spectrum has a shallower slope with index $``$ -1.15. Both the kinetic and potential energies have shallower slopes because of the modulation by the density spectrum. Figure 1b is a compensated plot of the velocity power spectrum multiplied by $`k^{5/3}`$ that highlights the excellent fit of the Kolmogorov power spectrum to the velocity power spectrum. We also plot the compensated velocity power spectrum for the 2563 resolution simulation for comparison. The inertial range of the turbulent flow in this lower resolution simulation is shorter but still close to the Kolmogorov power spectrum.

### 3.2 Density PDF and Core Mass Spectrum

The PDF provides important and complementary information to the power spectrum. The PDF of isothermal turbulent flows can be approximated by a log normal function (e.g. Vázquez-Semadeni, 1994; Padoan, Nordlund, & Jones, 1997; Scalo et al., 1998; Nordlund & Padoan, 1999). With a log normal PDF of mass density, most of the material concentrates in a small fraction of the total volume in the simulation. Figure 2 shows the best log normal fit to the PDFs of all the cells in the 2563 and 5123 simulations at $`t`$ = 0. The standard deviation of the best fit is 1.46 $`\pm `$ 0.06 and 1.42 $`\pm `$ 0.03, as defined in equation (6) in Padoan & Nordlund (2002). Using equation (9) in Padoan & Nordlund (2002), $`\sigma ^2`$ = ln(1 + $`M_{\mathrm{A},\mathrm{rms}}^2\eta ^2`$), we calculate $`\eta `$ to be 0.41 and 0.38, respectively, a little smaller than the value $`\eta `$ 0.5 reported in Padoan & Nordlund (2002) from numerical experiments.

Figure 3 shows the projected density (column density) images of the simulation, in logarithmic scale, along the x-, y-, and z-direction, which shows the clumpiness expected from a log normal PDF. The first row of plots are at $`t`$ = 0, just before gravitation is switched on. The second and third rows are at $`t`$ = 0.65 $`\tau _{\mathrm{ff}}`$ and $`t`$ = 1.31 $`\tau _{\mathrm{ff}}`$, respectively. At $`t`$ = 0 (Figures 3a-c), filamentary structures dominate the view and we see some higher density nodes on the filaments. Later (Figures 3d-f), gravity breaks the filaments into pieces and higher density clumps slowly increase their mass by accretion. The peak density at this moment is about 1000 times the average density. At $`t=1.31\tau _{\mathrm{ff}}`$ (Figures 3g-i), many cores with high column-densities have formed. The peak density is about four orders of magnitude above the average density and the total mass in cores is less than 9% of the total mass in the box. Qualitatively, this picture does not change during the subsequent simulation time (up to $`t=1.52\tau _{\mathrm{ff}}`$), except that more cores form and their peak densities grow. Most of the volume in the simulation corresponds to low density voids. The clumpiness can also be seen in Figure 4 which shows a volume rendering of the logarithmic density and the magnetic pressure of the turbulent box in 3D at $`t`$ = 1.31 $`\tau _{\mathrm{ff}}`$. An animation of the evolution of density and magnetic pressure of the 5123 simulation is included in the electronic version of this paper.

Padoan & Nordlund (2002) propose that the stellar IMF is a direct result of turbulent fragmentation of molecular clouds. As shown in Figure 3, the gas is compressed into 1D filaments and 2D sheets by supersonic turbulence. Therefore, the core mass distribution will be dependent on the jump conditions for isothermal shocks in a magnetized gas. The typical size of a dense core is comparable to the thickness of the post shock gas. Based on this, Padoan & Nordlund (2002) calculate by taking into account the scale-dependent Mach number that the mass distribution of cores formed inside a supersonic MHD turbulent cloud is

$$N(m)dlog(m)m^{3/(4\beta )}dlog(m)$$
(5)
where $`\beta `$ is the spectral power index. With our 5123 simulation, we have sufficient resolution to check the turbulent fragmentation prediction. We proceed as follows.

The clumpy structure inside the turbulent molecular cloud in our simulation is determined using the algorithm CLUMPFIND described by Williams, De Geus, & Blitz (1994). The algorithm is also described in Paper I. A “core” is defined as a gravitationally bound region. We define it to consist of a set of connected zones with average density exceeding the mean density expected for isothermal shocks, $`\rho >M_{\mathrm{rms}}^2\rho _\mathrm{o}`$, potential energy exceeding kinetic energy, $`|E_{pot}^{cell}|>E_{kin}^{cell}`$, and core mass exceeding the local Jeans mass, $`M_{core}>M_\mathrm{J}(\rho )`$. We use logarithmic density contours in CLUMPFIND to get a wide enough density range, as in Paper I. We believe the core definition chosen here is more rigorous than the simple gravitationally unstable overdensity criterion of Padoan et al. (2001). However, our definition identifies fewer cores. This definition of core will be applied to all results on cores in this paper. We also define “clumps” to be overdense regions that do not satisfy the other criteria. We choose the overdensity threshold for clumps somewhat arbitrarily to be $`\rho >0.1M_{\mathrm{rms}}^2\rho _o`$

In Figure 5, we plot the mass spectra of cores (dot-dashed line) and clumps (dashed line) in the 2563 simulation and the spectrum of cores (solid line) in 5123 simulation at $`t`$ = 0, normalized by the Jeans mass $`M_\mathrm{J}(\rho _o)`$. There are only 11 cores in the 2563 simulation. Therefore, we also plot the clump spectrum for comparison. Many of these clumps could be just density fluctuation created by the isothermal shocks that will be destroyed later. The 2563 clump spectrum and 5123 core spectrum both show a well-defined “universal” IMF (see the latest review by Kroupa, 2002). From Figures 1a and 1b, the velocity power spectral index $`\beta `$ -5/3. Substituting $`\beta `$ into equation (5), we expect the slope of the high-mass wing of the core mass spectrum will be -1.29. This spectral index is plotted as thin line in Figure 5 and is consistent with the high-density wings of all the spectra, as predicted by the turbulent fragmentation. Even the 2563 core spectrum appears consistent with the prediction. Note that the number of cores in Figure 5 is still too few to have an accurate fitting unless we relax the definitions of the cores to clumps.

## 4 Cloud Core Evolution

### 4.1 Delay of Core Collapse

Paper I and Paper II concluded that local collapse could not be prevented even when the turbulent velocity field carries enough energy to balance the gravitational contraction on global scales. Paper II found that the comparison between 1283 and 2563 MHD results did not show full convergence, with increasing resolution delaying local collapse (see Figure 6 in Paper II). This occurs because increasing the resolution produces thinner shocks with peak densities closer to the theoretical value. The resulting deeper potential wells accrete more mass and collapse more quickly. At the same time, the simulation can better follow the short-wavelength MHD Alfvén waves, which could create an isotropic pressure to counteract the gravitational collapse. From the 2563 simulations in Paper II, it appears that the MHD waves are able to delay the collapse but not to prevent it. In our new 5123 simulation, this imbalance between gravitation and Alfvén wave pressure is verified. Figure 6 combines our new results with the previous 643, 1283, and 2563 results. We plot the collapsed mass $`M_{}`$, defined as the sum of the core masses, versus time. Collapse still occurs in all cases. Comparison between the 2563 and 5123 models suggests that the 2563 result of Paper II is, in fact, converged. Remaining discrepancies are produced by the chaotic nature of the flow, as was shown by comparisons of different realizations of the same resolution model in Paper I.

### 4.2 Core Formation and Evolution

In the models, we find that self-gravitating cores form out of initially overdense clumps created by supersonic turbulence. We also find that dense clumps can be destroyed by passing shocks as well as by the merging of clumps to form larger and denser objects. The formation and destruction of cores is thus a complex dynamical process. The core formation and destruction process is shown in Figure 7 in terms of number of cores formed as a function of time. Instead of monotonically increasing, the core number drops sometimes because of merging or destruction.

Figures 8 to 10 show the evolution of the physical properties of cores after gravitation is turned on and before the Jeans criterion Truelove et al. (1997) is violated.

Figure 8a and 8b show the ratio of core mass to the mass that the local magnetic field can support, $`M/M_{\mathrm{cr}}`$. The ratio

$$\frac{M}{M_{\mathrm{cr}}}=\frac{M/\mathrm{\Phi }}{(M/\mathrm{\Phi })_{\mathrm{cr}}}=\frac{N}{FB_{\mathrm{los}}}\frac{\sqrt{G}}{c_\mathrm{\Phi }},$$
(6)
where $`N`$ is the column density; the geometrical correction factor $`F=2`$ for a uniform sphere and 3 for an isothermal disk. $`B_{\mathrm{los}}`$ is the line of sight magnetic field strength; $`G`$ is the gravitational constant; and $`c_\mathrm{\Phi }`$ $``$ 0.12 for a spherical cloud or 0.16 for an isothermal disk (Bourke et al., 2001). In these two plots, we choose $`F=3`$ and $`c_\mathrm{\Phi }=0.16`$ for the extreme case. The ratios exceed an order of magnitude at all time. If we assumed spherical cloud geometry, the ratios would be even higher. Therefore, all the cores are magnetically supercritical during the gravitationally collapsing phase, despite the action of gravitational fragmentation, which could be expected to reduce the mass to flux ratio (Mestel & Spitzer, 1956).

Figures 8c and 8d show the relationship between the central magnetic field strength $`B_c`$ and the central core density $`\rho _c`$ at $`t`$ = 0 and $`t`$ = 0.65 $`\tau _{\mathrm{ff}}`$ respectively. The power law of index 0.5 is also plotted. Interpretations of $`B_c\rho _c^{1/2}`$ can be that: 1) cores are magnetically subcritical (e.g. Mouschovias, 1991; Basu & Mouschovias, 1994, 1995), 2) cores are magnetically supercritical but the Alfvénic speed is about constant ($``$ 1) for the cores (e.g. Bertoldi & McKee, 1992; Crutcher, 1999), 3) cores are magnetically supercritical and dense cores tend to accrete mass along magnetic field lines and reduce their magnetic flux to mass ratio efficiently, even in the absence of ambipolar diffusion (Padoan & Nordlund, 1999b). The cores in the simulations are all magnetically supercritical and there is no ambipolar diffusion implemented in ZEUS-MP. Therefore, Figure 8 could be a result of case 2) or 3). Magnetically subcritical cores are no longer a unique explanation to the $`B_c\rho _c^{1/2}`$ relation.

In Figures 9a and 9b, we plot the specific angular momentum, $`j`$, of each core against the radius, R, of the core at $`t`$ = 0 and 0.65 $`\tau _{\mathrm{ff}}`$. The radius $`R`$ is determined from the maximum dimension of the core determined using CLUMPFIND. The majority of cores obey the power law relationship, $`jR^p`$, with index $`p`$ = 3/2, which is plotted for comparison. From observations, Goodman et al. (1993) concludes that $`j`$ scales roughly as $`R^{3/2}`$ based on the velocity gradients of 29 cores in dark clouds. They interpreted this relationship between specific angular momentum and core size to imply that cores are in approximate ”virial equilibrium”. With the isothermal equation of state and relatively weak magnetic field support in our simulation, the equilibrium in cores is basically between gravitational force and rotational angular momentum, producing rotationally supported cores. In their study of rotational properties of centrally condensed, turbulent, molecular cloud cores, Burkert & Bodenheimer (2000) conclude that the specific angular momentum is related to the radius by $`j\mathrm{\Omega }R^2`$, where $`\mathrm{\Omega }`$ is the angular speed of the core. Therefore, $`j=\mathrm{\Omega }R^2R^{3/2}`$ implies $`\mathrm{\Omega }R^{1/2}`$. This agrees reasonably well with observations that suggest $`\mathrm{\Omega }`$ scales roughly $`R^{0.4}`$ (Goodman et al., 1993; Barranco & Goodman, 1998). However, we cannot conclude just from Figure 9 that the cores are in Keplerian rotation, $`\mathrm{\Omega }(r)r^{0.5}`$, or in angular momentum conserved rotation, $`\mathrm{\Omega }(r)r^1`$. In addition, there is still a slight scattering of data in Figure 9. The rotation could indeed be Keplerian at the inner parts of the core but angular momentum conserved in the outer parts.

In Figure 10a, we plot the relationship between mean density and radius of each core at $`t`$ = 0 (circle) and $`t`$ = 0.65 $`\tau _{\mathrm{ff}}`$ (triangle), where mean density $`<\rho >=M_{core}/V_{core}`$, and $`V_{core}`$ is the volume of the core. The result agrees with that of Ballesteros-Paredes & Mac Low (2002) (see their Figure 9) that mean density is basically independent of radius. In Figure 10b, we plot the mass to radius relationship of each core. A power law relationship is seen, with $`M_{core}R^{2.75}`$ from $`t`$ = 0. Figures 10a and 10b would seem to be inconsistent because we would expect $`M_{core}R^3`$ if cores are spherical. If we have $`V_{core}R^{2.75}`$ from Figure 10, that means either that the geometry of the cores lies between a sphere and a disk, or that the cores are a mixture of prolate, triaxial, and oblate objects, as data points are scattered around the straight line. We investigate the shapes and orientation of cores in §4.4, and confirm this conjecture. The correlations seen in Figures 8 to 10 remain valid for the whole simulation even after Truelove et al. (1997) criterion is violated. They are also found at lower significance in the 2563 resolution simulation.

The number of cores found with CLUMPFIND from $`t`$ = 0.0 to $`t`$ = 2.25 basically increases with time, as shown in Figure 7. Some cores are destroyed by the supersonic turbulence in the simulation. At the same time, some other cores actually merge to form more massive cores through gravitation (e.g. Bonnell et al., 1997; Klessen & Burkert, 2000; Bonnell et al., 2001a, b). Figure 11 shows the merging of two cores, or, to be more specific, the accretion of a small clump by a massive clump. The three columns of subplots in Figure 11 show views along three different axes of the merging sequence. The time interval between two rows is 0.05 $`\tau _{\mathrm{ff}}`$. As the result of core mergers and accretion of surrounding material, the final cores possess much larger angular momentum than the initial clumps. This leads to the formation of disk-like cores during the simulation. This phenomenon is entirely absent in isolated star formation simulations, which consider only a single, gravitationally bound core.

### 4.3 Core Mass Spectrum Evolution

It is interesting to examine the effect of core mergers on the core mass spectrum. In Figure 12, we plot the evolution of the core mass spectrum from $`t`$ = 0 to the end of the simulation. In each panel, we plot a power-law with slope -1.29 for reference. As discussed in §3.2, at $`t`$ = 0, the core mass spectrum resembles closely a log normal distribution and the slope of the high-mass wing is consistent with the power law with index predicted by turbulent fragmentation theory (Padoan & Nordlund, 2002). There are some fluctuations of the core mass spectra at later time, but the slopes of the high-mass wing remain about the same until $`t1\tau _{\mathrm{ff}}`$. The slope becomes shallower at yet later times, but we do not have too much confidence in the correctness of the core mass spectra after $`t=0.65\tau _{\mathrm{ff}}`$, as the centers of the cores become unresolved. The hydrodynamical simulation by Klessen (2001), using SPH and sink particles, also shows shallower slope at the later stages of the simulation, which he attributes to the coalescence of cores, but again, fragmentation in the centers of cores is suppressed, in this case by the use of sink particles. Even if the trend of increasingly shallower slope in core mass spectra is qualitatively correct in our simulation, the rate of slope increase may be over-estimated. Our assumed periodic boundary condition means the cores have no place to disperse except merging in the later stage. The finite resolution in the simulation prevents cores from collapsing as far as real cores would, so that the simulated cores have larger cross sections for merging (Paper I). Nevertheless, for a dense cluster, the core merging rate may still be high. Recent millimeter wave observations of cloud cores also reveal mass spectra of clumps similar to the stellar IMF (e.g Motte, André, & Neri, 1998; Testi & Sargent, 1998), with slope of the high-mass wing of $``$ -1.1 to -1.5, closer to Salpeter -1.35 than that derived for gaseous clumps of -0.5 (Blitz, 1993), using the form of IMF definition in equation (5). Our simulations are in good agreement with these studies.

Because of the higher resolution of our 5123 simulation, we have many more cores to analyze than in Gammie et al. (2003, see their Figure 4) and we are able to demonstrate a core mass distribution similar to the general IMF. Note that not all material within a core will finally end up in the protostar. In §5, we will find out that more than half of the core mass resides outside the central region. Even material inside the central region may be evaporated or ejected away during star formation. In addition, the cores that we observed in the simulation may further fragment into binaries or even multiple systems if the resolution of the simulation can be increased further. The Truelove et al. (1997) criterion is eventually violated at the centers of our cores (see §5) so we cannot follow this evolution. We note that, despite the emphasis on artificial fragmentation because of the particular test problem chosen by Truelove et al. (1997), regions that are underresolved by their criterion may also show too little fragmentation if the geometry is less pathological. Competitive coagulation (e.g Silk & Takahashi, 1979; Lejeune & Bastien, 1986; Murray & Lin, 1996; Bonnell et al., 2001a, b) may reshape the final IMF when gravitation becomes more important. Therefore, we cannot reach a definite conclusion that the final stellar IMF is similar to the core mass spectrum.

### 4.4 Core shapes and Orientation

Gammie et al. (2003) recently studied gravitationally collapsed clumps in decaying MHD turbulence simulations on 2563 grids with different magnetic field strengths. In our 5123 driven turbulence simulation, $`v_A/c_s`$ = 1.5, which is close to the value for their run C. Using the same definition and symbols of principal axis lengths and ratios as in Gammie et al. (2003), we calculated the principal axes of the cores in our simulation at $`t=0,0.65\tau _{\mathrm{ff}}`$, and 1.31 $`\tau _{\mathrm{ff}}`$ and plotted the axis ratios in Figure 13. It is useful to divide the lower triangular part of the figure into 3 regions: a prolate region of $`b/a`$ = 0 to 0.33, a triaxial region of $`b/a`$ = 0.33 to 0.66, and an oblate region of $`b/a`$ = 0.66 to 1. The result from our higher resolution, driven run basically resembles that of Gammie et al. (2003, see their Figure 12). The majority of cores are prolate or triaxial in shape, even near the end of our simulation. We also see a small number of cores that become oblate late in the simulation, with large axis ratios. This implies that some cores become disk-like. Figure 14 shows the angles between the shortest core axis and the density weighted mean magnetic field in cores at $`t=0.65\tau _{\mathrm{ff}}`$ (dot-dash line) and 1.31 $`\tau _{\mathrm{ff}}`$ (solid line). We see a slight correlation in the angles between the shortest core axis and the density weighted local mean magnetic field at $`t=0.65\tau _{\mathrm{ff}}`$ and an even stronger correlation at 1.31 $`\tau _{\mathrm{ff}}`$, which is different from the result of Gammie et al. (2003). This difference could possibly be a result of higher resolution in our simulation.

## 5 Disk-like Cores

In Figure 13, we see a relatively small number of cores inside the oblate region. The cores are too small to study at $`t`$ = 0 but large enough to study at $`t=1.31\tau _{\mathrm{ff}}`$. Most of these cores have a well defined disk-like appearance that could provide some hints about protostellar disk formation. Before we begin the discussion of disk-like core, we must emphasize that for some cores, the density of the inner cells grows so high that that the local Jeans length is no longer resolved in their centers, violating the Truelove et al. (1997) criterion. Unfortunately, these cores are also the largest whose physical properties can be studied with reasonable statistical significance. Other models using SPH and sink-particles also find disk-like cores in the later stages (e.g. Bate, Bonnell, & Bromm, 2003). Although our cores are under-resolved at $`t=1.31\tau _{\mathrm{ff}}`$, the external disks of cores may still provide some useful information on the dynamics of core accretion. With the clear realization that even higher resolution will ultimately be needed to properly address the structure of these objects, we proceed with an analysis of the data we have in hand, bearing in mind that our 5123 model is the highest resolution driven, self-gravitating MHD simulation currently available.

### 5.1 General Properties of disk-like Cores

We select 12 cores inside or near the oblate region in Figure 13 with well defined disk-like appearance and tabulate their mass and diameters in Table 1. The cores selected are shown in the surface density map (Figure 15). Core numbering is arbitrary. The orientation of the disk-like cores is uncorrelated to the 3D grid in the simulation despite even though they accidentally appear correlated in Figure 15. We have examined the orientation of the shortest axis of each disk-like core and they are basically random with respect to the 3D grid. The corresponding masses and radii using the high-density scaling mentioned in §2 are also shown for easy comparison with observations. The masses of cores range from 0.1 to 2.5 M and the diameters of the disks are a few thousand AU. The accretion rates of the cores will be discussed later in this section. The 12 cores can be roughly categorized in two groups: (1) a cluster of seven cores that remain close to each other in a region with higher mean magnetic field, and (2) five cores that are basically isolated and with higher spatial velocity (see Table 1). Cores in category 1 have higher mass, central density, magnetic pressure, and accretion rates compared to cores in category 2. The mean magnetic field is an order of magnitude higher inside the cluster. The size of both grouped and isolated cores appears similar. By looking at these two categories of disk-like cores, we can see some manifestation of mass segregation in that low-mass cores are preferentially found at large radii from the cluster whereas massive cores sink towards the center (e.g. Bonnell & David, 1998).

The mass and sizes of the 12 disk-like cores shown in Table 1, using high-density scaling, are in the range of the recent observations of young stellar objects that possess disk-like envelopes or circumstellar disks (e.g. Zhang, Hunter, & Sridharan, 1998; Fuente et al., 2001; Hogerheide, 2001; Jayawardhana et al., 2001; Wiseman et al., 2001). Most of these observations suggest that the sizes of the circumstellar disks of many young stellar objects are just a few hundred AU and the sizes of the disk-like envelopes are a few thousand AU. In addition, most of the studied protostars inside the disk-like cores are low or medium mass stars of mass $``$ 1 M. The estimated mass of the circumstellar disk is only a small fraction of the central star. The more massive objects we find in our simulation are better thought of as disk-like envelopes rather than true protostellar disks.

### 5.2 Resolution study of disk-like Cores

In order to better understand how well the cores are being resolved in our simulation, we compare them with typical cores formed in other lower resolution simulations. In Figure 16, we show images of the projected column density of four cores from four simulations at the time when the total mass inside cores is $``$ 10% of the total box mass. In the 643 run (first row of Figure 16), there are four cores found by CLUMPFIND and we cannot identify any disk-like core (axis ratio $`>`$ 2:1). In the 1283 run, there are 24 cores found, and still no disk-like core is identified (second row in Figure 16). In the 2563 run, there are 38 cores determined and four are disk-like (third row in Figure 16). In the 5123 run, there are 67 cores and twelve are disk-like. In the last row of Figure 16, we show core 1 from the 5123 run. Note that the coordinates shown in Figure 16 are grid zones, not spatial distance. The spatial size of the core in 643 is the largest but the diameters of cores in 1283, 2563, and 5123 are about the same, indicating the convergence of core size with increasing simulation resolution. We can also see that the disk-like appearance of cores becomes more pronounced with higher resolution. Figure 17 shows the surface density profiles of the cores in Figure 16. The curves are offset for ease of comparison. We can see a converging density profile of the core with almost uniform density at the center and a power law distribution of the disk/envelope with higher resolution, as predicted by other single cloud core simulations. Core 1 has the largest spatial diameter and is selected for further study below.

### 5.3 Physical Properties of a Disk-like Core

In Figure 18, we show close-up column density images of core 1 projected along the x-, y-, and z-axes. At the end of the simulation, core 1 seems to be an isolated core as it is far from other cores. The closest neighbor is about half a box size ($`>`$ 0.14 pc) away from it. In fact, core 1 passes through a dense region at $`t1\tau _{\mathrm{ff}}`$, interacts with several cores, and accretes a large amount of material. The disk only appears clearly after this; the core remains isolated thereafter. Gravitational interaction among protostars in nascent star clusters is argued to tidally truncate or even disrupt accretion disks (e.g. Clarke & Pringle, 1991; Hall, Clarke, & Pringle, 1996; Scally & Clarke, 2001; Kroupa & Burkert, 2001). The disk-like cores in this simulation suggest the alternative that disks could survive at least some interactions. Core 1 is not massive ($``$ 0.457 M) but the extent of the accretion disk is the largest ($`>`$ 5000 AU) of all the cores in the simulation. In the z-direction plot of Figure 18, we can see some evidence of spiral structure in the disk.

In Figures 19a-c, the radial profiles of energies, velocities, and surface densities of core 1 are shown. All the physical quantities in the radial direction are the average values of zones binned at the same radius from the center of mass. We can see that the magnetic energy remains relatively small compared to the kinetic energy (Figure 19a). The radial velocity (dashed line), $`v_{rad}`$, of the disk shown in Figure 19b shows a typical picture of angular momentum transfer inside a viscous disk, in which friction tries to spin up the outer part of the disk and spin down the inner part (Lynden-Bell & Pringle, 1974). However, the largest source of viscosity is likely numerical in our model, so we cannot draw any quantitative conclusions from this figure.

From the choice of sound speed in the simulation, the unit velocity corresponds to 2 km s-1. In Figure 19b, the radial infall rate is $``$ 0.07 in model units, corresponding to $``$ 0.14 km s-1, on the order of what is observed or modeled (Williams et al., 1999; Myers, Evans, & Ohashi, 2000). Error bars that plot variances of radial and rotational velocities estimated at each radius are also shown. The disk is basically in a rotationally supported state. Note that the inward or outward radial velocity of the disk is not the velocity used in estimating the accretion rate of material to the core in Table 1, because, as we will show later in this section, material mostly accretes onto the core in the direction perpendicular to the disk after the disk is well formed.

The dot-dash line in Figure 19b is the theoretical Keplerian rotation curve ($`v_{rot}r^{0.5}`$). The actual core rotational curve (solid line), $`v_{rot}`$, shows typical Keplerian rotation in most parts of the disk. The angular momentum conserved rotational curve (dotted line) is also shown for comparison. The rotational speed of the core 1 is $``$ 0.25–0.35 within the radius of 7 zones, which corresponds to $``$ 0.13–0.18 rad Kyr-1 in angular velocity. It takes about 40,000 yr for the core to complete one full rotation. Figure 19c shows the surface density profile (solid line) along with the best fitting power law of surface density $`r^{2.4}`$ (dash line).

Observationally, most of the information about true protostellar disks is obtained by using the position-velocity (PV) diagram (Richer & Padman, 1991). Note that the information on rotation and infall obtained from PV diagrams is subject to large uncertainties if the source is poorly resolved. Figure 20 shows the PV diagram of core 1, observed in the x-direction along the y-axis (see Figure 19a). The contours in the PV diagram are logarithmic in column density, which is proportional to the observed intensity if the disk is optically thin and isothermal. In Figure 20a, the PV diagram is shown at full resolution. The relative position ($`\mathrm{\Delta }`$Y) in zones is plotted against the LSR velocity in model units. The shifted appearance of velocity in the PV diagram is the result of disk rotation. For qualitative comparison to observations, we convolve the map with a circular Gaussian beam having FWHM of three grid zones (Figure 20b). The dashed curve shows the distribution of the rotation velocity proportional to $`r^1`$, which is characteristic of angular momentum conserving rotation. The solid curve shows the distribution of the rotation velocity proportional to $`r^{0.5}`$, which is characteristic of Keplerian rotation. We can see that Keplerian rotation better describes most parts of the disk in this core, although the poor resolution of the PV diagram limits the strength of this conclusion. If we use this PV diagram to estimate the point mass of the central star, we would conclude that the central mass of core 1 is about 0.21 M, using the high-density scaling, which is the mass enclosed by the region of radius $`<`$ 3.4 zones at the center of this core. We find that the mass of the disk outside this central region is more than 53% of the total core mass (Table 1). An example of an object with properties similar to core 1, in high-density scaling, is IRAS 05413-0104 (Wiseman et al., 2001). IRAS 05413-0104 is a flattened disk-like core of average volume density 105 cm-3, with an envelope of mass $``$ 0.2 M, and a disk with diameter of $``$ 12,000 AU, roughly twice the size of core 1.

In Figure 21, we show the velocity and magnetic field components perpendicular and parallel to the disk for core 1. The velocity field on the disk plane clearly shows rotation (Figure 21b). The edge-on view of the velocity (Figure 21a) shows that material above and below the disk is mostly falling onto the core along the magnetic field lines, as expected when magnetic field is present.

The accretion rates of all cores at $`t=1.31\tau _{\mathrm{ff}}`$ are listed in Table 1. The accretion rate is calculated by integrating the material falling onto the whole area of the disk above and below. From Table 1, the accretion rates of the 12 cores are in the range of $`10^4`$$`10^5`$ M yr-1. If we assume the accretion rate is constant, the ages of the cores are about 15,000 to 50,000 years old, which is only $`0.14`$–0.46 $`\tau _{\mathrm{ff}}`$. Therefore, the accretion rates of cores must have dramatically increased near $`t=1.31\tau _{\mathrm{ff}}`$, as expected from a core collapse scenario. The rate of accretion onto core 1 is $``$ 9.72 $`\times 10^6`$ M yr-1. On average, the cores inside the cluster have accretion rate more than twice of those outside. The accretion rates of all 12 cores are basically at the high-end of the range of observed envelope infall rates for protostellar cores (Ostriker, 1998).

In §4.4, we noted a weak correlation between the core minor axis and the local magnetic field direction. Here in Figure 21c, the vertical component of magnetic field lines does show a bipolar appearance in core 1, but the magnetic field is bent, due to the movement of the core relative to the ambient material. The magnetic field component along the disk also shows rotating features, which indicates that the magnetic field is frozen and dragging along with the rotating disk (Figure 21d). This magnetic structure is common for all the disk-like cores. In Figure 22, we show the 3D structure of the magnetic field in core 6, which is helical as a result of core rotation. The angular momenta of cores are expected to be transferred to the ambient material by such helical Alfvén waves, which in return brake the rotation of the cores (e.g. Mouschovias & Paleologou, 1979, 1980; Basu & Mouschovias, 1994).

## 6 Conclusion and Discussion

Molecular clouds are magnetized and in a state of turbulent motion. It remains an open question whether the turbulence is sub-, trans-, or super-Alfvénic (e.g. Mac Low & Klessen, 2004). Incompressible hydrodynamic turbulence has been thoroughly studied, but compressible MHD turbulence is still a relatively young area, especially as applied to astrophysical problems. In the last decade, numerical studies of compressible turbulence in molecular clouds have revealed important results that may directly relate to the early stages of star formation. Such simulations still suffer from either lack of magnetic field, lack of self-gravity, or insufficient resolution. In this paper, we present the results of our latest studies on self-gravitating core formation in a super-Alfvénic turbulent molecular cloud.

A statistical analysis of the simulation after the turbulence has fully developed shows that the power spectra of velocity and magnetic field basically follow Kolmogorov scaling, even though the system is compressible and magnetized. This result is consistent with some other recent high resolution studies on super-Alfvénic turbulence (e.g. Müller & Biskamp, 2000; Cho, Lazarian, & Vishniac, 2003). The power spectra for density, kinetic, and potential energies show shallower power-law indexes.

We calculated the core mass spectrum in our high-resolution simulation including both MHD and self-gravity, and followed its evolution. At least at early times when the model is completely resolved, our core mass spectrum is consistent with two main predictions of the turbulent fragmentation theory advanced by Padoan & Nordlund (2002): that the slope of the velocity spectrum defines the slope of the high-mass wing of the core mass spectrum, and that the low mass end flattens and turns over. Padoan & Nordlund (2002) propose that this could account for the form of the stellar IMF. At later times in our model, the slope of the high-mass wing becomes shallower, possibly because of core coalescence or lack of resolution of fragmentation in central regions of cores. Our core encounter rate may well be overestimated because of the choice of periodic boundary conditions and grid resolution.

The simulation that we have presented in this paper is based on initial conditions that are both favorable and hostile to the formation of cores. In the case simulated, the molecular gas is magnetically supercritical with $`M/M_{\mathrm{cr}}`$ = 8.3. We expect a priori that gravitational collapse is unavoidable and will occur eventually. At the same time, the system is subjected to a strong turbulent driving force producing flows with rms thermal Mach number $`M_{\mathrm{rms}}`$ = 10. The turbulence is strong enough to formally support the entire region, so it could prevent or delay the collapse of cores or even destroy some cores in the process of collapsing. In Paper II, it was shown that in a simulation on a 2563 grid, with the same initial conditions, the magnetic and turbulent support could not prevent local collapse of cores in the molecular cloud. Cores form almost immediately after the gravitational force is turned on. One of the motivations for our higher resolution simulation—twice the spatial resolution and eight times the mass resolution of the 2563 model of Paper II—is to determine whether stronger Alfvén wave support may occur when shorter wavelengths are resolved in a simulation using the same initial conditions. We find that gravitational collapse still wins over magnetic and turbulent support. The 5123 simulation lasts about 2/3 of a full free-fall time, after the gravitation is turned on, before the first cell in the computation violates the Jeans criterion Truelove et al. (1997) that the local Jeans length be resolved, a factor of four longer than the 2563 simulation. A resolution study shows convergent trends with increasing simulation resolution of both the fraction of the box mass in bound cores (Fig. 5) and the core size (Fig. 17).

Our simulation began with the entire box magnetically supercritical. However, it is often argued that gravitational fragmentation will produce smaller and smaller regions that eventually become subcritical (Mestel & Spitzer, 1956). We found that gravitationally bound cores formed and collapsed until they were rotationally supported (Fig. 9), while still having masses at least an order of magnitude over the critical mass for magnetic support (Figs. 8a-b). We further found that the central magnetic field strength depends on the central density as $`B_c\rho _c^{1/2}`$, suggesting that observations of such relations do not necessarily point to subcritical cores.

Most of the cores appear to be prolate or triaxial in shape, consistent with simulation results from Gammie et al. (2003). However the correlation of the shortest axis to the density-weighted local magnetic field of the cores is found to be stronger in our simulation.

In our simulation, we observe interaction between cores during close encounters. The continuing accretion of large amounts of material, and merging of cores results in cores with high specific angular momentum that helps the formation of disk-like structure. The final outcome is highly unpredictable as both magnetic and turbulent forces are present during the gravitational collapse process. In our simulation, we have identified twelve disk-like cores with axis ratio $`>`$ 2:1 and large enough for statistical analysis. The sizes and mass of the twelve cores are in the range of young protostellar objects observed with flattened disk-like structure (Table 1). All the disk-like cores are rotationally supported. Magnetic pressure support is relatively unimportant for these cores, perhaps due to the large $`M/M_{\mathrm{cr}}`$ ratio that we chose as an initial condition. The surface density of most of the disk-like cores has a power law distribution, as observed in real cores, and as predicted in isolated core collapse simulations. Instead of being tidally truncated or disrupted, the core disks survive and flourish among strong interactions.

Most of the disk-like cores appear to be in Keplerian rotation from their PV diagram. Observationally, many protostellar disks are found to be in Keplerian rotation (e.g. Wiseman et al., 2001, and references therein). There are still some observed protostellar disks found to be in non-Keplerian rotation (e.g. Myers, Evans, & Ohashi, 2000). Wiseman et al. (2001) points out that the non-Keplerian rotation curve reflects a mass distribution that is not yet dominated by a central condensation. Also, a system could be in Keplerian rotation but a high optical depth would result in velocity maps that are dominated by the motion in the outer layers along the line of sight toward the center of the core. For the twelve rotationally supported disk-like cores, a large amount of matter is still accreting perpendicular to the disks (Figure 21a); converted to astrophysical units this corresponds to substantial accretion rates (see Table 1). When the mass of the disk keeps increasing, the gravitational stability of the rotational disk may be affected. Unfortunately, because of the limited resolution, we cannot follow the core collapse long enough to see consequences beyond the first hint of spiral structure in the disk.

There is as yet no direct detection of disks around high-mass ($`>`$ 1 M) protostars. Fuente et al. (2001) suggests that for stellar mass $`>`$ 5 M, destruction mechanisms like radiation pressure on dust or ionization by stellar ultraviolet emission may be responsible for the rapid erosion of the outer disk. However, these mechanisms would only operate at the later stages of protostellar evolution. Zhang, Hunter, & Sridharan (1998) point out that the lack of direct detections of disks around high-mass protostars may simply be due to the greater distances to high mass star formation regions, as well as their clustering, which further complicates the kinematics. Therefore, the existence of disks around high-mass protostars at the very early stage of core collapse remains uncertain. It is important to note that in a number of objects studied at high resolution with radio and millimeter-wave interferometry, flattened structures and Keplerian motions argue strongly for the existence of rotationally supported disks (e.g. Sargent & Beckwith, 1991; Koerner, Sargent, & Beckwith, 1993; Dutrey, 1996). Improvements in the resolution of both observations and simulations on the early stage of core collapse in molecular clouds are vital to understand the details of star formation inside magnetized turbulent molecular clouds.

M.-M. M. L. acknowledges partial support by NASA Astrophysical Theory Program grant NAG5-10103 and NSF CAREER grant AST 99-85392. F. Heitsch is supported by a Feodor-Lynen fellowship of the Alexander-von-Humboldt Foundation and by NSF grant AST-0098701. We thank R. S. Klessen, P. Padoan, and E.G. Zweibel for useful discussions. Computations were performed at the NCSA using the 256-processor SGI/Cray Origin 2000 through NRAC allocation QUS.