© ESO, 2015

1. Introduction

The ongoing Gaia-ESO Survey (GES) at the VLT is a gold-mine of knowledge concerning the formation and evolution of young star clusters (YSCs) and OB associations. The aim of the GES, which began on December 31 2011 and will be completed in ~5 years, is to obtain high quality, uniformly calibrated spectra of >10⁵ stars in the Milky Way (Gilmore et al. 2012; Randich et al. 2013). The GES targets include a significant number of YSCs. The radial velocity (RV) and the chemical composition is being derived for a large number of YSC members with high accuracy. This information, when combined with data from the Gaia mission (Perryman et al. 2001), will enable reconstruction of the full three-dimensional spatial distribution and kinematics of several YSCs.

The first GES YSC target was the Gamma Velorum cluster, a marginally bound cluster of low-mass pre-main sequence (PMS) stars, surrounding the massive Wolf-Rayet (WR) – O binary system named $\mathit{\gamma }{\mathrm{2}}^{}$ Velorum (HD 68273, WR11; Smith 1968; Schaerer et al. 1997). The Wolf-Rayet and O star components of $\mathit{\gamma }{\mathrm{2}}^{}$ Velorum have current masses ~9 ± 2 and 30 ± 2 M_⊙, respectively (De Marco & Schmutz 1999), which implies that the initial masses were ~35 and 31 M_⊙, respectively (Eldridge 2009). $\mathit{\gamma }{\mathrm{2}}^{}$ Velorum is the most massive member of the common proper motion Vela OB2 association, composed of 93 early-type stars, spread over 100 square degrees (de Zeeuw et al. 1999). The Gamma Velorum cluster was firstly identified by Pozzo et al. (2000), on the basis of the strong X-ray emission of the PMS stars, and then further investigated by Jeffries et al. (2009). The distance of the Gamma Velorum cluster is about 350 pc (Pozzo et al. 2000; van Leeuwen 2007; Millour et al. 2007; North et al. 2007). The PMS stars in the Gamma Velorum cluster have an estimated age ~5−10 Myr (or even more), which is only marginally consistent with the age of the $\mathit{\gamma }{\mathrm{2}}^{}$ Velorum binary (~5.5 Myr, Eldridge 2009). Jeffries et al. (2014) suggest, on the basis of the Lithium (Li) depletion observed among its M-dwarf members, that the Gamma Velorum cluster is older than 10 Myr.

The GES data (Jeffries et al. 2014) show that the Gamma Velorum cluster has a complex kinematic structure: it is composed of two kinematically distinct populations (named population A and B by Jeffries et al. 2014). Population A (hereafter pop. A) has a 2.15 ± 0.48 km s^-1 lower RV than population B (hereafter pop. B). The velocity dispersion of the two populations is also different: pop. A (pop. B) has a best-fitting velocity dispersion $\mathit{\sigma }{\mathrm{A}}_{}\mathrm{=}\mathrm{0.34}\mathrm{\pm }\mathrm{0.16}$ km s^-1 ($\mathit{\sigma }{\mathrm{B}}_{}\mathrm{=}\mathrm{1.60}\mathrm{\pm }\mathrm{0.37}$ km s^-1).

The colour-magnitude diagram suggests that either pop. A is younger than pop. B by 0.4 ± 0.6 Myr (if the two populations are at the same distance), or pop. A and pop. B are coeval but pop. A is closer to us by 4 ± 5 pc. The spectral analysis indicates a possible age difference, but in the opposite sense: pop. A appears to be marginally more Li depleted than pop. B, which implies that pop. A might be 1−2 Myr older. The age difference derived from Li depletion, combined with the analysis of the colour-magnitude diagram, hints that pop. A might be closer to us than pop. B by a few parsecs.

Finally, the spatial distribution also suggests differences between pop. A and B: the former is slightly more concentrated than the latter (see Fig. 11 of Jeffries et al. 2014), and the centroids of the two populations are offset by 4.4 ± 3.0 arcmin (0.4 ± 0.3 pc for an assumed distance of 350 pc). The centroid of pop. A is almost coincident with the position of $\mathit{\gamma }{\mathrm{2}}^{}$ Velorum, indicating a connection between pop. A and the WR-O binary system.

Jeffries et al. (2014) propose three possible scenarios for the formation of the Gamma Velorum cluster. In the first scenario, pop. A and pop. B are the remnants of a larger cluster. After the evaporation of the parent molecular gas, only pop. A remained bound and in virial equilibrium, while pop. B became unbound. This scenario fails to explain the difference in RV between the two populations and the different Li depletion. In the second scenario, $\mathit{\gamma }{\mathrm{2}}^{}$ Velorum formed in a supervirial OB association. Pop. A represents the component of the initial OB association which was gravitationally captured by $\mathit{\gamma }{\mathrm{2}}^{}$ Velorum (e.g. Parker et al. 2014), while pop. B remained unbound. This second scenario explains the velocity differences, but cannot easily explain the age differences between the two populations. In the third scenario, favoured by Jeffries et al. (2014), pop. A formed in a denser region, while pop. B formed in a less dense region: pop. A is still marginally bound, while pop. B expands in the Vela OB2 association. This scenario could explain an age difference between the two populations, as well as a different spatial position and a different RV, because pop. A and pop. B formed in two slightly different star formation episodes.

