Rotation in young massive star clusters

TL;DR

This study uses hydrodynamical simulations to show that massive star clusters commonly inherit significant rotation from their parent gas, supporting hierarchical assembly as a key formation process.

Contribution

It demonstrates through simulations that hierarchical merging leads to rotation in massive star clusters, a novel insight into their formation dynamics.

Findings

Massive star clusters inherit rotation from parent gas.

Clusters with mass > few hundred Msun show clear rotation.

Simulated clusters are highly elliptical and subvirial.

Abstract

Hydrodynamical simulations of turbulent molecular clouds show that star clusters form from the hierarchical merger of several sub-clumps. We run smoothed-particle hydrodynamics simulations of turbulence-supported molecular clouds with mass ranging from 1700 to 43000 Msun. We study the kinematic evolution of the main cluster that forms in each cloud. We find that the parent gas acquires significant rotation, because of large-scale torques during the process of hierarchical assembly. The stellar component of the embedded star cluster inherits the rotation signature from the parent gas. Only star clusters with final mass < few X 100 Msun do not show any clear indication of rotation. Our simulated star clusters have high ellipticity (~0.4-0.5 at t=4 Myr) and are subvirial (Q_vir<~0.4). The signature of rotation is stronger than radial motions due to subvirial collapse. Our results suggest that rotation is common in embedded massive (>~1000 Msun) star clusters. This might provide a key observational test for the hierarchical assembly scenario.

Tables

Table 1. Table 1: Initial conditions.

Run	$M$ [M⊙]	$R$ [pc]	$N_{\mathrm{gas}}$ $[\times {10}^6]$
A	$4.3\times {10}^4$	8.8	10
B	$4.3\times {10}^4$	8.8	10
C	$1.0\times {10}^4$	5.4	10
D	$1.7\times {10}^3$	3.0	2

Table 2. Table 2: Main properties of simulated star clusters.

Run	$t$	$t_{\ast }$	$M_{\ast }$	$r_{\mathrm{hm}}$	$Q_{\mathrm{vir}}$	$r_1,$ $r_2,$ $r_3$	$ϵ$	${\sigma }_x$	$⟨\left|v_x\right|⟩/\sigma {}_x$	$L_z$
	[Myr]	[Myr]	[M⊙]	[pc]		[pc], [pc], [pc]		[km s-1]		[M⊙ km2 s-1]
A	2.5	$0.7\pm 0.2$	865	0.012	0.31	0.037, 0.022, 0.009	0.75	6.9	0.82	$2.8\times {10}^{15}$
A	4.0	$1.3\pm 0.6$	7968	0.085	0.36	0.55, 0.39, 0.32	0.42	8.9	0.77	$1.6\times {10}^{17}$
B	2.5	$0.7\pm 0.2$	1060	0.008	0.22	0.024, 0.013, 0.007	0.71	6.5	0.85	$2.7\times {10}^{15}$
B	4.0	$1.6\pm 0.6$	11874	0.049	0.36	0.70, 0.35, 0.33	0.53	11.8	0.73	$1.8\times {10}^{17}$
C	2.5	$0.8\pm 0.1$	202	0.006	0.20	0.018, 0.011, 0.004	0.75	2.9	0.77	$2.0\times {10}^{14}$
C	4.0	$1.3\pm 0.7$	2057	0.070	0.32	0.65, 0.47, 0.40	0.39	5.0	0.76	$1.2\times {10}^{16}$
D	2.5	$0.3\pm 0.1$	15	$3\times {10}^{-4}$	0.11	0.0010, 0.0005, 0.0003	0.71	–	–	$5.4\times {10}^{11}$
D	4.0	$0.7\pm 0.5$	277	0.074	0.16	0.51, 0.41, 0.34	0.51	3.5	0.70	$6.5\times {10}^{14}$

Table 3. Table 3: Comparison of runs C and D with the additional test realizations.

Run	$t$	$t_{\ast }$	$M_{\ast }$	$r_{\mathrm{hm}}$	$Q_{\mathrm{vir}}$	$r_1,$ $r_2,$ $r_3$	$ϵ$	${\sigma }_x$	$⟨\left|v_x\right|⟩/\sigma {}_x$	$L_z$
	[Myr]	[Myr]	[M⊙]	[pc]		[pc], [pc], [pc]		[km s-1]		[M⊙ km2 s-1]
D	4.0	$0.7\pm 0.5$	277	0.074	0.16	0.51, 0.41, 0.34	0.51	3.5	0.70	$6.5\times {10}^{14}$
D1	4.0	$0.9\pm 0.6$	402	0.28	0.36	1.21, 0.58, 0.41	0.66	1.6	0.75	$3.5\times {10}^{15}$
D2	4.0	$0.8\pm 0.5$	379	0.37	0.31	1.88, 0.54, 0.48	0.75	2.1	0.69	$2.1\times {10}^{15}$
D3	4.0	$1.1\pm 0.6$	775	0.18	0.44	0.93, 0.68, 0.58	0.37	2.4	0.75	$6.1\times {10}^{15}$
D4	4.0	$1.1\pm 0.6$	529	0.11	0.40	0.51, 0.44, 0.40	0.22	2.5	0.72	$2.8\times {10}^{15}$
C	4.0	$1.3\pm 0.7$	2057	0.070	0.32	0.65, 0.47, 0.40	0.39	5.0	0.76	$1.2\times {10}^{16}$
C1	4.0	$1.4\pm 0.5$	1850	0.064	0.36	0.53, 0.41, 0.16	0.69	5.8	0.78	$1.6\times {10}^{16}$

