A Simple Law of Star Formation
We show that supersonic MHD turbulence yields a star formation rate (SFR) as low as observed in molecular clouds (MCs), for characteristic values of the free-fall time divided by the dynamical time, , the alfvénic Mach number, , and the sonic Mach number, . Using a very large set of deep adaptive-mesh-refinement simulations, we quantify the dependence of the SFR per free-fall time, , on the above parameters. Our main results are: i) decreases exponentially with increasing , but is insensitive to changes in , for constant values of and . ii) Decreasing values of (stronger magnetic fields) reduce , but only to a point, beyond which increases with a further decrease of . iii) For values of characteristic of star-forming regions, varies with by less than a factor of two. We propose a simple star-formation law, based on the empirical fit to the minimum , and depending only on : . Because it only depends on the mean gas density and rms velocity, this law is straightforward to implement in simulations and analytical models of galaxy formation and evolution.
Subject headings:ISM: kinematics and dynamics — stars: formation — magnetohydrodynamics — turbulence
Stars are the result of the gravitational collapse of cold gas, but on large scales the conversion of cold gas into stars takes much longer than a collapse time. Early explanations for this low star formation rate (SFR) involved the magnetic support against the collapse of giant molecular clouds (GMCs), while the more recent scenario of turbulent fragmentation relies on the support of GMCs by supersonic turbulence (Mac Low & Klessen, 2004). However, numerical simulations of star formation in turbulent clouds have generally failed to yield SFRs as low as observed values111After submitting this Letter, we became aware of a closely related paper by Federrath and Klessen, submitted to the Astrophysical Journal (arXiv:1209.2856), presenting a large number of star formation simulations that are complementary to ours. (see § 4).
This work shows that supersonic MHD turbulence results in a SFR as low as observed in GMCs, for characteristic values of the free-fall time divided by the dynamical time, , the alfvénic Mach number, , and the sonic Mach number, (defined in § 2). Using a very large set of adaptive-mesh-refinement (AMR) simulations of driven MHD turbulence, including self-gravity and sink particles, we quantify the dependence of the SFR per free-fall time, , on , confirming the results of our previous uniform grid simulations in the case of a very weak mean magnetic field.
The SFR per free-fall time is the efficiency factor of a theoretical Schmidt-Kennicut law,
often used in models and simulations of star formation in galaxies. We find that decreases exponentially with increasing , has no dependence on alone, and varies by less than a factor of two with variations in within reasonable values for star-forming regions. We therefore propose a new empirical law of star formation that depends only on .
2. Numerical Simulations
|Runs for study of numerical convergence:|
|Runs for study of variations due to initial conditions:|
The simulations are carried out with the AMR code Ramses (Teyssier, 2002), modified to include random turbulence driving and sink particles, and optimized for supercomputers with large multi-core nodes through an OpenMP/MPI hybrid layout. As in Padoan & Nordlund (2002, 2004, 2011), we adopt periodic boundary conditions, isothermal equation of state, and solenoidal random forcing in Fourier space at wavenumbers ( corresponds to the computational box size). We choose a solenoidal force to make sure that dense, star-forming structures are naturally generated in the turbulent flow, rather than directly imposed by the external force.
The simulations are first run without self-gravity for approximately five dynamical times, starting with uniform density and magnetic fields, and random velocity field with power only at wavenumbers . The driving force keeps the rms sonic Mach number, , at the approximate values of either 10 or 20, characteristic of GMCs on scales of order 10 pc ( is the three-dimensional rms velocity, and is the isothermal speed of sound). For the simulations with , we have chosen two different values of the mean magnetic field, corresponding to rms alfvénic Mach numbers of and 5, where , and is the Alfvén velocity corresponding to the mean magnetic field and the mean density. For the larger Mach number, , we have adopted three different values of the mean magnetic field, corresponding to , 5, and 1.25.
In the supersonic regime, the turbulent dynamo is suppressed (Haugen et al., 2004), and the magnetic energy growth saturates at approximately 2% of the kinetic energy of the turbulence, when the mean magnetic field is negligible and does not characterize the system (Federrath et al., 2011). On the other hand, if the mean magnetic field is not negligible, as in our simulations and in real GMCs, it is readily amplified, well beyond the dynamo saturation level, by stretching and compression in the turbulent flow, and thus its saturated rms value depends directly on its mean value, as found in simulations of Kritsuk et al. (2009), and in equation (20) of Padoan & Nordlund (2011). Therefore, the value of , defined with the mean magnetic field, characterizes the saturated regime and can be taken as a fundamental parameter of the system.
The final snapshots of these five turbulence simulations are used as initial conditions for the star-formation simulations with self-gravity and sink particles. Each initial condition is used for 5 to 7 simulations with different numerical values of the gravitational constant, corresponding to different values of , resulting in 28 star-formation simulations (plus 6 larger runs used to test numerical convergence, and 11 smaller runs to test both convergence and variations of with initial conditions – a total of 45 star-formation simulations). The dynamical time is defined as , where is the size of the computational domain, and , where is the mean density in the computational volume. The simulations were run until the star formation efficiency value of SFE 0.2-0.3, where is the mass in the sink particles, and the total mass (gas plus sinks), except the simulations with very low SFR, that were run until SFE 0.05-0.1.
The five initial turbulence simulations have a root grid of 128 computational cells and 5 AMR levels, thus a maximum resolution equivalent to 4,096 cells. In all runs, the first AMR level is created if the density is 2.5 times larger than the mean. Subsequent AMR levels are created when the gas density increases by a factor of 4 relative to the previous level. In contrast to Kritsuk et al. (2006) and Schmidt et al. (2009), we do not refine based on velocity or pressure gradients. However, we verified that the tails of the density PDFs are sufficiently converged in the turbulence simulations used as initial conditions (details will be given elsewhere). When self-gravity is included, the root grid is reduced to 64 cells, while the number of AMR levels is increased from 5 to 8, to achieve a maximum resolution equivalent to 16,384 cells.
Sink particles are created in cells where the gas density is larger than times the mean density (much larger than the maximum density reached in the initialization runs without self-gravity), if the following additional conditions are met at the cell location: i) The gravitational potential has a local minimum value, ii) the three-dimensional velocity divergence is negative, and iii) no other previously created sink particle is present within an exclusion radius, ( in these simulations, where is the highest-resolution cell size). These conditions are similar to those in Federrath et al. (2010). We have verified that they avoid the creation of spurious sink particles in regions where the gas is not collapsing.
The numerical parameters characterizing the simulations are given in Table 1. The simulation names contain the rms sonic Mach number, and the numerical values of the mean magnetic field and of the gravitational constant. In Table 1, we give three non-dimensional parameters, , , and , that fully characterize the simulations. The ratio expresses the relative strength of gravity and turbulence, like the virial parameter, , used in previous works. Although formally equivalent to , for example through the usual expressions valid for a uniform sphere, we prefer to use as the third parameter, in line with the choice of the Mach numbers (ratios of crossing times as well) as the first two parameters.
In Table 1, we also give the ratio between the minimum Jeans length at each level, (corresponding to the highest density before further refinement), and the computational cell size, , at that level. Due to the refinement by a factor of two in size for every increase by a factor of four in density, the minimum number of cells per Jeans length, , is the same for every level in a given simulation. The last column of Table 1 shows the values of obtained from the simulations.
At the bottom of Table 1, we also list the simulations used to test numerical convergence and to study the variations of with initial conditions. The names of these extra runs begin with their root grid size, either or computational cells. Projected density and sink particles from one of these runs, 128M10B2G50, are shown in Figure 1, at a time when SFE.
Figure 2 shows the time dependence of the SFE for the two sets of simulations with and and 1.25. In some runs with high , the SFR tends to grow with time. To extract a representative value, we measure as the slope of the least-squares fit of the SFE versus , in the interval , as in Padoan & Nordlund (2011). The simulations with the lowest values of were not integrated up to SFE=0.2. However, they do not exhibit a systematic growth of the SFR with time, and they are run for many free-fall times (up to 9 in some cases), so their average is nonetheless well defined.
The values of from all simulations are plotted in Figure 3, as a function of . A striking result, stressed by the linear-log scale of Figure 3, is the approximately exponential dependence of on (for fixed values of and ), showing that supersonic MHD turbulence is very effective at slowing down the star formation process, which seems to be controlled primarily by the ratio of the characteristic times of gravity and turbulence.
Another important result is the lack of a dependence of on , for fixed values of and . Furthermore, the dependence on is not very strong, considering the very large range of values of the simulations. This is in part due to the complex dependence of on . Down to , decreases with decreasing (increasing magnetic field strength), while for further decreasing , starts to increase.
Figure 3 shows that a simple exponential law, where is only controlled by the ratio , could describe all simulation results within a factor of two. For even weaker magnetic fields than considered here, would certainly be larger, as it should approach the values found in simulations without magnetic fields (Padoan & Nordlund, 2011). However, observations in star forming regions are consistent with (Padoan & Nordlund, 1999; Lunttila et al., 2008, 2009; Heyer & Brunt, 2012), so it is likely that approaches its minimum value, which is approximately fit by a simple exponential law, , shown as a dashed line in Figure 3. We therefore propose a new empirical law of star formation, based on this fit to the minimum :
where accounts for mass loss through jets and winds, during the formation of a star. This law only depends on the mean gas density and the rms velocity of a star-forming region, so it is easily implemented in analytical models and simulations of galaxy formation and evolution.
Equation (2) is valid as long as the fraction of the system mass, , available to form stars is larger than the characteristic stellar mass, , which we can express as . If predicted by equation (2) is too small to satisfy this condition, then , simply because the total system mass is too small. Adopting the characteristic stellar mass predicted by our turbulent fragmentation model (Padoan & Nordlund, 2002; Padoan et al., 2007), , which corresponds to approximately three times the mean Bonnor-Ebert mass, , evaluated for an external pressure equal to the ram pressure of the turbulence, and using the minimum derived here, we find the condition can be expressed as .
As shown by Table 1, the weakest gravity we have simulated corresponds to , giving and 38.6 for the runs with and 20 respectively. In good agreement with these values, in each of the runs with , we get only a single star, so the SFR is practically zero, while in the three runs with we get approximately 10-20 stars per free-fall time.
To test the numerical convergence of , we have rerun three of the simulations with and , using both twice smaller and twice larger root grid size ( and computational cells), and the same number of AMR levels, corresponding to maximum resolutions equivalent to 8,192 cells and 32,768 cells. These three larger runs were then rerun again using a 10 times larger value of the threshold density above which sink particles are created (10 times the mean density, instead of ). As shown in Table 1, the values of in the runs with and 1.09 are practically independent of numerical resolution, showing that the correct SFR can be derived even when Truelove’s condition (Truelove et al., 1997) is violated, and that the value of times the mean density (adopted in the main parameter study), above which sink particles are allowed to form, is sufficiently high to avoid dense regions that are not collapsing.
In the case of , is approximately 20% lower in the -root-grid simulations than in the corresponding -root-grid run. However, this low SFR is obtained from an average over many dynamical times. During the first dynamical time, the SFR is virtually independent of resolution. The corresponding simulation with root-grid size has the same value of as the largest simulations, showing that the difference from the -root-grid run is not due to lack of convergence, but more likely to random differences between runs after a long time, as expected in chaotic systems. The values of from the -root-grid simulations can be seen in Figure 3.
To investigate the effect of the initial conditions, we have carried out 10 simulations with twice smaller root grid size ( computational cells), and the same number of AMR levels, five of them with , and the other five with , starting from different initial conditions of fully-developed turbulence, with and . As shown in Table 1, the maximum difference in from such small samples is approximately a factor of two. The mean values of , plus and minus the rms values, are plotted as error bars in Figure 3.
The values of predicted by equation (2) are consistent with those found in GMCs. As discussed in Padoan & Nordlund (2011), accounting for the full uncertainty in the lifetime of the Class II phase ( Myr), the results of Evans et al. (2009) implies -0.12. Using the GMC sample by Heyer et al. (2009), we derive (1-), resulting in -0.12 based on our equation (2), with , in nice agreement with the values derived by Evans et al. (2009) (or and -0.19 if we double the cloud masses). Similar values are found also for the Orion Nebula Cluster, -0.09, the object with the smallest error bars in Krumholz & Tan (2007). Based on the values of in the sample of Murray (2011) we would predict an average value of , while, based on the free-free emission, Murray finds an average value of .
The AMR simulations with and presented here are consistent with our previous uniform grid simulations (Padoan & Nordlund, 2011) with and (the values of given in the plot of Figure 2 are derived from based on the case of a uniform sphere, while those reported in Figure 7 of Padoan & Nordlund (2011) were based on the total mass in the computational box, and should be multiplied by a factor of 2 for the comparison).
The SFR model of Padoan & Nordlund (2011) adopted the simplifying assumption that, in the MHD case, the effective postshock ratio of gas to magnetic pressure, , and the density pdf do not depend on the mean magnetic field strength, for reasonable values of . The complex, but weak dependence of on found in the simulations justifies that approximation. However, by neglecting the dependence of on both and , the model yields a direct dependence on , not found in our AMR simulations. Equation (20) of Padoan & Nordlund (2011), relating to and , the ratio of gas to magnetic pressure computed with the average magnetic field strength, , was derived under the assumption of flux freezing, neglecting the postshock thermal pressure and dynamical alignment of magnetic and velocity fields. Padoan & Nordlund (2011) noticed that when is averaged at high densities, it becomes independent of , so they took , which also provided a good fit to the density pdfs. However, they did not test the dependence of on .
It may be more physical to retain the dependence of on even at high densities, because the postshock magnetic pressure tends to balance the ram pressure of the turbulence causing the shocks. Assuming instead of , we obtain practically independent of , as found in our AMR simulations. This assumption yields a density variance in MHD given by , for , while the approximation in Padoan & Nordlund (2011) gave (this issue will be addressed in a future study).
Although this correction eliminates the dependence on , the model does not predict the subtle dependence of on , as it does not account for the increasing anisotropy and dynamical alignment of magnetic and velocity fields, with decreasing . Therefore, the simple empirical law derived above is a better choice for practical applications.
The models by Krumholz & McKee (2005) and Hennebelle & Chabrier (2011) do not include the effect of magnetic fields. As shown by our previous uniform-grid simulations (Padoan & Nordlund, 2011), even with a very weak mean magnetic field, corresponding to , is approximately three times smaller than in the non-magnetized case. The AMR simulations presented here show that, with more realistic, larger magnetic field values, decreases even further. Analytical models of star formation neglecting magnetic fields should not yield low values of consistent with observations of star-forming regions. If they do, they are inconsistent with the numerical experiments.
The new AMR simulations of star formation presented in this work show that i) decreases exponentially with increasing , but is insensitive to changes in (in the range ), for constant values of and . ii) Decreasing values of (increasing magnetic field strength) reduce , but only to a point, beyond which increases with a further decrease of . iii) For values of characteristic of star-forming regions, varies with by less than a factor of two. We propose a simple law of star formation depending only on , based on the empirical fit to the minimum : .
This law shows that MHD turbulence is very effective at slowing down the star-formation process and can explain the low average SFR in molecular clouds. Because it only depends on the mean gas density and rms velocity of a star-forming region, the star-formation law we propose is straightforward to implement in simulations and analytical models of galaxy formation and evolution. Future work should test its effect in simulations where the feedback of SNs is accounted for, but the process of star formation is not spatially resolved and needs to be modeled.
- Evans et al. (2009) Evans, II, N. J., Dunham, M. M., Jørgensen, J. K., Enoch, M. L., Merín, B., van Dishoeck, E. F., Alcalá, J. M., Myers, P. C., Stapelfeldt, K. R., Huard, T. L., Allen, L. E., Harvey, P. M., van Kempen, T., Blake, G. A., Koerner, D. W., Mundy, L. G., Padgett, D. L., & Sargent, A. I. 2009, ApJS, 181, 321
- Federrath et al. (2010) Federrath, C., Banerjee, R., Clark, P. C., & Klessen, R. S. 2010, ApJ, 713, 269
- Federrath et al. (2011) Federrath, C., Chabrier, G., Schober, J., Banerjee, R., Klessen, R. S., & Schleicher, D. R. G. 2011, Physical Review Letters, 107, 114504
- Haugen et al. (2004) Haugen, N. E. L., Brandenburg, A., & Mee, A. J. 2004, MNRAS, 353, 947
- Hennebelle & Chabrier (2011) Hennebelle, P. & Chabrier, G. 2011, ApJ, 743, L29
- Heyer et al. (2009) Heyer, M., Krawczyk, C., Duval, J., & Jackson, J. M. 2009, ApJ, 699, 1092
- Heyer & Brunt (2012) Heyer, M. H. & Brunt, C. M. 2012, MNRAS, 420, 1562
- Kritsuk et al. (2006) Kritsuk, A. G., Norman, M. L., & Padoan, P. 2006, ApJ, 638, L25
- Kritsuk et al. (2009) Kritsuk, A. G., Ustyugov, S. D., Norman, M. L., & Padoan, P. 2009, in Astronomical Society of the Pacific Conference Series, Vol. 406, Numerical Modeling of Space Plasma Flows: ASTRONUM-2008, ed. N. V. Pogorelov, E. Audit, P. Colella, & G. P. Zank, 15
- Krumholz & McKee (2005) Krumholz, M. R. & McKee, C. F. 2005, ApJ, 630, 250
- Krumholz & Tan (2007) Krumholz, M. R. & Tan, J. C. 2007, ApJ, 654, 304
- Lunttila et al. (2008) Lunttila, T., Padoan, P., Juvela, M., & Nordlund, Å. 2008, ApJ, 686, L91
- Lunttila et al. (2009) —. 2009, ApJ, 702, L37
- Mac Low & Klessen (2004) Mac Low, M.-M., & Klessen, R. S. 2004, Reviews of Modern Physics, 76, 125
- Murray (2011) Murray, N. 2011, ApJ, 729, 133
- Padoan & Nordlund (1999) Padoan, P. & Nordlund, Å. 1999, ApJ, 526, 279
- Padoan & Nordlund (2002) Padoan, P. & Nordlund, Å. 2002, ApJ, 576, 870
- Padoan & Nordlund (2004) —. 2004, ApJ, 617, 559
- Padoan & Nordlund (2011) —. 2011, ApJ, 730, 40
- Padoan et al. (2007) Padoan, P., Nordlund, Å., Kritsuk, A. G., Norman, M. L., & Li, P. S. 2007, ApJ, 661, 972
- Schmidt et al. (2009) Schmidt, W., Federrath, C., Hupp, M., Kern, S., & Niemeyer, J. C. 2009, A&A, 494, 127
- Teyssier (2002) Teyssier, R. 2002, A&A, 385, 337
- Truelove et al. (1997) Truelove, J. K., Klein, R. I., McKee, C. F., Holliman, J. H., Howell, L. H., & Greenough, J. A. 1997, ApJ, 489, L179