In this paper, we focus on the third aforementioned scenario, investigating its implications in detail, by means of an N-body model of the Gamma Velorum cluster. The main idea is the following: simulations of turbulent molecular clouds suggest that star clusters form from the merger of sub-clusters (Bonnell et al. 2003, 2011; Bate 2009; Girichidis et al. 2011; Kruijssen et al. 2012). Furthermore, observations of star clusters embedded in molecular clouds indicate that the process of star formation is highly sub-structured (e.g. Cartwright & Whitworth 2004; Gutermuth et al. 2009; Sánchez & Alfaro 2009; André et al. 2010; De Marchi et al. 2013; Gouliermis et al. 2014; Kuhn et al. 2014; Kirk et al. 2014). The merger of sub-clusters can also explain some intriguing dynamical properties of YSCs (e.g. a fast mass segregation, Parker et al. 2011; Fujii et al. 2012; Parker et al. 2014). Thus, we propose that the two populations of the Gamma Velorum cluster originate from two different sub-clusters, born in the same parent molecular cloud. The velocity shift between two sub-clusters is expected to be of the same order of magnitude as the turbulent velocity in the parent molecular cloud. The observed turbulent motions of Galactic molecular clouds are of the order of ~1−2 km s^-1 (e.g. McKee & Ostriker 2007, and references therein), significantly larger than the sound speed of cold gas (~0.2 km s^-1). Furthermore, stars in the Orion Nebula Cluster (Tobin et al. 2009) and molecular gas in star-forming regions such as Perseus and Ophiuchus (André et al. 2007; Kirk et al. 2007; Rosolowsky et al. 2008) exhibit large-scale linear velocity gradients (≲ 1 km s^-1 pc^-1). Velocity gradients (~0.2−2 km s^-1 pc^-1) have also been found in simulations of molecular clouds and YSC formation (e.g. Offner et al. 2009). Such velocity offsets would naturally explain the observed RV difference of 2 km s^-1 between pop. A and pop. B. We investigate this possibility by means of direct-summation N-body simulations.

This paper is organized as follows. In Sect. 2, we discuss the simulation method. In Sect. 3, we present our results. In Sect. 4, we discuss the main advantages and drawbacks of our simulation method, and we propose future simulations and observational tests.

2. Methods and simulations

In this paper, we simulate the collision between two different sub-clusters, by means of direct-summation N-body simulations. The simulations were performed with the starlab public software environment (Portegies Zwart et al. 2001), in the version modified by Mapelli et al. (2013; see also Mapelli & Bressan 2013). starlab integrates the dynamical evolution of a star cluster, by means of a predictor-corrector Hermite scheme (implemented in the kira routine), and follows the stellar and binary evolution of the star cluster members (implemented in seba).

Initially, the two sub-clusters are described as two Plummer spheres (Plummer 1911). The single stars and the primary members of binaries were randomly generated following a Kroupa initial mass function (IMF, Kroupa 2001) with minimum and maximum mass 0.1 and 150 M_⊙, respectively. The secondary members of the binaries were randomly generated from a uniform distribution between 0.1 M_⊙ and m₁ (where m₁ is the mass of the primary member). We enable stellar and binary evolution at solar metallicity (consistent with Spina et al. 2014), as described in Mapelli & Bressan (2013). In the following, we call pop. A and pop. B the population of the first and the second simulated YSCs, respectively.

We ran several simulations with different initial conditions. The most relevant runs and the corresponding parameters are summarized in Table 1. The most important initial parameters appear to be (i) the number of particles in each star cluster (N_A and N_B, respectively) and the corresponding initial masses (M_A and M_B, respectively); (ii) the Plummer scale radius of each star cluster (r_A and r_B, respectively); (iii) the relative velocity (Δv) and the initial distance (D) between the centres-of-mass of the two clusters; (iv) the binary fraction f_bin; (v) the deviation of the star cluster from the virial equilibrium. The last condition is particularly important, because several observed YSCs show indication of non-virial velocities (e.g. Cottaar et al. 2012). To quantify the deviation from virial equilibrium, we define the Q parameter as Q ≡ K/ |W| (where K and W are the kinetic and the potential energy of the system, respectively). For a system in virial equilibrium Q = 0.5. Systems with Q> 0.5 (i.e. supervirial systems) might be the outcome of the evaporation of the parent gas. In fact, the gas component is not included in our models, but simulating a stellar population with initial velocities above the virial value mimics the fact that the potential well changed very fast in the recent past, i.e. that the parent gas evaporated rapidly.

We test two extreme values for the initial mass: a total mass (pop. A+pop. B) of ~150 M_⊙ (runs 12 and 13) and of ~800 M_⊙ (the other runs). The mass of the clusters simulated in runs 12 and 13 is similar to the current total mass of the Gamma Velorum cluster (where for current total mass we mean the mass within the 0.9 deg² area surveyed by GES), while a total mass of ~800 M_⊙ is likely closer to the initial mass of the system. The last column of Table 1 shows the most massive primordial binary system in each run. A binary as massive as $\mathit{\gamma }{\mathrm{2}}^{}$ Velorum is likely to form only in the most massive systems amongst those simulated here (~800 M_⊙). The large scatter in the values of M_A and M_B listed in Table 1 (given a fixed value of N_A and N_B, respectively) is due to stochastic fluctuations, because we randomly generate the masses of star particles according to a Kroupa IMF.

In run 21 the two sub-clusters do not collide. They are initially located at D = 5 pc (Table 1), and their separation increases (negative Δv). This run has been performed to check whether the properties of the Gamma Velorum cluster are consistent with a mere superposition between two sub-clusters (without physical interaction).

Fig. 1 Position of the simulated stars in the xy plane at t = 0, 2.7, 4.8 and 6.9 Myr in run 1. Red crosses: cluster 1 (corresponding to population A), black open circles: cluster 2 (corresponding to population B). The centre of the frame is the centre of mass of cluster 1 (pop. A).

Fig. 2 From top to bottom and from left to right: distribution of the RVs of stars in run 1, 2, 4, 5, 7, 11, 13, 15 and 16 at t = 4.8 Myr (see Table 1). The simulated RVs include randomly generated observational uncertainties. Solid red line: pop. A (cluster 1); dotted black line: pop. B (cluster 2); dashed blue line: sum of the two components. We define RV as the velocity along the y axis (i.e. we assume that the line-of-sight coincides with the direction that maximizes the velocity shift between the two clusters). Green shaded histogram: observed RV distribution (the same as in Fig. 6 of Jeffries et al. 2014).

