STELLAR KINEMATICS OF YOUNG CLUSTERS IN TURBULENT HYDRODYNAMIC SIMULATIONS

The majority of stars form in clusters of 100 or more members (Lada & Lada 2003). The dynamics and distribution of stars within clusters have important ramifications for planet formation, stellar multiplicity, mass segregation, cluster disruption, and any other property of stars and planets that can be influenced by interactions between stars in the dense environment of a young cluster. For example, Adams & Laughlin (2001) show that close encounters with neighboring stars can disrupt the outer parts of planetary systems (however, see Megeath et al. 2008). At an earlier stage, protoplanetary disks in dense clusters are subject to dynamical disruption (e.g., Spurzem et al. 2009), photoevaporation (e.g., Adams et al. 2004), and contamination with short-lived radionuclides (e.g., Gounelle et al. 2009; Takigawa et al. 2008; Ouellette et al. 2007). On larger scales, radiation from neighboring stars in a cluster is an important feedback effect that can modify the stellar initial mass function (Fatuzzo & Adams 2008; Offner et al. 2009), while gravitational interactions among stars lead to mass segregation (Bonnell & Davies 1998; McMillan et al. 2007; Allison et al. 2009). All of these phenomena depend critically on the kinematics of newly formed stars, since this determines the frequency of star–star encounters. Kinematics also influence the fraction of star clusters that remain bound rather than dissolving into the field (e.g., Baumgardt & Kroupa 2007).

Observations of the kinematics of young stars are challenging because the stars spend much of their youth enshrouded by gas, making spectroscopy difficult. However, measurements have been possible in a few regions, most notably the Orion Nebula Cluster (ONC), a particularly interesting and well-studied young star cluster containing thousands of optically revealed stars as well as traces of remnant molecular gas. The ONC is observed to be centrally concentrated and mass segregated, with the most massive stars residing in the center (e.g., Hillenbrand & Hartmann 1998). Radial velocity measurements of stars within the ONC indicate a substantial stellar velocity gradient and significant correlation between the local gas and star velocities (Fűrész et al. 2008; Tobin et al. 2009). As a result, Tobin et al. (2009) suggest that the ONC is kinematically young and not dynamically relaxed, which they interpret as requiring that the parent cloud be in a state of dynamical collapse (Hartmann et al. 2001; Bonnell et al. 2003), rather than evolving in quasi-equilibrium (Tan et al. 2006; Krumholz & Tan 2007; Huff & Stahler 2007). Proszkow et al. (2009) use N-body simulations to demonstrate that collapsing clusters with subvirial initial conditions can indeed generate such observational signatures.

However, it is far from clear how to extrapolate from observations of stellar kinematics to the conditions within the host giant molecular cloud (GMC), particularly since not all components of a GMC share the same velocity dispersion. Instead, the dispersion of the line-of-sight (centroid) gas velocity decreases with density (e.g., Blitz & Stark 1986; Heyer et al. 2001). Gas velocity dispersions in prestellar cores are comparable to the sound speed and are much less than the global cloud virial velocity (Goodman et al. 1993; André et al. 2007; Kirk et al. 2007; Rosolowsky et al. 2008). Furthermore, the dispersion of the bulk velocities of the cores themselves is found to be subvirial (André et al. 2007), indicating the relative motion between cores is small. Molecular gas in star-forming regions such as Perseus and Ophiuchus also exhibits large-scale linear velocity gradients similar to those observed in the Orion stellar population (André et al. 2007; Kirk et al. 2007; Rosolowsky et al. 2008). Perseus, in particular, resembles an elongated filament with a gradient of ∼0.3 km s⁻¹ pc⁻¹ along its length, similar to the ONC.

These observations of dense gas suggest that one cannot simply assume that, because the stars are subvirial, the parent cloud must be too. Indeed, for gas the opposite appears to be the case. All of these subvirial signatures for dense components coexist with large-scale velocity dispersions that are consistent with virial balance (e.g., Blitz et al. 2007; Heyer et al. 2009; Bolatto et al. 2008). Simulations of turbulent, virialized clouds are able to reproduce the observed kinematics of dense gas (Padoan et al. 2001; Walsh et al. 2004; Offner et al. 2008b), while simulations of clouds undergoing subvirial collapse have trouble producing quiescent cores like those we observe (Offner et al. 2008b). Nor do clouds appear to convert large fractions of their mass into stars in a single dynamical time, as expected for a subvirial collapse (e.g., Klessen & Burkert 2000; Bonnell et al. 2003). Instead, observations both of individual star-forming regions (Evans et al. 2009) and entire galaxies (Krumholz & Tan 2007) show that, even for gas as dense as the cloud that formed the ONC, ≲5% of the mass per dynamical time can form stars.

