Summary by Aaron Bray and Gongjie Li
The kinematics of newly formed star clusters are interesting both as a probe of the state of the gas clouds from which the stars form, and because they influence planet formation, stellar mass segregation, cluster disruption, and other processes controlled in part by dynamical interactions in young clusters. However, to date there have been no attempts to use simulations of star cluster formation to investigate how the kinematics of young stars change in response to variations in the properties of their parent molecular clouds. In this Letter, we report the results of turbulent self-gravitating simulations of cluster formation in which we consider both clouds in virial balance and those undergoing global collapse. We find that stars in these simulations generally have velocity dispersions smaller than that of the gas by a factor of ∼5, independent of the dynamical state of the parent cloud, so that subvirial stellar velocity dispersions arise naturally even in virialized molecular clouds. The simulated clusters also show large-scale stellar velocity gradients of ∼0.2–2 km s−1 pc−1 and strong correlations between the centroid velocities of stars and gas, both of which are observed in young clusters. We conclude that star clusters should display subvirial velocity dispersions, large-scale velocity gradients, and strong gas–star velocity correlations regardless of whether their parent clouds are in virial balance, and, conversely, that observations of these features cannot be used to infer the dynamical state of the parent gas clouds.
Since most stars form in clusters with greater than 100 members, the kinematics of these dense environments is of substantial interest not only in understanding cluster formation, but also general properties of planet formation, the stellar binary fraction, mass segregation, and accretion of clusters onto larger systems, to name just a few. However, dynamical measurements of these incipient clusters can be difficult because they remain obscured by gas during the early stages of their evolution. Luckily, however, one system with thousands of visible stars — the Orion Nebula Cluster (ONC) — has been well-studied and provides a suitable comparison for simulations of clusters’ kinematic properties.
Fürész et al. (2008) and Tobin et al. (2009) both offer examples of the kinematic studies that are possible with the ONC. Those papers extracted spectra of 1215 and 1613 stars, respectively.
Shown above, in figures 10 and 3, respectively, they both find substantial north-south velocity gradients in the ONC, and also that the stellar velocity distribution was well-matched to the 13CO tracer gas. To the extent that the gas and star velocity distributions do not trace each other, Tobin et al. argues that binary motion and the inclusion of a few ONC non-members are likely responsible. In conclusion, Fürész et al. and Tobin et al. suggest that these observations most likely indicate that the ONC is a subvirial gaseous filament undergoing global collapse.
Furthermore, recent N-body simulations have shown that subvirial stellar populations models can reproduce observations of the ONC. Allison et al. (2009) show that a clumpy, subvirial cluster undergoes rapid mass segregation — fast enough that the ONC’s mass distribution need not be primordial, but could rather be of dynamic origin. Proszcow et al. (2009), meanwhile, considers subvirial, aspherical systems and conclude that the kinematic observations of the ONC, such as the velocity gradient can be substantially reproduced given the correct viewing angle.
However, other work has shown that not all components of GMCs share the same kinematics, so Offner et al. propose that is far from clear that subvirial stellar populations are indicative of subvirial gas. Moreover, the other recent simulations have been purely N-body, and this particular system clearly would benefit from a hydrodynamic treatment.
These simulations use a gravitohydrodynamic (GHD), adaptive mesh refinement (AMR) code (c.f., Truelove et al. 1998). This code takes into account both the self-gravity of the system (gas + stars) and the hydrodynamics of the gas. A hyperbolic solver is used to locally solve the three hydrodynamic equations (i.e., conservation of mass, momentum, and energy; Eq. 1-3) sequentially in each dimension using update operators. That is,
Qn+2= LxLyLz[LzLyLx (Qn)]
Meanwhile, an elliptic solver is used to iterate over the entire numerical grid to solve Poisson’s gravitational equation (Eq. 4). The density at a given node is calculated using the weighted mean of neighboring cells, while the Laplacian at a given node is calculated based on neighboring nodes, as shown in the figure at right.
An AMR code works by creating and eliminating grid cells on the fly, in order to provide better resolution in the densest environments while not wasting computational power in sparse environments. This particular AMR code works in roughly three steps. First, once all densities and other scalars has been calculated, it flags all cells in need of refinement. Second, it analyzes the complete set of flagged cells to determine where, in fact, refined grids will be inserted (or removed). Two factors need to be balanced. One the one hand, a perfectly efficient code (grid efficiency = 1.0) would insert a new grid only in those new cells that needed it. However, these is a significant overhead cost to creating and maintaining each refined grid, so having many, small, disconnected refinements is actually a waste of computational power, as well. Thus, it is often preferable to refine large regions. This code is optimized with a so-called grid efficiency of ~ 0.7. Lastly, a buffer of refined cells is placed around all refined regions, so that no two cells have more than one level of refinement difference between them.
This suite of simulations has three different runs — demarcated U1, D1, and D2 — that differ based on the kinematics of the gas. As summarized in Table 1 (below), all three runs are initialized with driven turbulent motion (with gravity turned off) in order to obtain a steady-state in which the virial parameter α = 5 σ2R/GM = 1.67. In simulation U1, energy injection is turned off at t = 0, mimicking the global collapse suggested by Tobin et al. (2009). In the D1 simulation, energy injection continues (although declining), maintaing a constant Mach number and virial parameter.
This corresponds to the physical case where stell ar feedback prevents collapse. Finally, in simulation D2, energy is injected at a constant rate, so that the virial parameter increases over time. These simulations are all run for the mean free-fall timescale (tff) of the GMC, and star formation is modeled by a sink particle that replaces the gas whenever the Jeans criterion is reached on the finest refinement level of AMR code. Table 1 shows both the initial conditions, as well as the final Mach number and virial parameters for each of the simulations.
3.1. Stellar Velocity Dispersion
Figure 1 shows the stellar velocity dispersion as a function of time, for each of the three simulations. The top plot is normalized to the sound speed in the gas, and the bottom plot is normalized to the velocity dispersion of the gas.
Other than the fact that D2 takes longer to begin star formation due to the higher turbulence, the overall velocity dispersions are remarkably similar. The final dispersions all fall within 0.1 < σ*/σgas < 0.3; that is, the stars are substantially kinematically colder than the gas in all simulation runs. These dispersions correspond to α ~ 0.1 for the entire system (stars + gas) in all three simulations, and to α ≤ 1.1 for only stars (assuming all the gas disperses). Thus, only in the case when the gas virial parameter is significantly supervirial (α = 3.4) does the stellar virial parameter reach α = 1, and then, only barely.
3.2. Projected Stellar Velocities
The author used the final time of the D2 simulation, which corresponds to the final virial parameter being 3.4 for the most detailed examination on the projected stellar velocities. The stellar velocities are projected along the three cardinal sight lines (x, y, z). The authors binned nine closely located stars as a group and error bars indicating the velocity dispersion in each bin. The diamonds in Figure 2 shows the projected stellar velocity for each star. This plot (Figure 2) looks similar to the plots for ONC in Proszkow et al. 2009.
Fitting the line-of-sight velocities with the following function,
V = V0 +∆VxX + ∆VyY; |∆V| = Sqrt (∆V^2*X + ∆V ^2*Y)
the authors calculate the velocity gradients. The authors estimate the velocity gradients in different scales by considering the velocity gradients in boxes with different sizes (full box; half boxes, and quarter boxes), as shown in Figure 3. The x axis corresponds to the size of the boxes and the vertically aligned dots correspond to boxes with the same size. For a fixed box size, the authors compared the velocity gradients at three different times — 1.0tff, 0.9tff, and 0.8tff — with the ONC data. The magnitude of the gradients increases in both the simulations and the ONC for smaller boxes due to smaller statistics, where stars separated by several shock fronts may easily acquire significantly different velocities. Proszkow et al. (2009) find similar radial velocity gradients to the simulation for the ONC. The authors noted that the velocity gradients are not necessarily a hallmark of subvirial turbulence or a strongly collapsing cloud. Instead, large-scale velocity gradients appear to be a signature of the turbulent power in large-wavelength modes.
3.3. Gas Velocities
The line-of-sight 13CO gas and star velocities are shown in Figure 4. The intensity is summed over all y values for x and v. Similar to that seen the star and gas velocity maps, the simulated image shows a strong correlation between gas and star velocities. Averaging and excitation of 13CO serve to improve the correlation. They improve the correlation because summing over y-direction, the velocity is density weighted, and because excitation of 13CO only in relatively dense gas that is closely associated with forming stars. However, the low density gas carries the bulk of the kinetic energy due to the anti-correlation between velocity and density. The subthermal excited lower density gas that dominates the dispersion.
Some notes on the procedure from Offner et al. (2008): The procedure of modeling the 13CO is included in Offner et al. (2008), where the authors generate a position–position–velocity (PPV) data cube from simulation and compute the emissivity as a function of density assuming that the gas is in statistical equilibrium, that it is optically thin, and that radiative pumping by line photons is negligible.
- The simulations represent the initial conditions of stars in clusters prior to dynamic evolution and cloud dispersal.
- Subvirial stellar velocity dispersions arise naturally from clouds in virial equilibrium.
- There is still a strong similarity between the dense gas traced by (super)virialized 13CO and the subvirial stellar velocities.
- Trends in the stellar velocities as a function of position may be indicative of the dominance of large-scale turbulent modes, rather than evidence of global collapse or cloud rotation.
As an extension, the early dynamical evolution of the star clusters is investigated by (Allison et al. 2010) based on this simulation and observational result that the star clusters form subvirial (cool). The authors notice that the warmer and less substructured a cluster is initially, the longer timescale the mass segregation takes. Specifically, defining the virial ratio as Q = T/U, the total energy is E=−ηGM2/R(1− Q), and R0/Rf= η0/ηf2(1 − Q0). R0/Rf must be as large as possible to cause the maximum degree of collapse. This implies that 2(1 − Q0) needs to be large (Allison et al. 2010).