The formation and survival of young embedded star clusters

Most if not all stars form in clustered environments (see pictures to the right), ranging from low-density, low-number, loosely bound associations, to extreme star bursts forming hundreds of massive star clusters in confined areas (now called star cluster complexes). With the exception of the extreme star formation events, we see that after only 10-20 Myr the majority of stars are dispersed into the field. The destruction of the young star clusters has been dubbed 'infant mortality'. The best candidate to explain this early destruction is the expulsion of the remaining gas (i.e. gas which was not converted into stars), by stellar feedback (through radiation and stellar winds and finally once the first super-novae explode).

Giant molecular clouds (GMCs) turn only a few to a few tens of per cent of their gas into stars, referred to as the Star Formation Efficiency (SFE - the ratio between mass converted into stars and total mass of the star forming region). Therefore, the potential of the embedded star cluster (star cluster which is still within its gas-cloud) is dominated by the remaining gas. The removal of this gas by feedback will leave the stars of the cluster in an out-of-equilibrium state. In most theoretical studies, the embedded star cluster is assumed to be in virial equilibrium and forms a spherically symmetric distribution of stars. In this case all studies agree that, to obtain a bound remnant star cluster, the SFE has to be extremely high - of order 30 per cent.

But, recent observations of embedded star forming regions show that stars do not form in spherically symmetric distributions, nor in virial equilibrium. Instead they form in clumps and filaments, and most likely dynamically cool.

Star formation in filamentary structures are studied by various groups using state-of-the-art SPH codes on very fast super-computers (see example to the right). These simulations are very time-consuming (many months/years of CPU-time) on the one hand and therefore cannot provide a statistical significant parameter search. On the other hand, these simulations lack the exact collisional stellar dynamical treatment of the newly formed stars. Furthermore, it is almost impossible to isolate the influence one particular parameter has on the outcome of the whole simulations.

In this project we try to assess the influence of each ingredient to the outcome separately by using very fast highly idealised models. For this we make use of the state-of-the-art direct N-body code N-body 6. This code provides us with the exact stellar dynamical treatment of the young stars and allows us to perform dozens or hundreds of simulations with the same parameter set to obtain statistical significant results. Finally, we can use the code in a way that we can introduce one or more parameter each time and separately. Therefore, our models are providing a completely complementary insight into the star cluster formation problem. Our initial conditions are a fractal distribution of stars (fractal dimension D=1.6; an example is shown in Fig.2) embedded in a smooth analytical background representing the residual gas of the embedded cluster. All of our model start with an overall SFE = 0.2. In terms of the old simulations this means that none of our initial distributions should form a bound cluster. To mimic very fast and maximal destructive gas expulsion we remove the background potential instantaneously.

The Erasure of Substructure

In a first experiment we asked how fast can we erase the substructure of our initial distribution of stars. We use small spherical clumps of stars instead of a fractal distribution and follow the dynamics counting how many clumps remain with time. The results are shown in Fig.3. We see that substructure is erased very fast, i.e. exponentially, but maybe not fast enough.

The Local Stellar Fraction (LSF)

Instead of the overall SFE, which does not vary with time even though we have a highly dynamical system, we find that the LSF (local stellar fraction = an instantaneous SFE, measured at the varying half-mass radius at exactly the time of gas-expulsion) is a better description to predict the final bound mass fraction f_bound, i.e. how many stars of the embedded cluster remain bound to each other after gas-expulsion. This is particularly true as soon as the initial fractal distribution had enough time to settle into a virial equilibrium, i.e. when we have a late gas expulsion. We are able to quantify this relation as:

where f_bound is the fraction of stars bound after gas expulsion and 0.14 is an off-set, determined through our simulations. Below this off-set we do not see any surviving cluster. The results are shown in Fig.4. We still see a lot of scatter in the results which may be due to a second parameter.

The chaotic regime

We dub the time of early gas-expulsion as the 'chaotic' regime, as our proposed parameter, the LSF no longer predicts the bound-mass fraction, Here we clearly see that the survival of our star cluster has to depend on a second parameter. We suspect that the virial ratio at time of gas expulsion plays an important role in our experiments.

A new time-scale

The young embedded star clusters are not in virial equilibrium, even though in most of our simulations the stars have virialised velocities inside the embedded cluster. But they are not forming a spherically symmetric distribution. What we see in our simulations is an initial collapse into the centre and an oscillation of the star distribution afterwards until we reach a final spherically symmetric and virialised distribution. Also the instantaneous virial state of our stellar distribution does oscillate.

Measuring time in initial crossing times of the system is not a correct measure of the dynamical state. We choose a new time-scale T_Q, which measures time in oscillations of the instantaneous virial parameter Q_f. With this new time-scale we can determine at which instant our distribution has exactly Q_f=0.5 and when do we see extremes in the virial parameters.

The importance of the virial ratio