The complicated relationship between the kinematics of dense and diffuse gas in molecular clouds suggests that it would be profitable to perform a systematic investigation of how the kinematics of newborn stars relate to the properties of the clouds from which they formed. Such an investigation is also of broader interest, since, as noted above, the kinematics of young clusters have a strong influence on planet formation, mass segregation, and other aspects of star and star cluster evolution. Thus, it is important to develop first-principles predictions for the expected motion of newborn stars as a function of parent cloud properties, a subject that has not been investigated prior to this work. In this Letter, we use turbulent, self-gravitating, adaptive mesh refinement (AMR) hydrodynamics simulations of star cluster formation to study the properties of the star clusters that form in clouds in a variety of dynamical states. In particular, we examine the stellar velocity dispersion, projected velocities, and velocity gradients. We compare the simulation results with recent observations of the ONC and discuss the relevance to current modeling of star clusters. In Section 2 we describe our simulation methods and parameters. In Section 3 we present our results, and we summarize our findings in Section 4.

For our study, we perform three simulations of an isothermal, non-magnetized molecular cloud using the ORION AMR code (e.g., Truelove et al. 1998; Klein 1999). We reexamine the driven and decaying turbulent simulations described in Offner et al. (2008a) and perform an additional simulation with higher resolution and different turbulent conditions.

We generate the initial state of the runs by driving turbulent motions in the gas for two crossing times without gravity. At t = 0, a turbulent steady state has been achieved and the gravitational energy is comparable to the turbulent energy in the box. The cloud virial parameter is then

$$\begin{equation} \alpha = \frac{5 \sigma ^2R}{GM } = 1.67 \end{equation} \tag{ 1 }$$

(Bertoldi & McKee 1992). In the simulation we refer to as U1, energy injection is halted after the initial driving phase and the turbulence decays on a dynamical timescale, initiating a global collapse of the type posited by Tobin et al. (2009). In the simulation D1, we continue turbulent driving so that the Mach number and virial parameter are approximately constant; this is a very rough representation of the behavior expected in a cloud where feedback drives turbulence and prevents global collapse (e.g., Li & Nakamura 2006; Nakamura & Li 2007). The third simulation, D2, begins with the same physical initial conditions, but we adopt a driving pattern with k = 1..2 and drive the turbulence at a constant energy injection rate. As gravitational motions become significant, the net kinetic energy in the box increases and the virial parameter rises. The parameters for the three runs are summarized in Table 1.

Table 1. Simulation Parameters^a

Run	N30 b	Δxminc	Driving Method	kmin..kmax	${\mathcal M}_f$d	αfe
U1	1283	200	No driving	3..4	3.9	1.1
D1	1283	200	Constant $\mathcal {M}$	3..4	8.5	1.7
D2	5123	100	Constant dE/dt	1..2	12.0	3.4

Notes. ^aAll simulations have the same initial L₀, ρ₀, T₀, ${\mathcal M}$. ^bNumber of cells on the base grid. ^cMinimum cell size in AU. ^dFinal three-dimensional Mach number at 1t_ff. ^eFinal virial parameter at 1t_ff.

Download table as:  ASCIITypeset image

Isothermal self-gravitating gas is scale-free, so it is possible to scale the simulation values to observed cloud properties using the characteristic Jeans length and Mach number (see scaling relations in Offner et al. 2008b). However, the main quantity of interest, the ratio of the stellar dispersion to the gas dispersion, is independent of the physical values we adopt. To estimate physical stellar velocities, we adopt a gas temperature of T = 20 K and a gas density of ρ = 9.74 × 10⁻²¹ g cm³, which are similar to the conditions observed in Orion (Johnstone et al. 2001, 2006). This corresponds to a box length of 2 pc and a cloud mass of 1150 ${M}_{\mathord \odot }$.

We run the simulations with gravity for one mean-density free-fall time, t_ff. We introduce a sink particle whenever the Jeans condition is exceeded on the finest AMR level (Krumholz et al. 2004). These particles represent stars, serving as an indicator of the stellar cluster velocity dispersion imparted by the initial turbulent conditions.

Our simulations neglect radiative feedback, which likely has negligible influence on the initial stellar velocity. The simulations also neglect magnetic fields, which may affect the details of the simulated turbulence and exact distribution of gas densities. However, strong magnetic fields are more likely to further constrain the motion of the dense gas and forming stars rather than enhance it (relative to non-magnetized gas).