Equations

$$I_{\text{rot}}=\sum _i{m}_i\left[\begin{array}{ccc}{x}_i^2& x_i{y}_i& x_i{z}_i\\ x_i{y}_i& y_i^2& y_i{z}_i\\ x_i{z}_i& y_i{z}_i& z_i^2\end{array}\right],$$

$$r_i=\beta \sqrt{\frac{5{\lambda }_i}{M_{\ast }}},i=1,2,3,$$

$$t_{\mathrm{rlx}}\sim 0.5\mathrm{Myr}{\left(\frac{M_{\ast }}{8000{\mathrm{M}}_{\odot }}\right)}^{1/2}{\left(\frac{r_{\mathrm{hm}}}{0.1\mathrm{pc}}\right)}^{3/2}.$$

Full text

1 Introduction

Hydrodynamical simulations of turbulent molecular clouds show that star clusters form from the hierarchical merger of several sub-clumps (e.g. Bonnell et al. 2003; Bate 2009; Federrath et al. 2010; Bonnell et al. 2011; Dale & Bonnell 2011; Girichidis et al. 2011). These sub-clumps, composed of stars and gas, develop preferentially along dense filaments, and then aggregate into dense regions at the junction between filaments.

Whether real-life star clusters develop hierarchically or monolithically is still matter of debate. The simulations by Banerjee & Kroupa (2015) support a monolithic formation or a prompt ($<1$ Myr) assembly for the young dense star cluster NGC 3603. In contrast, the Hubble Space Telescope data of 30 Doradus indicate the existence of two distinct stellar populations, which might be the signature of an ongoing merger (Sabbi et al., 2012). Overall, the properties of simulated star clusters deserve to be further investigated and compared against observations, to probe the robustness of the hierarchical formation scenario.

In this paper, we discuss a set of high-resolution simulations of star cluster formation from turbulent molecular clouds. In our simulations, star clusters form from the hierarchical assembly of several clumps. We show that most star clusters, especially the most massive ones, develop rotation as a consequence of the torques arising from hierarchical merging.

Rotation was found to be common in globular clusters (Pryor et al., 1986; van Leeuwen et al., 2000; Anderson & King, 2003; Pancino et al., 2007; Anderson & van der Marel, 2010; Bellazzini et al., 2012; Bianchini et al., 2013; Fabricius et al., 2014; Kimmig et al., 2015; Lardo et al., 2015; Lee, 2015). Several theoretical studies investigate the origin and discuss the importance of rotation for models of globular clusters (Bekki, 2010; Mastrobuono-Battisti & Perets, 2013; Vesperini et al., 2014; Bianchini et al., 2015; Gavagnin et al., 2016; Mastrobuono-Battisti & Perets, 2016). In particular, Vesperini et al. (2014) show that globular cluster rotation might arise as a consequence of violent relaxation in the tidal field of the host galaxy, while Gavagnin et al. (2016) suggest that it may be the result of the merger between proto-globular clusters.

Signatures of rotation have been claimed also in few young (e.g. R136 in the Large Magellanic Cloud, Hénault-Brunet et al. 2012) and intermediate-age star clusters (e.g. NGC 1846, Mackey et al. 2013, GLIMPSE-CO1, Davies et al. 2011). Other observations highlight the existence of velocity gradients in young star forming regions, which might be connected with rotation (e.g. André et al. 2007; Rosolowsky et al. 2008; Tobin et al. 2009; Cottaar et al. 2015; Tsitali et al. 2015; Rigliaco et al. 2016).

In agreement with other recent theoretical studies (e.g. Lee & Hennebelle 2016a, b), our results show that rotation should be very common in the first stages of star cluster formation, and suggest future observational campaigns should look for this signature. This paper is organized as follows. In Section 2, we describe the methodology; our results are presented in Section 3; in Section 4 we discuss our results, and in Section 5 we draw our conclusions.

2 Methods

We generate models of turbulence-supported molecular clouds. The clouds are initially spherical, with homogeneous density. They are seeded with supersonic turbulent velocities and marginally self-bound (see Hayfield et al. 2011; Mapelli et al. 2012, 2013; Mapelli & Trani 2016). To simulate interstellar turbulence, the velocity field of the cloud is generated on a grid as a divergence-free Gaussian random field with an imposed power spectrum $P(k)$, varying as $k^{-4}$. This yields a velocity dispersion $\sigma{}(l)$, varying as $l^{1/2}$, chosen to agree with the Larson (1981) scaling relations.