3. Results

We simulate the collision between two different YSCs, during and after the first approach. Figure 1 shows the time evolution of run 1 in the xy plane. In all runs, the initial velocity shift is entirely along the y axis. In run 1, the two clusters are initially at a distance of 5 pc along the y axis (which is less than the typical radius of a giant molecular cloud in the Milky Way, e.g. Dame et al. 2001). They collide at time t ~2.5 Myr. The YSCs simulated in the other runs have similar evolution, but the collision time is different, since it depends on the relative initial distance, velocity and masses.

Detailed information about all runs (including those shown in Figs. 1 and 2) can be found in Table 2. In calculating the RVs in Table 2 (as well as in Fig. 2), we have assumed that the line-of-sight is along the y axis. Thus, the RV is the velocity along the y axis, i.e. the axis along which we generated the velocity shift Δv between the centre-of-mass of the two clusters¹. Therefore, the plane of the sky is the xz plane of the simulations.

In Table 2, we consider as reference times t ~2.5, 5 and 7 Myr since the beginning of the simulations. We stress that we select these three reference snapshots, since the radial distribution of the simulated stars, and especially the shift along the line-of-sight between pop. A and pop. B, are not consistent with the observations if t ≪ 2.5 Myr or t ≫ 7.5 Myr (see Sect. 3.2). On the other hand, the RV distribution does not change significantly at different times (see Table 2), provided that t ≥ 2.5 Myr (i.e. provided that the two sub-clusters already collided). Thus, the main constraints on the time elapsed since the collision come from the radial distribution, while the RV distribution is less affected. Finally, we recall that the time indicated in Table 2 does not necessarily correspond to the age of the stars in pop. A and pop. B: the tabulated time is the time elapsed since the beginning of the simulations, but the two sub-clusters might have formed before t = 0. Evaluating the age of the Gamma Velorum cluster, which is rather uncertain (see the Introduction), is beyond the aims of this paper.

3.1. Radial velocity and proper motions

Figure 2 shows the distribution of the RVs for nine selected runs at t = 4.8 Myr (for other runs, see Table 2). The simulated RVs shown in Fig. 2 include observational uncertainties. These were randomly drawn from the distribution reported in Fig. 2 of Jeffries et al. (2014). The zero-point of the simulated RVs is the RV velocity of the centre-of-mass of the total system (pop. A+pop. B stars). In Fig. 2, we also show the observed RV distribution (the same as in Fig. 6 of Jeffries et al. 2014). The observed RV distribution has been shifted to match the simulated ones.

In our simulations, we exactly know which stars belong to the first and to the second cluster, respectively, and we can derive their intrinsic velocity dispersions². Columns 7 and 8 of Table 2 show the intrinsic velocity dispersions (${\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\sigma }\end{array}}_{\mathrm{A}}$ and ${\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{\sigma }\end{array}}_{\mathrm{B}}$) of the simulated pop. A and pop. B, along the line-of-sight. Column 9 of Table 2 ($\begin{array}{c}{\mathrm{˜}}_{}\\ \mathrm{\Delta }\end{array}$RV) shows the intrinsic velocity shift of the two simulated populations along the line-of-sight, i.e. the difference between the mean velocity of pop. A and that of pop. B along the line-of-sight.

On the other hand, we want to fit the total RV distribution of each simulation in the same way as Jeffries et al. (2014), for a better comparison with the data of the Gamma Velorum cluster. Thus, we apply a fit with two Gaussian components to the total pop. A+pop. B distribution, for all runs. The fitting procedure is carried out on binned data, and adopting a procedure equivalent to the one described in Sect. 4 of Jeffries et al. (2014, i.e. statistically correcting for the binary fraction and for the observational uncertainties a posteriori). From this fitting procedure, we obtain the velocity dispersions of the two simulated populations along the line-of-sight ($\mathit{\sigma }{\mathrm{A}}_{}$ and $\mathit{\sigma }{\mathrm{B}}_{}$), the velocity shift of the two simulated populations along the line-of-sight (ΔRV), and the fraction of estimated pop. A members with respect to the total number of simulated stars (f_A). The results are shown in Cols. 3−6 of Table 2. The intrinsic velocity dispersions generally agree with the ones derived from this fitting procedure, but they might also disagree by a factor of ~2 (especially for pop. B), when the distribution of simulated RVs deviates too much from a Gaussian distribution, and when it is more difficult to disentangle pop. A from pop. B.

We adopt two criteria to check whether a simulation is consistent with the observations of the Gamma Velorum cluster. Namely, we decide that a simulation (at a given time) is not consistent with the observations if (i) at least one of the six reference quantities ($\mathit{\sigma }{\mathrm{A}}_{}$, $\mathit{\sigma }{\mathrm{B}}_{}$, ΔRV, f_A, the distance between the centre of mass of pop. A and pop. B along the line-of-sight Δy, and the distance between the centre of mass of pop. A and pop. B in the plane of the sky Δxz) differs from the observed one by more than 3σ; or (ii) P_KS< 0.05 (where P_KS is the probability that the chance deviation between the observed RV distribution and the simulated RV distribution is expected to be larger, according to a Kolmogorov-Smirnov test, last column of Table 2).

Run 4 (where Q_A = 3) and run 7 (where Q_B = 12.5) do not satisfy the first criterion (because of the large values of $\mathit{\sigma }{\mathrm{A}}_{}$ and $\mathit{\sigma }{\mathrm{B}}_{}$, respectively). While P_KS is < 0.05 in runs 11 (Δv = 3.5 km s^-1) and 16 (f_bin = 0.75). From Fig. 2 and Table 2, it is apparent that some of our models match the kinematics of the Gamma Velorum cluster quite well. Furthermore, we are able to put constraints on the most relevant initial parameters. In particular, if the initial velocity shift is Δv> 3 km s^-1 (runs 10 and 11), then ΔRV is larger than the observed one by more than 1σ, for the entire simulation. The tidal forces exerted during the interaction between pop. A and pop. B are not sufficient to quench this velocity significantly. If Δv ≥ 3.5 km s^-1, P_KS is very low (< 0.05).