In this section, we analyze the distribution and velocity information of the stars that form in each simulation. The simulation with decaying turbulence, U1, produces only 20 stars. D2 with the highest amount of turbulence and highest resolution forms 74 stars. In our simulations we do not resolve close binaries, and so in some cases the stellar properties represent the net mass and momentum of a stellar system. At the final time, the runs have efficiencies of ∼10%–15%. Because we do not include protostellar winds, we neglect the ejection of mass from prestellar cores, and thus likely overestimate both the stellar masses and the overall efficiency by a factor of 3–4 (Enoch et al. 2008; Rathborne et al. 2009). However, this shift in stars' masses is unlikely to affect their kinematics.

3.1. Stellar Velocity Dispersion

Figure 1 shows the three-dimensional stellar velocity dispersion as a function of time. The larger driving scale in D2 works to delay the onset of collapse so that stars form much later than in the D1 and U1 calculations. For U1 and D1, the stellar velocity dispersion is sonic to transonic. Once a sufficient number of stars are formed, the velocity dispersion trend appears to be mostly uncorrelated with the simulation Mach number. In the undriven case, the stellar velocity dispersion declines slightly, while remaining nearly fixed in the simulation with constant Mach number driving. In D2, the dispersion decreases slightly with increasing Mach number. However, in all cases, the stellar dispersion is significantly less than the global gas velocity dispersion: σ_*/σ_gas ≃ 0.2. This corresponds to final virial parameters of ∼0.10, 0.09, and 0.10 for U1, D1, and D2, respectively, using the total mass in stars and gas. If all the remaining gas were to disperse, leaving only stellar mass, then the virial parameter of the stars would be measured to be ∼0.8, 0.6, and 1.1 for the runs U1, D1, and D2, respectively.

Zoom In Zoom Out Reset image size

This indicates that only if the initial gas cloud is significantly supervirial does the virial parameter of the stellar component become greater than unity, and that star clusters forming in turbulent virialized clouds naturally begin with subvirial velocities. Consequently, the subvirial nature of the stellar velocity dispersion is not necessarily an indication of collapse. Indeed, Proszkow et al. (2009) find that they are able to reproduce the kinematic signature of Orion assuming an initial virial ratio of ∼0.1 (for stars plus gas), the value produced by simulation D1, in which the gas cloud is turbulent and virialized. Similarly, Allison et al. (2009) assume an initial cluster virial parameter of 0.8, much less subvirial than we find for our simulated clouds with virial parameters of α ≃ 1. If the ratio of the stellar dispersion to the gas dispersion influences how quickly the cluster undergoes mass segregation, then most clusters with substantially subvirial initial stellar velocities dispersions may segregate on short timescales.

3.2. Projected Stellar Velocities

In this section, we investigate the stellar velocities projected along the line of sight to the cluster. For the most detailed examination, we use the final time of the D2 simulation, which has the best statistics. Figure 2 shows the velocities projected along the three cardinal sight lines. The velocities are somewhat clustered with only a couple instances of large deviations. We have further binned the projected velocities into groups containing eight stars as a function of position and we have overlaid error bars indicating the velocity dispersion in each bin. These plots look similar to Figure 9 in Proszkow et al. (2009) of the ONC, which shows evidence of a velocity gradient across the cluster.

Zoom In Zoom Out Reset image size

We derive the stellar velocity gradient by fitting the line-of-sight velocities with the function V = V₀ + ∇V_xX + ∇V_yY. Although the data may not be exactly linear in nature (e.g., the ONC exhibits a linear gradient mainly in the cluster center), we can characterize the change in the velocities across the cloud. Since we use periodic boundary conditions the velocity bins tend to wrap around the domain. Thus, we can better compare to the data by dividing each projection into smaller boxes and investigating the velocity trend in each box separately. In Figure 3, we show the velocity gradients, |∇V| = (∇V²_x + ∇V²_y)^1/2, along each projection for the entire domain and for boxes half and a quarter of the domain. We include gradients at three different times and plot ONC data from Tobin et al. (2009) fitted in the same manner. To limit statistical fluctuations, we only include gradients fitted to four or more stars. The full and half-box fits have median velocity gradient magnitudes of 1 km s⁻¹ pc⁻¹ and 2.1 km s⁻¹ pc⁻¹, respectively. 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. For comparison, Proszkow et al. (2009) find that the ONC has a radial velocity gradient of 1.1 km s⁻¹ pc over half of its length, which corresponds to an overall gradient of ∼0.5 km s⁻¹ pc. This value is similar to what we find for our simulated virialized clusters. Thus, velocity gradients are not necessarily a hallmark of a subvirial turbulence or a strongly collapsing cloud. Instead, large-scale velocity gradients appear to be a signature of the fact that both our simulations and real molecular clouds have most of their turbulent power in large-wavelength modes.