The clouds are simulated with the smoothed-particle hydrodynamics (SPH) code gasoline (Wadsley et al., 2004), upgraded with the Read et al. (2010) optimized SPH (OSPH) modifications, to address the SPH limitations outlined by Agertz et al. (2007). We include radiative cooling in all our simulations. The radiative cooling algorithm is the same as that described in Boley (2009) and in Boley et al. (2010). According to this algorithm, the divergence of the flux is $\nabla{}\cdot{}F=-(36\,{}\pi{})^{1/3}\,{}s^{-1}\sigma{}\left({\rm T}^{4}-{\rm T}^{4}_{\rm irr}\right)\,{}(\Delta{}\tau{}+1/\Delta{}\tau{})^{-1}$, where $\sigma{}=5.67\times{}10^{-5}$ erg cm*-2* s*-1* K*-4* is the Stefan’s constant, ${\rm T}_{\rm irr}$ is the incident irradiation, $s=\left(m/\rho{}\right)^{1/3}$ and $\Delta{}\tau{}=s\,{}\kappa{}\,{}\rho{}$, for the local opacity $\kappa{}$, particle mass $m$ and density $\rho{}$.

D’Alessio et al. (2001) Planck and Rosseland opacities are used, with a 1 $\mu{}$m maximum grain size. Such opacities are appropriate for temperatures in the range of a few Kelvins up to thousands of Kelvins. In our simulations, the irradiation temperature is $T_{\rm irr}=10$ K everywhere. The only feedback mechanism we account for is compressional heating of gas. We include no treatment for photoionisation, stellar winds or supernovae, to keep our analysis as simple as possible. The impact of stellar winds on star cluster formation was found to be relatively minor (Dale et al., 2013; Dale et al., 2015). Moreover, several studies suggest that star formation proceeds on a timescale which is much faster than the timescale for photoionising radiation or supernovae to remove the gas (Kruijssen, 2012; Kruijssen et al., 2012; Dale & Bonnell, 2011; Dale et al., 2012; Iffrig & Hennebelle, 2015). The impact of feedback on our results will be assessed in a follow-up study.

In our simulations, star formation is modelled through the sink-particle technique. Sink particles have been implemented according to the criteria described in Bate et al. 1995 (hereafter B95). We adopt a sink accretion radius $r_{\rm acc}=2\,{}\tilde{\epsilon{}}$, where $\tilde{\epsilon}$ is the softening length ($\tilde{\epsilon}=10^{-4}-10^{-3}$ pc, depending on the resolution of the simulation, see e.g. Mapelli et al. 2012). A gas particle is considered a sink candidate if its density is above the threshold $\rho{}_{\rm th}=10^{-17}$ g cm*-3*, corresponding to $n_{\rm th}\sim{}10^{7}$ atoms cm*-3* (we checked that the sink mass function does not depend on this choice significantly). If gas particles inside the accretion radius of the sink candidate satisfy B95 criteria111B95 criteria for converting gas particles into sink particles require i) that the thermal energy of particles inside $r_{\rm acc}$ is $E_{\rm th}\leq{}0.5\,{}E_{\rm g}$, where $E_{\rm g}$ is the magnitude of the gravitational energy of the particles; ii) that $E_{\rm th}/E_{\rm g}+E_{\rm r}/E_{\rm g}\leq{}1$, where $E_{\rm r}$ is the rotational energy of particles; iii) that the total energy of particles inside $r_{\rm acc}$ is negative., then the candidate becomes a sink particle. A similar procedure is followed to decide whether a particle which has already become a sink will accrete gas particles222According to B95 criteria a gas particle within $r_{\rm acc}$ will be accreted by the sink particle if i) the gas particle is bound to the sink; ii) the specific angular momentum of the particle about the sink is less than required to form a circular orbit at $r_{\rm acc}$; iii) the gas particle is more tightly bound to the considered sink particle than to other sink particles.. The whole procedure is clearly affected by the resolution of the simulation. We checked that this dependence has no effect on the results of this paper (i.e. the feature of rotation), by running the same simulation with lower resolution.

We simulate four molecular clouds, with initial mass $M$ spanning from 1700 to $4.3\times{}10^{4}$ M⊙, and initial radius $R$ in the $3.0-8.8$ pc range (see Table 1). $R$ and $M$ were chosen so that the initial density of the cloud is always the same $\rho{}\sim{}10^{-21}$ g cm*-3* (i.e. $n\sim{}250$ cm*-3* when assuming molecular weight $\mu{}=2.46$). Thus, the free fall time-scale of the simulated clouds is $t_{\rm ff}=[3\,{}\pi{}/(32\,{}G\,{}\rho{})]^{0.5}\sim{}2$ Myr, where $G$ is the gravity constant. Each run was integrated for $\sim{}2\,{}t_{\rm ff}=4$ Myr.