More interestingly, if the simulated pop. B is close to virial equilibrium (Q_B ≤ 2, runs 5 and 6), then its velocity dispersion (≤ 0.7 km s^-1) is much smaller than the observed one (1.60 ± 0.37 km s^-1), and the RV distributions of the two populations appear almost completely separated (Fig. 2). On the other hand, if Q_B ≥ 12.5, then the simulated velocity dispersion of pop. B is too large with respect to the observed one ($\mathit{\sigma }{\mathrm{B}}_{}\mathrm{~}\mathrm{2.8}$ km s^-1 for Q_B = 12.5 in run 7, i.e. 3σ higher than the observed value). Thus, the model can reproduce the observed RV distribution only if pop. B is significantly hotter than a virial system, but not too hot (2 <Q_B< 12.5). This implies that pop. B is unbound, likely as a consequence of the evaporation of the parent molecular gas.

Other constraints can be put on the virial ratio of pop. A: if M_A ~ 500 M_⊙ and Q_A ≥ 3.0, then the velocity dispersion of the simulated pop. A (~1 km s^-1, run 4) is too large to be consistent with the observed one (0.34 ± 0.16 km s^-1). In contrast, values of Q_A< 2.0 are consistent with the data, indicating that pop. A is approximately in virial equilibrium.

The binary fraction f_bin is very important for the comparison with observations. In fact, unresolved binaries tend to broaden the RV distribution (we recall that in Fig. 2 we plot even the simulated binaries as unresolved, for a better comparison with the data). In most simulations, we assume f_bin = 0.46 for analogy with the results of Jeffries et al. (2014), but we run two test cases with f_bin = 0 (run 15) and 0.75 (run 16), respectively. The effect of f_bin is particularly important for pop. A. If f_bin = 0 (run 15), then the scatter in the RVs of pop. A is much smaller: a larger value of Q_A is requested to reproduce $\mathit{\sigma }{\mathrm{A}}_{}$ if f_bin = 0 (run 15), but the simulations still fail to reproduce the high-velocity tail of the RV distribution. In contrast, if f_bin = 0.75 (run 16), the RV distribution is broader, and a lower value of Q_A is requested to match the data.

The characteristic radii of the two populations are more important for the spatial distribution than for the RV. On the other hand, the characteristic radii have also an impact on the RV, and there is partial degeneracy with the effect of the virial ratio. For example, if pop. B is initially spread over a larger radius (r_B = 2 pc, run 18), for the same initial mass M_B, then its velocity dispersion will be lower. Thus, the observations can be matched only for a larger value of Q_B, but this affects the disruption timescale of pop. B.

Finally, the mass of the two populations is another important quantity, and there is partial degeneracy with the virial ratio. For example, run 12 has the same parameters as run 1, but the mass of the two clusters is lower by a factor of ~5. The RV distributions of pop A. and pop. B partially overlap in run 1, while they are well separated in run 12. The result is that run 12 does not match the data (especially $\mathit{\sigma }{\mathrm{B}}_{}$ and ΔRV, Table 2). In case of a lower initial total mass, larger values of Q_A and Q_B are requested, to obtain a best matching with the data. For example, in our simulations, if M_A + M_B ~ 150 M_⊙ instead of ~800 M_⊙, we need Q_A ≥ 4 and Q_B ≥ 8 (as in run 13, Fig. 2) to reproduce the Gamma Velorum data.

The values adopted in run 13 match the current mass of the Gamma Velorum cluster quite well, but present an issue: they can hardly account for the formation of the $\mathit{\gamma }{\mathrm{2}}^{}$ Velorum WR-O binary. In Table 1, we report the masses of the members of the most massive binary in each run. The most massive binaries in runs 12 and 13 are well below 10 M_⊙. Even if we account for stochastic fluctuations, the formation of a ~(30, 30) M_⊙ binary in clusters with initial mass ~100 M_⊙ is very unlikely. Thus, we can argue either that the formation of $\mathit{\gamma }{\mathrm{2}}^{}$ Velorum is not related to the rest of the Gamma Velorum cluster (which is quite unlikely), or that the initial mass of the clusters was a factor of ≳ 5 higher, and then the partial cluster dissolution due to gas evaporation has led to a smaller total mass. Our runs with virial ratio >0.5 account for this effect of gas evaporation and ongoing cluster dissolution. Furthermore, our simulated YSCs evolve dynamically and progressively spread over an area larger than the one observed by the GES. This partially explains the difference between the mass estimated by the observations (within the 0.9 deg² area surveyed by GES) and the total mass of the simulated YSCs in run 1 (see the discussion in Sect. 3.2).

Fig. 3 Proper motions in the plane of the sky for run 1 at t = 4.8 Myr. In the main panel, red crosses: pop. A, black open circles: pop. B. In the marginal histograms, red (black) histogram: pop. A (pop. B). On the y axis of the marginal histograms: number of stars per bin, normalized to the total number of stars in each cluster (N∗/Ntot). The simulated proper motions shown in this figure do not include observational uncertainties.