We suspected that we have a second important parameter to take into account, namely the actual virial ratio of the stars at exactly the time of gas-expulsion, which we call Q_f, to distinguish with the virial ratio at which the embedded cluster is born (Q_i). We find that we can generalise Eq.1 in the following sense to again be able to predict the results of our simulations:

which reproduces Eq.1 within the error-bars. Still, for low LSF and Q_f values this relation raises very steeply so that with the remaining scatter, we still see in our results, the prediction of the final bound mass is impossible.

For the remaining scatter in our results we suspect that we have to take dQ/dt into account and the stored energy in small self-gravitating sub-clumps. We could exclude the clumpiness parameter C = MST / s. This parameter can predict the amount of scatter but not if a simulation falls above or below the trend (more below).

But this investigation is postponed to later as we already solved the problem to a better degree than any observational measurement, made in the near future, can distinguish (i.e. it may be possible to measure with some precision a half-light radius but the calculation of the actual virial ratio to a precision required to distinguish the outcome according to our relation is virtually impossible).

A Quick Theoretical Insight

Please find a short description of the theory written in this PDF document. The fitting curves are shown in figure 8.

Influence of the background

The remaining gas in our simplified simulations id modeled as a static background potential, which then is removed to mimic gas-expulsion. We investigated if the shape of this potential has any influence on the outcome of our results. We varied the potential from an uniform sphere to a highly concentrated Plummer distribution but our results showed that the answer of the star distribution (fractal distribution with fractal dimension 1.6) is almost completely independent to the shape of the gas potential. Simulations with the same initial conditions of the stars formed similar denser configurations before gas expulsion. This result was on one hand reassuring that our simplified models are quite insensitive to the choice of static potential and on the other hand allowed us to probe different regions of the LSF, as concentrated backgrounds have more remaining gas within the star distribution than e.g. an uniform sphere.

The Clumpiness parameter

The clumpiness parameter C describes how sub-structured a distribution of stars is. As long as C < 0.8 we have a clumpy, filamentary structure and C can be related with the fractal dimension of the distribution and if C > 0.8 we have a smooth distribution and C related directly to the shape (i.e. concentration) of the distribution. We find that in our simulations C changes from the fractal into the smooth regime at about the same time we change from the 'chaotic' into the well behaved regime. This happens at about 1.5 initial crossing times of the system or t_Q approx 2 (see below).

The Snowball method

To analyse the data and determine the bound mass fraction we developed a new tool, we call the 'Snowball method'. The method will be described in detail in a paper, which we are preparing, to appear in a refereed journal.

The influence of the IMF

Not truncated IMF

As stars do not have the same mass, we included now in our models a standard Kroupa IMF (initial mass function) for stars between 0.1 and 120 M_sun. Even though there is observational and theoretical evidence that the IMF in low-mass systems is truncated, not just stochastically (i.e. a 50 M_sun star cluster cannot harbour a 100 M_sun star and 30 M_sun star is highly unlikely to be found in a 500 M_sun star cluster), but that there is a certain maximum star mass for a given mass of the star cluster, we first included only a stochastic truncation into our models as we find it quite hard to determine which should be the mass of our embedded cluster, as long as it is still in a fractal state. We will postpone the investigation of truncated and varied IMFs to a later stage of the project.

We see in our results, so far, that the internal evolution of the star distribution is faster. Especially for dense systems, we find that the bound fractions measured late are well below the trend we determined for equal mass particles. A closer look to our results then showed that this is due to the fast post-gas-expulsion evolution of our low-N systems. Once we realised this fact, we started to measure the bound fraction directly after the gas-expulsion and we again recovered the trend shown in Eq.2.

Heavily mass-segregated models

As a first step into the investigation of mass-segregation we performed models which are heavily mass-segregated, i.e. we place the most massive star closest to the centre and the least massive one furthest out. If this is a valid assumption, is debatable and we will investigate also cases in which we have mass-segregated sub-clumps (see reasoning for truncated IMF). What we see is that, as predicted from stellar dynamics, the most massive stars form immediately a tight binary in the central area. The results then are a bit more counter-intuitive. Even though this hard binary ejects many small stars from the system, a system with that binary survives better than a system in which this binary gets ejected due to a three-body encounter. In this sense we get different results than our collaborators in the UK. We are still in the process to verify our results.

With many more simulations at hand, we see that the spread around our theoretical curve is much larger than in the equal mass particle simulations. If our theoretical model still holds will be subject to investigations, which are in process.

Advanced SPH models

The direct N-body-SPH code Hybrid-Seren is still not ready to produce results as the binary routines are not completely functional yet. But the description paper is published with R. Smith as co-author. In the meantime we make use of the AMUSE software package which combines multiple purpose codes with a Phyton envelope. With this we could combine a direct N-body code with an SPH code to perform the first few test simulations (without feedback).