Cloud A and B are two different realizations (with two different random seeds for the turbulence) of the same molecular cloud, with $M=4.3\times{}10^{4}$ M⊙ and $R=8.8$ pc. Clouds C and D have $M=10^{4}$ M⊙ and $1.7\times{}10^{3}$ M⊙ ($R=5.4$ and 3.0 pc), respectively. Clouds A, B and C were simulated with $10^{7}$ equal-mass gas particles, while the low-mass cloud D was simulated with $2\times{}10^{6}$ particles. Thus, the particle mass is always \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}{}4\times{}10^{-3} M⊙.

3 Results

Star formation begins at time $t=1.3\,{}{\rm Myr}\sim{}0.6\,{}t_{\rm ff}$ in clouds A and B, at $t=1.5\,{}{\rm Myr}\sim{}0.7\,{}t_{\rm ff}$ in cloud C, and at $t=1.9\,{}{\rm Myr}\sim{}0.9\,{}t_{\rm ff}$ in cloud D. The top panel of Figure 1 shows the projected density of gas in simulation A at $t=2.5$ Myr (i.e. $\sim{}1.2\,{}t_{\rm ff}$). The cloud has developed the well-known filamentary structures (Bonnell et al., 2003; Bonnell et al., 2011), as an effect of the interplay between self-gravity and turbulence (Moeckel & Burkert, 2015; Federrath, 2016). The overall behaviour is similar in the other simulated clouds.

The bottom panel of Figure 1 is a zoom of the region where the main star cluster is assembling. The gas velocity maps reported in Fig. 2 indicate that this region rotates. Rotation is still apparent in the gas even at the end of the simulation ($t=4$ Myr). In our simulations, rotation arises from the torques exerted by gas filaments and clumps while they merge with the main proto-cluster structure.

The rotation of the parent molecular gas leaves a strong imprint on the star cluster. Fig. 3 shows the velocity field of the main star cluster in simulations A, B, and C at time $t=2.5$ Myr. Each point in the Figure is a sink particle. The star clusters are projected in the $xy$ plane, that is defined as the plane perpendicular to the angular momentum vector (which is thus aligned along the $z$ axis). Simulation D is not shown because the main star cluster is not well defined at $t=2.5$ Myr (it has a total mass of only $\sim{}15$ M⊙). All star clusters in Fig. 3 show a clear signature of rotation in $v_{\rm x}$ and $v_{\rm y}$.

Similarly, Fig. 4 shows the velocity field of the main star cluster in simulations A, B, and C, at time $t=2.5$ Myr, when projected in the $yz$ plane. With respect to this projection, the velocity along the $x$ axis ($v_{\rm x}$), shown by the colour-coded map, can be interpreted as the line-of-sight velocity. Again, the rotation of the three star clusters is apparent. Moreover, this projection shows that the star clusters are extremely flattened. The reason is that they form from a disk-like rotating gas structure, as shown in Figures 1 and 2. In the first stages of their assembly, the star clusters still retain memory of their parent gas disc.

At this stage ($t=2.5$ Myr) the star clusters are still extremely small and compact. The main structural and kinematic parameters of the simulated star clusters are listed in Table 2. The total bound mass to the star clusters at $t=2.5$ Myr ranges from only $\sim{}15$ M⊙ up to $\sim{}1060$ M⊙ for runs D and B, respectively. The half mass radius $r_{\rm hm}$ is smaller than $\sim{}0.012$ pc. The virial ratio $Q_{\rm vir}=K/|W|$ (where $K$ and $W$ are the kinetic and potential energy of the stars, respectively333In the definition of $Q_{\rm vir}$ we do not include the contribution of gas particles, as suggested by Kruijssen (2012) and Farias et al. (2015).) is $\sim{}0.2-0.3$, indicating that our protoclusters are markedly subvirial. At time $t=2.5$ Myr since the beginning of the simulation, the average age of stars444The age of a star is defined as the time elapsed since the formation of the sink particle. in the main star cluster is $t_{\ast}\sim{}0.7-0.8$ Myr in runs A, B, and C, and $t_{\ast}{}\sim{}0.3$ Myr in run D, with a non-negligible age spread ($\sim{}0.1-0.2$ Myr).

The ellipticity ($\epsilon$) allows us to quantify the flattening of the simulated star clusters. To estimate the ellipticity, we first calculate the rotational inertia matrix of sink particles with respect to their centre of gravity:

[TABLE]

where $i$ is the index of sink particles. The three eigenvalues of this matrix ($\lambda_{1},\lambda_{2},\lambda_{3}$) give the star cluster size in three orthogonal directions:

[TABLE]

where $M_{\ast}$ is the total mass of sinks. The factor $5$ comes from the assumption that the sink mass is uniformly distributed in space, and a correction factor $\beta\geq 1$ accounts for the fact that the mass distribution might not be uniform but rather centrally concentrated (see Lee & Hennebelle 2016a). The values of $r_{1},\,{}r_{2},\,{}$ and $r_{3}$ shown in Table 2 were calculated assuming $\beta{}=1$ for simplicity.