Figure 3 shows the proper motions in the plane of the sky (the xz plane, according to our convention) for run 1. Pop. A and pop. B are not significantly decoupled, if we consider the proper motions. On the other hand, since pop. B is more supervirial than pop. A, the distribution of proper motions has a larger dispersion for pop. B ($\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{\sigma }}_{\mathit{x,}\hspace{0.17em}\mathrm{B}}\end{array}\mathrm{=}\mathrm{0.96}$ km s^-1, $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{\sigma }}_{\mathit{z,}\hspace{0.17em}\mathrm{B}}\end{array}\mathrm{=}\mathrm{1.00}$ km s^-1, where $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{\sigma }}_{\mathit{x,}\hspace{0.17em}\mathrm{B}}\end{array}$ and $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{\sigma }}_{\mathit{z,}\hspace{0.17em}\mathrm{B}}\end{array}$ are the intrinsic velocity dispersions of pop. B along the x and the z axis, respectively) than for pop. A ($\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{\sigma }}_{\mathit{x,}\hspace{0.17em}\mathrm{A}}\end{array}\mathrm{=}\mathrm{0.46}$ km s^-1, $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{\sigma }}_{\mathit{z,}\hspace{0.17em}\mathrm{A}}\end{array}\mathrm{=}\mathrm{0.47}$ km s^-1, where $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{\sigma }}_{\mathit{x,}\hspace{0.17em}\mathrm{A}}\end{array}$ and $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{\sigma }}_{\mathit{z,}\hspace{0.17em}\mathrm{A}}\end{array}$ are the intrinsic velocity dispersions of pop. A along the x and the z axis, respectively). Furthermore, if the plane of the sky is slightly different from the assumed one (i.e. if the relative velocity vector between the two YSCs is not perfectly aligned to the line-of-sight), we expect to observe a further difference between pop. A and pop. B. Gaia data will be crucial to test this feature.

3.2. Spatial distribution

Fig. 4 Normalized radial stellar density in run 1 (left-hand panels) and in run 13 (right-hand panels) at different times. The simulated cluster is assumed to be at a distance of 350 pc. Solid red line: pop. A; dotted black line: pop. B. The density n is calculated as stars per arcmin2, within concentric annuli. The normalization n0 was chosen so that the total area below the histogram of pop. A is one. The error bars account for Poisson uncertainties.

In the data published by Jeffries et al. (2014), pop. A appears to be marginally more concentrated than pop. B, and the two centroids are offset by ~4.4 ± 3.0 arcmin (~0.45 ± 0.31 pc). Our simulations can naturally explain the radial offset and the different centroid of the two populations of the Gamma Velorum cluster: in our models, pop. A and pop. B come from two different sub-clusters, which might be born with different concentration and are expected to have different centroids.

Figure 4 shows the normalized radial density of stars in pop. A and pop. B in runs 1 and 13, at different times. Pop. B is slightly less concentrated than pop. A. Furthermore, the concentration of pop. B diminishes with time, as a result of the fact that the second (less massive) cluster is being tidally disrupted by the first (more massive) cluster (see also Fig. 1). The probability that the chance deviation between the spatial distribution of pop. A and that of pop. B is expected to be larger, according to a Kolmogorov-Smirnov test, is P_KS< 10^-5 and < 0.1 for run 1 and run 13, respectively. Thus, the two simulated populations follow different spatial distributions. In the other runs, we observe the same trend, even if the relative difference between pop. A and pop. B changes (depending on the initial conditions). The simulated spatial distributions (Fig. 4) are qualitatively similar to the observed radial spatial distribution of the Gamma Velorum cluster, as shown in Fig. 11 of Jeffries et al. (2014), even if we recall that these two figures cannot be directly compared³.

Column 11 of Table 2 shows that there is an offset of 0−1.5 pc between the centres-of-mass of the two simulated clusters in the xz plane (i.e. in the assumed plane of the sky) at 0 ≤ t/ Myr ≤ 8, depending on the initial conditions. For most runs, the offset is consistent with the observed one (~0.45 ± 0.31 pc).

Finally, Table 2 also shows that there is an offset of 0−18 pc between the centres-of-mass of the two simulated clusters along the y axis (i.e. along the line-of-sight) at 0 ≤ t/ Myr ≤ 8. The colour-magnitude diagram of the Gamma Velorum cluster is consistent with an offset of 4 ± 5 pc (along the line-of-sight), if we assume that pop. A and pop. B are coeval (or with an even larger offset, if we assume that pop. A is slightly older, as indicated by the Li depletion, Jeffries et al. 2014). Thus, the simulated offset along the line-of-sight is fairly consistent with the one indicated by the observations.

Fig. 5 Position of simulated stars belonging to pop. A (red crosses) and pop. B (black open circles), projected in the plane of the sky (the xz plane of the simulation), at three different times (t = 2.5, 4.8 and 6.9 Myr since the beginning of the simulation) in run 1. The centre of the frame is the centre of mass of the entire cluster (pop. A + pop. B). We note that pop. B expands much faster than pop. A given its supervirial condition. Stars belonging to pop. B can be found at >10 pc from the centre of pop. A.

Figure 5 shows the position of simulated stars belonging to pop. A (red crosses) and pop. B (black circles), projected in the plane of the sky (the xz plane of the simulation), at three different times (t = 2.5, 4.8 and 6.9 Myr since the beginning of the simulation) in run 1. It is apparent that pop. B expands faster than pop. A, because of its supervirial condition. We note that stars of pop. B can be found at >10 pc from the centre of pop. A. Recently, Sacco et al. (2015) found that 15 stars, in the direction of the NGC 2547 cluster, i.e. about two degrees (~10 pc) south of $\mathit{\gamma }{\mathrm{2}}^{}$ Velorum, have the same properties as the pop. B members of the Gamma Velorum cluster. To compare our simulations with the results of Sacco et al. (2015), we select a box in our run 1 (at t = 4.8 Myr) with approximately the same size and the same location (with respect to the centre of the cluster) as the area analyzed by Sacco et al. (2015). We count 16 pop. B objects (6 of which are binaries) in this box. This is in excellent agreement with the result of Sacco et al. (2015), who find 15 pop. B members in the field of NGC 2547. This result supports our scenario, and strengthens the evidence that pop. B is strongly supervirial. Figure 5 also shows that our simulated star clusters are more extended than the area reported by Jeffries et al. (2014, ~0.9 deg²). It is reasonable to expect that the outer rim of the Gamma Velorum cluster is much larger than the area observed by Jeffries et al. (2014).