Zoom In Zoom Out Reset image size

We find that the linear fits of the projected velocities fluctuate somewhat as a function of time, but they tend to fall within 0.2–2 km s⁻¹ pc⁻¹ on the largest scales. These values are similar to gradients in the gas centroid velocity observed in regions such as Perseus and Ophiuchus (André et al. 2007; Rosolowsky et al. 2008).

3.3. Gas Velocities

Observations of clusters frequently find that the star and gas velocities are correlated (e.g., Fűrész et al. 2008). We can reproduce such observations by making simulated observations of the gas in particular molecular tracers. Following the procedure of Offner et al. (2008b), we generate a position–position–velocity cube of the gas intensity in the ¹³CO (J = 1 → 0) molecular line along one line of sight in the D2 run. We solve for the line emission, assuming that the gas is in statistical equilibrium and that radiative pumping is negligible.⁴ We adopt a cloud distance of 400 pc and 26'' beam size. To model the telescope resolution, we smear each velocity channel with the Gaussian beam.

Figure 4 shows the line-of-sight ¹³CO gas and star velocities; following Fűrész et al. (2008, their Figure 10), to produce this figure from the three-dimensional position–position–velocity cube, we have integrated along one direction on the plane of the sky, so the intensity shown in a given (x, v) pixel is the sum over all y values for that x and v. Our simulated image shows a strong correlation between gas and star velocities, similar to that seen in the star and gas velocity maps presented by Fűrész et al. (2008). For the projection shown in Figure 4, the mean stellar velocity is ∼−0.8 km s⁻¹, somewhat to the left of the mean ¹³CO gas velocity at 0.003 km s⁻¹. Most of the disagreement in the histogram occurs because not all regions traced by ¹³CO have formed stars. Projections along other directions show similar correlations between gas and star velocities.

Zoom In Zoom Out Reset image size

The strong correlation between gas and star velocities might seem surprising, given that the gas and stars have such different velocity dispersions. The explanation is that the anti-correlation between density and velocity found in turbulent gas serves to artificially enhance the gas–star correlation. Figure 4 is constructed by summing over the y direction, so each pixel represents a density-weighted average velocity. The dispersion, however, is a density-weighted average of velocity squared, which, since velocity and density are anti-correlated, is much more heavily weighted toward low-density, high-velocity gas than what is shown in Figure 4. This low-density gas is also that which is least correlated with the stars. The bias is exacerbated because ¹³CO is thermally excited only in relatively dense gas that is closely associated with forming stars. It is subthermally excited in the lower density gas that dominates the dispersion. Thus, both averaging and excitation serve to improve the correlation between ¹³CO emission and stellar velocities, while at the same time obscuring the low-density gas that carries the bulk of the kinetic energy.

In this Letter, we have investigated the kinematic properties of stars forming in simulated turbulent self-gravitating clouds. These calculations represent the initial conditions of stars in clusters prior to dynamic evolution and cloud dispersal. We demonstrate that subvirial stellar velocity dispersions arise naturally from clouds in virial equilibrium. Thus, long-lived clouds satisfying a quasi-steady-state are not ruled out by observations of clusters with small stellar velocity dispersions. We find that gas virial parameters of order unity produce a star-to-gas velocity dispersion ratio of approximately 0.2.

We find that turbulent initial conditions easily generate stellar velocity gradients of the magnitude observed in star-forming regions. Using a planar fit, we demonstrate typical gradients of 1 km s⁻¹ pc⁻¹, similar to the observed gradient in the ONC. We also show that for young clusters there is a strong similarity in the simulations between the dense gas traced by ¹³CO and the star velocities. Thus, 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 suggested by some authors.

Support for this work was provided by the US Department of Energy at the Lawrence Livermore National Laboratory under contract B-542762 (S.S.R.O., C.E.H.), the Alfred P. Sloan Foundation (M.R.K.), NASA/JPL through the Spitzer Theoretical Research Program (M.R.K.), and the National Science Foundation through grants AST-0807739 (M.R.K.) and AST-0901055 (S.S.R.O.). Computational resources were provided by the NSF San Diego Supercomputing Center through NPACI program grant UCB267.