We thus define $\epsilon=1-\sqrt{\lambda_{3}/\lambda_{1}}$, where $\lambda_{1}$ and $\lambda_{3}$ are the maximum and minimum eigenvalue, respectively. $\epsilon{}=0$ means that the star cluster is spherical, while $\epsilon\sim{}1$ means that the star cluster is nearly flat. For the protostar clusters we find $\epsilon{}\sim{}0.71-0.75$ at $t=2.5$ Myr. Thus, the simulated star clusters are extremely flattened during the early stages of their assembly.

At $t=4$ Myr (i.e. $\sim{}2\,{}t_{\rm ff}$), the signature of rotation is still strong in runs A, B, and C (Figs. 5 and 6), even if the ellipticity has diminished to $\epsilon{}\sim{}0.4-0.5$ (Table 2). In these late evolutionary stages, the star clusters have grown by hierarchical assembly to a size of $\approx{}0.5$ pc (Table 2). The half-mass radius in runs A, B and C is close to $r_{\rm hm}\sim{}0.05-0.1$ pc, and the virial ratio tends to increase ($Q_{\rm vir}\sim{}0.3-0.4$), even if the star clusters are still subvirial. The total final star cluster mass ranges from $\sim{}300$ M⊙ to $\sim{}12000$ M⊙ in runs D and B, respectively.

At $t=4$ Myr since the beginning of the simulation, the average age of stars in the main star cluster is $t_{\ast}\sim{}1.3-1.6$ Myr in runs A, B, and C, and $t_{\ast}{}\sim{}0.7$ Myr in run D, with a significant age spread (the standard deviation being $\sim{}0.5-0.7$ Myr). Thus, most simulated stars are still very young objects at the end of the simulation (Table 2).

Figure 7 shows the rotation curve of star clusters in runs A, B, C, and D at the end of the simulation, compared with the one-dimensional velocity dispersion ($\sigma{}_{\rm x}$). The rotation curves of runs A and B (which are two different random realizations of the same cloud) are very similar to each other, indicating that rotation is a common outcome in star clusters formed from molecular clouds with mass $M\sim{}4\times{}10^{4}$ M⊙. Rotation is apparent even in the star cluster of run C, whose mass is $\sim{}1/5$ of the mass of star clusters in runs A and B.

The star cluster of run D is the only one which does not show a clear signature of rotation. It assembles later and has a total mass of only $\sim{}300$ M⊙, much smaller than the other star clusters. This indicates that the signature of rotation depends on the total stellar mass and thus on the initial gas mass that is converted into stars. We expect rotation to be negligible for star clusters with mass M_{\ast}{}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}{} few $\times{}10^{2}$ M⊙.

Stochastic fluctuations might affect this result, because run D is the less massive cloud we simulated and was sampled with less particles ($N_{\rm gas}=2\times{}10^{6}$) than the other runs. Thus, we check the importance of stochastic fluctuations by running four additional realisations of run D, which are discussed in Appendix A. Stochastic fluctuations are found to be significant: star clusters with mass ranging from $\sim{}300$ to $\sim{}800$ M⊙ form in the five realisations of run D. Our main conclusion still holds, because there is no signature of rotation in those realisations of run D where the final star cluster mass is \lower 2.15277pt\hbox{;\buildrel<\over{\sim};}{}500 M⊙. Only one realisation, with star cluster mass $\sim{}800$ M⊙, shows some hints of rotation. This confirms that rotation is expected to be negligible for star clusters with mass M_{\ast}{}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}{} few $\times{}10^{2}$ M⊙.

Proszkow et al. (2009) point out that the collapse of subvirial clusters can produce kinematic signatures similar to rotation, with a red-shifted and a blue-shifted component. The main difference between the signature of genuine rotation and that of subvirial collapse is that the latter has no preferential plane (and ideally zero angular momentum vector), while the former is maximum in the plane perpendicular to the angular momentum vector.

Indeed, we expect that the effect of subvirial collapse is present even in our simulated star clusters, because they are subvirial and elongated (Table 2). The residuals of blue-shifted and red-shifted structures in the colour-coded map of the velocity along the angular momentum vector ($v_{\rm z}$, right-hand panels of Fig. 5) are the demonstration of this. However, these residuals are much weaker than the signature of rotation in the plane perpendicular to the angular momentum vector (left-hand and central panels of Fig. 5). Moreover, Fig. 8 compares $\bar{v}_{\rm x}(y)$ (i.e. a component of the velocity perpendicular to the angular momentum vector) with $\bar{v}_{\rm z}(y)$ (i.e. the component of the velocity parallel to the angular momentum vector). A velocity gradient (from red-shifted to blue-shifted) is clearly visible only in the $\bar{v}_{\rm x}(y)$ component. Thus, genuine rotation is stronger than the effect of subvirial collapse in our runs A, B, and C.

4 Discussion and caveats

Our simulations indicate that natal rotation must be a common feature in massive embedded star clusters. However, there are several possible issues. First, we might wonder whether the simulations are affected by any bias, which makes rotation stronger. The absence of stellar feedback in our simulations is a limitation for several reasons. Feedback is expected to quench star formation at some point. This affects the final mass of the stellar cluster. The local star formation efficiency (SFE) in our runs is \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}{}50 per cent; we expect SFE to be significantly lower, if feedback becomes important on a time shorter than (or comparable to) the star formation timescale.