3.3. Collision or line-of-sight superposition?

Our models are able to reproduce most properties of the Gamma Velorum cluster, provided that pop. B is significantly supervirial. However, we may wonder whether we really need a physical collision between the two sub-clusters. Is it possible to explain the properties of the Gamma Velorum cluster with a mere line-of-sight superposition between pop. A and pop. B, without requiring any interactions between the two sub-clusters?

To answer this question, we have performed run 21, in which the two sub-clusters form close to each other (5 pc) but increase their separation rather than collide. The initial properties of run 21 are the same as those of run 1 (which matches the properties of the Gamma Velorum cluster particularly well), apart from the sign of the relative velocity. The main kinematic properties of run 21 (i.e. the values of $\mathit{\sigma }{\mathrm{A}}_{}$, $\mathit{\sigma }{\mathrm{B}}_{}$, and ΔRV) are very similar to those of run 1, indicating that the velocity field of the two populations is not severely affected by the collision. The main reason is that the relative velocity of the two sub-clusters is large with respect to their internal kinematics. On the other hand, the line-of-sight distance between the two sub-clusters becomes too large (>19 pc, 3σ larger than the observed value) at t ≳ 7.1 Myr.

The same conclusion can be reached on the basis of a back-of-the-envelope calculation. The GES data indicate that there is an offset ΔRV = 2.15 ± 0.48 km s^-1 between the RV of pop. A and pop. B. Furthermore, pop. A might be closer to us by 4 ± 5 pc. If we assume that pop. A and pop. B formed in the same place and then drifted away from each other, pop. B is now at a distance d_AB from pop. A: ${\mathit{d}}_{\mathrm{AB}}\mathrm{=}\mathrm{10}\hspace{0.17em}\mathrm{pc}\hspace{0.17em}\left(\frac{\mathrm{\Delta }\mathit{RV}}{\mathrm{2}\hspace{0.17em}\mathrm{km}\hspace{0.17em}{\mathrm{s}}^{-1}}\right)\hspace{0.17em}\left(\frac{{\mathit{t}}_{\mathrm{B}}}{\mathrm{5}\hspace{0.17em}\mathrm{Myr}}\right)\mathit{,}$(1)where t_B is the age of pop. B. Thus, even if pop. A and pop. B formed in the same place 5 Myr ago, d_AB is already 1 σ larger than the observed displacement between pop. A and pop. B. If pop. B formed farther from us than pop. A, then we would expect a larger observed displacement between pop. A and pop. B. On the other hand, if pop. A and pop. B formed in the same place, they should have interacted between each other before drifting away. For the same reason, if pop. B formed closer to us, it should have interacted with pop. A before drifting away. This favours the scenario of a gravitational encounter between the two sub-clusters. This argument is weakened by the fact that our measurement of the line-of-sight displacement between pop. A and pop. B relies on the colour-magnitude diagram and is quite inaccurate. Parallax measurements by Gaia will give invaluable hints on this point.

Do the simulations show any other quantitative difference between a simple line-of-sight superposition and a gravitational interaction between pop. A and pop. B? A physical collision leaves an imprint on the structure of the two sub-clusters, and especially of the lighter one. In our simulations, pop. B is tidally perturbed and stretched by the gravitational encounter with the more massive pop. A. The tidal perturbation can be quantified by comparing the half-mass radius projected in the plane where the tidal force is maximum (r_h(yz) or r_h(xy)), with the half-mass radius projected in the plane where the tidal force is minimum (r_h(xz)). Figure 6 shows that r_h(yz) is significantly larger than r_hm(xz), indicating that pop. B is severely stretched by the interaction. Again, this might be checked with forthcoming parallax measurements by Gaia.

In summary, a collision scenario and a mere line-of-sight superposition show no significant differences in the velocity dispersions and in the RV distribution. The only way to distinguish between these two scenarios is to measure the line-of-sight displacement between the two populations and any possible tidal deformation of pop. B.

Fig. 6 Half-mass radius of pop. B as a function of time in run 1. rh(xz): half-mass radius projected in the xz plane. rh(yz): half-mass radius projected in the yz plane (where the tidal force exerted by pop. A is maximum).

3.4. Relaxing the assumption of line-of-sight superposition

It is very unlikely that the two sub-structures of a YSC are observed exactly along the line-of-sight. How much can we relax the assumption of line-of-sight superimposition and still match the observed data? To check this hypothesis, we focus on our fiducial run (run 1), and change the line-of-sight by rotating the yz plane about the x-axis by an angle θ. We stop rotating when at least one of our reference quantities for the comparison with the data ($\mathit{\sigma }{\mathrm{A}}_{}$, $\mathit{\sigma }{\mathrm{B}}_{}$, ΔRV, f_A, Δy and Δxz) differs from the best observed value by more than 3σ. The first parameter for which the discrepancy becomes >3 σ is the difference between the centroids (Δxz). At an angle θ = 12° between the original y axis and the new y axis (y′), we obtain $\mathit{\sigma }{\mathrm{A}}_{}\mathrm{=}\mathrm{0.48}$ km s^-1, $\mathit{\sigma }{\mathrm{B}}_{}\mathrm{=}\mathrm{1.64}$ km s^-1, ΔRV = 1.99 km s^-1, f_A = 0.40, Δy′ = 7.07 pc and Δxz′ = 1.50 pc (the best value for the observations is Δxz = 0.45 ± 0.31 pc, see Table 2). Thus, we can conclude that configurations in which the two sub-clusters are superimposed or slightly offset (by an angle θ ≤ 12°) reasonably match the data. Larger angles are possible, but only for different collision times. Figure 7 shows the RV distribution (along y′), the distribution of proper motions along z′, and the projected positions (in the xy′ and xz′ planes) if run 1 is rotated by θ = 12° about the x axis.