Also, we expect that the impact of feedback on rotation depends on its timescale. The first supernovae occur when the age of the massive stars is t_{\ast}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}{}3 Myr, depending on the stellar mass. Similarly, stellar winds by massive stars are particularly strong in the last stages of their life (Chen et al., 2015). In our simulations, rotation is already well developed when $t_{\ast}{}\sim{}1-2$ Myr (see Table 2). Thus, the onset of rotation should not be affected by stellar winds and supernovae. In contrast, photoionization and outflows from protostars affect the evolution of gas even at earlier stages. In particular, outflows from young stars are very common (see e.g. Bally 2016 for a review). They inject energy and momentum into their surroundings. This might counteract the inflow of gas toward the star cluster, reducing the amount of rotation. These possible feedback effects suggest that the minimum star cluster mass to have significant rotation might be higher than what we find in our simulations. Thus, it is essential to include feedback in follow-up simulations.

Another crucial ingredient that is missing in our simulations is the magnetic field. The effect of magnetic field is known to be crucial on the scale of single cores, where magnetic braking can remove the angular momentum of the cloud. A very efficient magnetic braking can even prevent the formation of rotationally supported circumstellar discs (Galli et al., 2006; Lizano & Galli, 2015). Studies on the effects of magnetic fields on star cluster scales (e.g. Federrath 2016) do not address the issue of rotation. However, Lee & Hennebelle (2016a) include a weak ($3-8\,{}\mu{}{\rm G}$) magnetic field in their simulations and still see a clear feature of rotation in the gas, even if they do not discuss the rotation curve of stars. Forthcoming studies of star cluster rotation including both stellar feedback and magnetic fields are needed to address this topic in detail.

An additional caveat we must mention is the importance of initial conditions. Several properties of the simulated star clusters (e.g. mass function and spatial distribution of stars) are known to depend on the initial conditions, the initial density profile of gas and the initial spectrum of turbulence being among the crucial ingredients (Girichidis et al., 2011). Thus, it is important to consider various initial conditions and to check their impact on star cluster rotation.

An important question raised by our results is how long the rotation signature lasts in star clusters. The expulsion of the residual gas, the galactic tidal field, and the dynamics of stars (e.g. two-body relaxation, close encounters between stars and binaries, mass segregation, core collapse, adiabatic cluster expansion) all might affect and possibly quench rotation. Rotation was found to accelerate the dynamical evolution of a star cluster through the transfer of angular momentum outward: mass segregation tends to be faster, and the core collapse occurs at earlier times (Kim et al., 2002; Kim et al., 2004, 2008). Under some circumstances, the gravo-gyro catastrophe might occur: the most massive stars, segregated in the cluster core, speed up and rotate faster than the rest of the cluster (Hachisu, 1979; Einsel & Spurzem, 1999; Kim et al., 2004). Another important consequence of rotation is that stars with high angular momentum (in the tail of the Maxwellian velocity distribution) are more likely to escape from the star cluster with respect to low angular-momentum stars (Agekian, 1958; Shapiro & Marchant, 1976). This leads to a removal of angular momentum from the star cluster, progressively erasing rotation.

At the end of our simulations, the gas is nearly exhausted in the central regions of the main star clusters, with a local SFE $>0.5$, in agreement with previous studies (Bonnell et al., 2011; Kruijssen, 2012). Thus, the star clusters are mainly gas-free and their evolution is driven by two-body encounters between stars. The two-body relaxation timescale can be expressed as (Portegies Zwart et al., 2010)

[TABLE]

Since $t_{\rm rlx}$ is short, the star cluster will evolve very fast because of two-body relaxation and, possibly, binary-single star encounters. Two-body relaxation causes the angular momentum to diffuse outward, removing it from the system (Kim et al., 2002). Moreover, encounters between single stars and binary systems tend to isotropise the distribution of angular momentum (Mapelli et al., 2005). This might be one of the reasons why the star clusters tend to become rounder at the end of the simulations. Another possible explanation for the loss of flattening is the accretion of new gas and stellar clumps, with a different orientation of angular momentum.

Given these considerations, we expect that stellar dynamics tends to suppress the natal rotation. Thus, it is of foremost importance to investigate the impact of stellar dynamics on the newly born star clusters. However, the adopted tree code (employing a second-order leap-frog scheme) is not particularly suited to investigate close stellar encounters in detail. Only few studies attempt to trace the dynamical evolution of star clusters born from hydrodynamical simulations of molecular clouds (Fujii & Portegies Zwart, 2015; McMillan et al., 2015; Fujii, 2015; Fujii & Portegies Zwart, 2016), because the problem is a numerical challenge (e.g. Pelupessy & Portegies Zwart 2012; Hubber et al. 2013). In a forthcoming study, we will trace the dynamics of stars after gas expulsion, adopting a suitable direct-summation $N-$body code.

Comparing our results against observations of young star clusters would be extremely useful to constrain the hierarchical model of star cluster formation. Unfortunately, a direct comparison between our simulations and the available data is complex for several reasons. First, our simulated star clusters are extremely young embedded objects: in Tab. 2 we show that $t_{\ast}{}\sim{}1-2$ Myr at $t=4$ Myr since the beginning of the simulations. $t_{\ast}$ is defined as the average time elapsed since the formation of sink particles. The ‘true’ stellar age might be even smaller than $t_{\ast}$, because the sink particle algorithm cannot distinguish between pre-stellar cores and stars. Moreover, our simulated star clusters are very dense ($>10^{4}$ M⊙ pc*-3*) and compact ($r_{\rm hm}<0.1$ pc). Only few star clusters with comparable age and density exist in the Milky Way and in the Magellanic Clouds (e.g. Portegies Zwart et al. 2010). Moreover, the relatively low line-of-sight velocities ($\sim{}5-10$ km s*-1*) require a very high spectral resolution. Undetected binaries significantly complicate the analysis. For these reasons, observations are rather scanty.

The youngest massive star cluster for which evidence of rotation was found is R136 in the Large Magellanic Cloud (Hénault-Brunet et al., 2012). Our simulated rotation curves (Fig. 7) are qualitatively similar to that of R136 (Hénault-Brunet et al., 2012). However, our simulated star clusters are less massive555The most massive star cluster in our simulations, run B, has a maximum mass of $\sim{}10^{4}$ M⊙, while R136 has a current mass of $\sim{}10^{5}$ M⊙ (Andersen et al., 2009). and much more compact than R136 (they have not yet started the expansion phase after gravitational collapse). Moreover, the sample of Hénault-Brunet et al. (2012) consists of only 36 apparently single O-type stars, of which only 16 are within 5 pc. Future integral field spectrographs on 30m-class telescopes (e.g. HARMONI at E-ELT) are needed to make a difference, thanks to their spatial and spectral resolution.

Other young star clusters and star forming regions (e.g. $\rho{}$ Ophiuchi, Rigliaco et al. 2016; Cha I, Tsitali et al. 2015; NGC 1333, Foster et al. 2015; IC 348 Cottaar et al. 2015) show velocity gradients that might be interpreted as radial collapse of subvirial structures, or rotation, or a combination of both. The masses of most such clusters ($\sim{}100-500$ M⊙) are similar to that of the simulated star cluster in run D, where we do not see a clear signature of rotation.

To significantly increase the sample of ‘real-life’ star clusters we can compare our results with, we need to simulate more massive clusters, and/or to evolve them for (at least) several tens of Myr. This can be done only switching to a direct-summation N-body simulation, which treats stellar dynamics in the proper way (this point will be addressed in a follow-up study).

In this respect, our results might be useful to interpret the origin of rotation in the old globular clusters of the Milky Way (e.g. Bellazzini et al. 2012), but it must be kept in mind that our simulated star clusters are extremely different from globular clusters, as to mass, star formation history, and age. To understand whether gas torques can produce the rotation observed in globular clusters, we must focus on larger star cluster masses, and we must account for the extreme conditions of the early Universe (e.g. Trenti et al. 2015; Ricotti et al. 2016). Moreover, it is essential to trace the dynamical evolution of simulated star clusters for several Gyr, to constrain what fraction of the pristine rotation can survive to present time.

An alternative approach to compare our simulations with observational data is to look for rotation in the velocity field of dense gas in massive cluster-forming regions. For example, Liu et al. (2015) present Atacama Large Millimeter Array (ALMA) data of the rotating massive molecular clump G33.92+0.11, which contains two central massive ($\sim{}100-300$ M⊙) molecular cores. Liu et al. (2015) interpret the velocity field of G33.92+0.11 as possibly dominated by infall motions on large scale (\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}{}1 pc) and by rotational motion on the small scale (\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}{}0.3 pc). Their figure 6 is quite reminiscent of our Fig. 2. Other proto-clusters show marginally centrifugally supported accretion flows (e.g. Galván-Madrid et al. 2009; Baobab Liu et al. 2010; Cesaroni et al. 2011), suggesting that the central $<1$ pc regions of molecular clouds might be dominated by angular momentum. Future ALMA observations will provide invaluable information about rotation in molecular clouds and proto-star clusters.

5 Summary

We ran SPH simulations of turbulence-supported molecular clouds, and we studied the kinematics of the main cluster that forms in each cloud. In agreement with previous studies (e.g. Bonnell et al. 2003; Bonnell et al. 2011; Girichidis et al. 2011), we find that the star clusters assemble hierarchically, by the merger of several gaseous and stellar clumps. The gas component shows significant rotation in the region where the star cluster assembles (Fig. 2). This rotation feature arises from large-scale torques in the cloud and from angular momentum conservation.

The stellar component of the embedded star cluster inherits the rotation from the parent gas (Figs. 3–6). Rotation is apparent in the simulated star clusters (Fig. 7), if their final stellar mass is M_{\ast}{}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}{}1000 M⊙.