The probability of observing two sub-clusters offset by an angle θ ≤ 12° is ~0.022. On the other hand, we recall that while the centroid of pop. A is quite well constrained, the centroid of pop. B is more uncertain, especially if pop. B is much more extended (as suggested by Sacco et al. 2015). The uncertainty on the centroid of pop. B might further relax the request of partial line-of-sight superposition.

Fig. 7 Properties of run 1 (at t = 4.8 Myr) if rotated by an angle θ = 12° about the x axis. Top left-hand panel: RV distribution (symbols are the same as in Fig. 2); top right-hand panel: distribution of proper motions along z′ (symbols are the same as in the marginal histograms of Fig. 3, in particular pop. A and pop. B correspond to the red and black histogram, respectively); bottom left-hand panel: positions in the xy′ plane (symbols are the same as in Fig. 1, in particular pop. A and pop. B correspond to red crosses and black open circles, respectively); bottom right-hand panel: positions in the xz′ plane (symbols are the same as in Fig. 5, in particular pop. A and pop. B correspond to red crosses and black open circles, respectively).

4. Conclusions

Several hints indicate that YSCs form from the merger of sub-clusters, born from the same molecular cloud. In this paper, we have simulated the collision between two sub-clusters. We have shown that the collision product reproduces some interesting features of the Gamma Velorum cluster, recently observed by the GES. The GES data (Jeffries et al. 2014) indicate that the Gamma Velorum cluster is composed of two kinematically decoupled populations: pop. A and pop. B. The two populations have a marginally different radial concentration and slightly offset centroids (by ~0.5 pc).

Our simulations can naturally explain the RV offset between the two populations of the Gamma Velorum cluster, as well as their different intrinsic velocity dispersions (see Fig. 2). Our simulations also account for the different concentration and the different centroid of the two populations (Figs. 4 and 5). The GES data suggest that pop. A is older by 1−2 Myr, based on the Li depletion. Our simulations are consistent with a small age difference between the two populations. In the scenario of a collision between two sub-clusters, the two populations can either be perfectly coeval or have formed in two slightly different star formation episodes. The latter hypothesis is supported by observations of similar age differences in other star forming regions (even if the question of age spread in star forming regions is still debated, see Jeffries et al. 2011, and references therein).

We predict that the dispersion of the distribution of proper motions is broader for pop. B than for pop. A, even if the relative velocity vector between the two sub-clusters is aligned with the line-of-sight, as a result of the fact that pop. B is supervirial (Fig. 3). Furthermore, if a component of the relative velocity vector is normal to the line-of-sight, we expect a difference in the average proper motions of pop. A and pop. B (e.g. Fig. 7). This will soon be tested with Gaia data.

Furthermore, our simulations can help us reducing the allowed parameter space for the initial conditions of the Gamma Velorum cluster. For example, models with a low binary fraction (f_bin ~ 0) can hardly account for the observed RV data. Similarly, models in which the initial relative velocity between the clusters is ΔRV > 3 km s^-1 are in disagreement with the observed RVs. Star cluster models with an initial total mass ≤ 150 M_⊙ can reproduce most of the observed features, but cannot account for the formation of the massive $\mathit{\gamma }{\mathrm{2}}^{}$ Velorum binary system. In contrast, models with a larger mass (≥ 500 M_⊙) can host massive binary systems, similar to the $\mathit{\gamma }{\mathrm{2}}^{}$ Velorum binary.

Finally, the simulations suggest that pop. B is not in virial equilibrium: supervirial models (with virial ratio 2 <Q_B< 12.5) are in better agreement with the GES data. The virial ratio Q is an essential ingredient of our models, and strongly affects the results: it is impossible to match the observed RV distribution of the Gamma Velorum cluster, without allowing Q_B to assume a value ≫ 0.5. This result might indicate that the gas from the parent molecular cloud evaporated very fast, leaving the stellar population out of virial equilibrium. Thus, pop. B is about to dissolve in the field. Furthermore, the physical meaning of Q is connected with the star formation efficiency in the cloud. A larger value of Q corresponds to a lower star formation efficiency, and consequently to a stronger effect of gas expulsion (e.g. Tutukov 1978; Hills 1980; Lada et al. 1984; Baumgardt & Kroupa 2007). Thus, we can infer that pop. A formed with higher star-formation efficiency, managed to survive gas evaporation, and preserved virial equilibrium, while pop. B was formed in a less efficient star formation episode, and is not going to survive gas evaporation. A possible scenario is that pop. B formed from a less dense molecular cloud core, where star formation efficiency was intrinsically lower. Another possible scenario is that pop. B formed later than pop. A (as indicated by the Li depletion, Jeffries et al. 2014), and that the gas surrounding pop. B was evaporated by the radiation field of pop. A, quenching the star formation episode of pop. B quite abruptly.

Since it is supervirial, pop. B expands faster than pop. A. Simulated stars belonging to pop. B lie at >10 pc from the centre of pop. A at t ≳ 4 Myr since the beginning of the simulation. Sacco et al. (2015) have recently found that some stars located in the NGC 2547 cluster have the same properties as those of pop. B. Since the NGC 2547 cluster lies about two degrees (~10 pc) south of $\mathit{\gamma }{\mathrm{2}}^{}$ Velorum, this result supports our scenario and strengthens the evidence that pop. B is strongly supervirial.