Our simulated star clusters form with large ellipticity ($\epsilon{}\sim{}0.70-0.75$ at $t\sim{}1.2$ $t_{\rm ff}$) and are markedly subvirial ($Q_{\rm vir}=0.2-0.3$ at $t\sim{}1.2$ $t_{\rm ff}$). During their evolution, the simulated star clusters tend to become rounder ($\epsilon{}\sim{}0.4-0.5$ at $t\sim{}2$ $t_{\rm ff}$) and to approach virial equilibrium ($Q_{\rm vir}=0.3-0.4$ at $t\sim{}2$ $t_{\rm ff}$), but their rotation signature is still apparent (Fig. 7). We find that rotation is more important than radial motions caused by subvirial collapse (Fig. 8). This result is a key test to probe the hierarchical formation scenario of star clusters, and might be useful to interpret the observed rotation signature in young massive star clusters (R136, Hénault-Brunet et al. 2012) and in old globular clusters (e.g. Bellazzini et al. 2012).

Acknowledgments

I thank the anonymous referee for their invaluable comments, which helped me improving the manuscript significantly. I warmly thank Francesco Palla for his invaluable suggestions. I also thank Germano Sacco, Elena Sabbi, Hauyu Baobab Liu, Giacomo Beccari, Angela Bragaglia, Eugenio Carretta, Enrico Vesperini, Anna Lisa Varri and Emanuele Ripamonti for useful discussions. The simulations were performed with the GALILEO cluster at CINECA (through CINECA Award N. HP10CRWH91) and with the SCIGHERA cluster at INAF-Osservatorio Astronomico di Padova. I acknowledge the CINECA Award N. HP10CRWH91 for the availability of high performance computing resources. I thank the authors of gasoline, especially J. Wadsley, T. Quinn and J. Stadel. I thank Tristan Hayfield for allowing me to use his implementation of cooling for gasoline. To analyze simulation outputs, we made use of the software tipsy666http://www-hpcc.astro.washington.edu/tools/tipsy/

tipsy.html and tipgrid777http://www.astrosim.net/code/doku.php?id=home:code:analysistools:misctools. MM acknowledges financial support from the Italian Ministry of Education, University and Research (MIUR) through grant FIRB 2012 RBFR12PM1F, from INAF through grant PRIN-2014-14 (Star formation and evolution in galactic nuclei), and from the MERAC Foundation.

Appendix A Impact of stochasticity on low-mass clouds

Stochastic fluctuations might affect our results, especially for the low-mass cloud (run D). To check the influence of stochasticity we ran four additional realizations of run D, i.e. four additional molecular clouds with $M=1.7\times{}10^{3}$ M⊙, $R=3$ pc, and $N_{\rm gas}=2\times{}10^{6}$, but with different random seeds for the turbulence. We call these four realizations runs D1, D2, D3, and D4.

Figure 9 is the same as Fig. 6 but for runs D1, D2, D3, and D4. The colour-coded map of Fig. 9 shows the component of the stellar velocity along the $x$ axis at $t=4$ Myr. With respect to the $yz$ plane, the velocity along the $x$ axis ($v_{\rm x}$), shown by the colour-coded map, can be interpreted as the line-of-sight velocity. Only run D3 shows some signature of rotation in Figure 9.

Fig. 10 is the same as Fig. 7 but for runs D1, D2, D3, and D4, i.e. it shows the rotation curve of the four test runs at $t=4$ Myr. Again, only the rotation curve of run D3 indicates that the star cluster rotates. The rotation curve of run D2 is more suggestive of a velocity gradient, whereas the rotation curves of runs D1 and D4 are very similar to that of run D (Fig. 7).

Table 3 (which is the same as Table 2 but for the test runs) gives some useful hints to interpret Figures 9 and 10. The total mass of a star cluster ($M_{\ast}$) varies considerably from one realization to the other, indicating that stochastic fluctuations are indeed significant. In particular, the star cluster formed in run D is also the less massive one among the five different realizations of the same cloud. In contrast, the star cluster in run D3 is by far the most massive one: its mass ($\sim{}800$ M⊙) is nearly three times larger than that of the star cluster in run D. This result confirms that the amount of natal rotation strongly correlates with the total mass of the star cluster. Moreover, this comparison also confirms that star clusters with final mass $<$ few $\times{}100$ M⊙ do not show clear indication of rotation. Rotation becomes significant only in star clusters with mass \lower 2.15277pt\hbox{;\buildrel>\over{\sim};}{}1000 M⊙.

Finally, we expect stochastic fluctuations to be less important in run C, because the number of particles (the mass of the cloud) is a factor of $\sim{}5$ ($\sim{}6$) larger than it is in run D. To check this, we ran a second realisation of run C, called run C1. Table 3 shows the comparison between runs C and C1. The mass of the main star cluster differs by only $\sim{}10$ per cent between run C and run C1. The other properties of runs C and C1 shown in Table 3 are also similar. The rotation signature of the main star cluster of run C1 is apparent from Fig. 11, consistent with what we found for run C.