The main intrinsic limits of our simulations are the following: we have not included the parent molecular gas (the effects of gas evaporation are indirectly modelled by assuming a virial ratio >0.5), and we do not account for the Galactic tidal field directly. We will refine our simulations, by including the missing ingredients, in a forthcoming paper. Our model requires a fine-tuning: we need to assume that most of the RV offset between the two sub-clusters is aligned along the line-of-sight (± 12°; for larger angles, the distance between the centroids differs by > 3σ with respect to the measured value). If the direction of the relative velocity vector between the two sub-clusters is significantly different from the line-of-sight, then we would expect to see a larger spatial offset between the two populations. If the Gamma Velorum cluster represents the case of a fortuitous alignment between velocity offset and line-of-sight, then we would expect to observe a number of “twin” sub-clusters in young star formation complexes, showing a larger spatial offset than the Gamma Velorum cluster. The forthcoming data by the Gaia mission will likely shed light on this issue, by providing accurate proper motion and parallax measurements.

Our results show that the Gamma Velorum cluster is an ideal test-bed to check different scenarios for the formation of young star clusters in the Milky Way. The RV precision achieved by the GES allowed us to make an unprecedented comparison between the kinematics of simulated and observed clusters. Similar kinematic data of other young star clusters and associations in the GES sample will be essential to shed light on the formation of star clusters in the local Universe.

Acknowledgments

We thank the anonymous referee for the careful reading of the manuscript and for the helpful comments. Based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under programme ID 188.B-3002. These data products have been processed by the Cambridge Astronomy Survey Unit (CASU) at the Institute of Astronomy, University of Cambridge, and by the FLAMES/UVES reduction team at INAF/Osservatorio Astrofisico di Arcetri. These data have been obtained from the Gaia-ESO Survey Data Archive, prepared and hosted by the Wide Field Astronomy Unit, Institute for Astronomy, University of Edinburgh, which is funded by the UK Science and Technology Facilities Council. This work was partly supported by the European Union FP7 programme through ERC grant number 320360 and by the Leverhulme Trust through grant RPG-2012-541. We acknowledge the support from INAF and Ministero dell’ Istruzione, dell’ Università’ e della Ricerca (MIUR) in the form of the grant “Premiale VLT 2012”. The results presented here benefit from discussions held during the Gaia-ESO workshops and conferences supported by the ESF (European Science Foundation) through the GREAT Research Network Programme. M.M. acknowledges financial support from the Italian Ministry of Education, University and Research (MIUR) through grant FIRB 2012 RBFR12PM1F, from INAF through grant PRIN-2011-1 and from CONACyT through grant 169554. The authors acknowledge the CINECA Award N. HP10B3BJEW, HP10CLI3BX, HP10CXB7O8, HP10C894X7, HP10CGUBV0, HP10CP6XSO and HP10C3ANJY for the availability of high performance computing resources and support.

References

All Tables

All Figures

	Fig. 1 Position of the simulated stars in the xy plane at t = 0, 2.7, 4.8 and 6.9 Myr in run 1. Red crosses: cluster 1 (corresponding to population A), black open circles: cluster 2 (corresponding to population B). The centre of the frame is the centre of mass of cluster 1 (pop. A).
In the text

Fig. 2 From top to bottom and from left to right: distribution of the RVs of stars in run 1, 2, 4, 5, 7, 11, 13, 15 and 16 at t = 4.8 Myr (see Table 1). The simulated RVs include randomly generated observational uncertainties. Solid red line: pop. A (cluster 1); dotted black line: pop. B (cluster 2); dashed blue line: sum of the two components. We define RV as the velocity along the y axis (i.e. we assume that the line-of-sight coincides with the direction that maximizes the velocity shift between the two clusters). Green shaded histogram: observed RV distribution (the same as in Fig. 6 of Jeffries et al. 2014).

In the text

Fig. 3 Proper motions in the plane of the sky for run 1 at t = 4.8 Myr. In the main panel, red crosses: pop. A, black open circles: pop. B. In the marginal histograms, red (black) histogram: pop. A (pop. B). On the y axis of the marginal histograms: number of stars per bin, normalized to the total number of stars in each cluster (N∗/Ntot). The simulated proper motions shown in this figure do not include observational uncertainties.

In the text

Fig. 4 Normalized radial stellar density in run 1 (left-hand panels) and in run 13 (right-hand panels) at different times. The simulated cluster is assumed to be at a distance of 350 pc. Solid red line: pop. A; dotted black line: pop. B. The density n is calculated as stars per arcmin2, within concentric annuli. The normalization n0 was chosen so that the total area below the histogram of pop. A is one. The error bars account for Poisson uncertainties.

In the text

Fig. 5 Position of simulated stars belonging to pop. A (red crosses) and pop. B (black open circles), projected in the plane of the sky (the xz plane of the simulation), at three different times (t = 2.5, 4.8 and 6.9 Myr since the beginning of the simulation) in run 1. The centre of the frame is the centre of mass of the entire cluster (pop. A + pop. B). We note that pop. B expands much faster than pop. A given its supervirial condition. Stars belonging to pop. B can be found at >10 pc from the centre of pop. A.

In the text

	Fig. 6 Half-mass radius of pop. B as a function of time in run 1. rh(xz): half-mass radius projected in the xz plane. rh(yz): half-mass radius projected in the yz plane (where the tidal force exerted by pop. A is maximum).
In the text

Fig. 7 Properties of run 1 (at t = 4.8 Myr) if rotated by an angle θ = 12° about the x axis. Top left-hand panel: RV distribution (symbols are the same as in Fig. 2); top right-hand panel: distribution of proper motions along z′ (symbols are the same as in the marginal histograms of Fig. 3, in particular pop. A and pop. B correspond to the red and black histogram, respectively); bottom left-hand panel: positions in the xy′ plane (symbols are the same as in Fig. 1, in particular pop. A and pop. B correspond to red crosses and black open circles, respectively); bottom right-hand panel: positions in the xz′ plane (symbols are the same as in Fig. 5, in particular pop. A and pop. B correspond to red crosses and black open circles, respectively).

In the text