The Long-term Evolution of Star Clusters Formed with a Centrally Peaked Star Formation Efficiency Profile

In Shukirgaliyev et al. (2017), we performed a set of direct N-body simulations to study the response of star clusters to instantaneous gas expulsion. We have considered clusters with different birth stellar masses (the stellar mass of a cluster at the time of instantaneous gas expulsion, M_⋆ = 3k, 6k, 10k, 15k, 30k ${M}_{\odot }$), different global SFEs (the star-forming gas mass fraction converted into stars before gas expulsion, SFE_gl = 0.10, 0.13, 0.15, 0.20, 0.25), and different mean volume densities (quantified by the ratio of half-mass radius to Jacobi radius, r_h/R_J = 0.025, 0.05, 0.07, 0.1, of the gas-free cluster). We have assumed that star clusters form according to the local-density-driven cluster formation model of Parmentier & Pfalzner (2013; i.e., stars have a density profile steeper than that of the residual gas at the time of gas expulsion). However, while Parmentier & Pfalzner (2013) infer the density profile of the residual gas and of the embedded cluster starting from that of the initial gas, we infer the density profile of the initial and residual gas from the assumed Plummer density profile of the embedded cluster. Therefore, in all model embedded clusters, stars follow a Plummer profile at the time of gas expulsion. Initial and residual gases do not, however, and are different for different SFEs (see Figure 2 of Shukirgaliyev et al. 2017). We recover the total and residual gas density profiles using the cluster formation model of Parmentier & Pfalzner (2013; see their Equations (19)–(20); for more details, see Section 2 of Shukirgaliyev et al. 2017). We have used the program mkhalo from the falcON package (McMillan & Dehnen 2007) to generate the initial conditions of our N-body models assuming that embedded clusters are in virial equilibrium with the total gravitational potential (i.e., stars + gas) at the time of gas expulsion. The considered model clusters are on a circular orbit in the Galactic disk plane at a Galactocentric distance of R_G = 8 kpc. The Galactic tidal field is given analytically by an axisymmetric three-component Plummer–Kuzmin model (Miyamoto & Nagai 1975) with the parameters as given in Just et al. (2009; their Equation (32) and their Table 1). We assume that all stars reach the zero-age main sequence at the time of instantaneous gas expulsion for simplicity, so the time span since gas expulsion equals the cluster age in our simulations. The initial mass function (IMF) of Kroupa (2001) with initial star mass limits of m_low = 0.08 M_⊙ and m_up = 100 M_⊙ has been applied. For more information about the initial conditions and the violent relaxation phase of cluster evolution, we refer our readers to Shukirgaliyev et al. (2017). By "violent relaxation" we mean the evolution of gas-free clusters from a supervirial state into a new state of quasi-equilibrium following instantaneous gas expulsion. We assume that the violent relaxation ends when model clusters stop (violently) losing their mass in response to gas expulsion as a result of an initial state of nonequilibrium. We note, however, that our model clusters may not have fully revirialized at that time (i.e., the end of violent relaxation, as we define it, does not coincide exactly with the time of cluster revirialization). In this paper, we refer to the end of violent relaxation as t = 30 Myr after instantaneous gas expulsion.

In this study, we continue our existing set of direct N-body simulations up to the full dissolution of the model star clusters in the tidal field of the Galaxy. We use only the Shukirgaliyev et al. (2017) models with r_h/R_J = 0.05 and that survive as bound clusters after violent relaxation (i.e., SFE_gl ≥ 0.15). We completed this initial model set with newly run SFE_gl = 0.17 models for some birth masses and with SFE_gl = 0.15 models for birth masses higher than what was considered in Shukirgaliyev et al. (2017; i.e., M_⋆ = 60k ${M}_{\odot }$ and 100k ${M}_{\odot }$, equivalent to N_⋆ = 105,554 and 174,257 stars, respectively, for a Kroupa 2001 IMF with m_low = 0.08 M_⊙ and m_up = 100 M_⊙). All model clusters considered in this study have evolved from the time of instantaneous gas expulsion until full dissolution. Therefore all of them bear the information about their formation conditions and violent relaxation. The full parameter space covered by our N-body simulations is provided in Table 1. All of our high-resolution, direct N-body simulations were performed using the ϕgrape-gpu code developed by Berczik et al. (2011) with single-stellar evolution recipes of sse by Hurley et al. (2000; for more details about ψ-grape-gpu see also Berczik et al. 2013). We have run our simulations on high-performance computing clusters: jureca (Jülich Supercomputing Centre 2016, Jülich, Germany), laohu (NAOC/CAS, Beijing, China), kepler (ARI/ZAH, Heidelberg, Germany), and the GPU part of bwForCluster: mls&wiso (Heidelberg Unversity, Heidelberg, Germany).

Table 1.  Model Cluster Parameters

N⋆	$\tfrac{{M}_{\star }}{{M}_{\odot }}$	SFEgl	nrnd	$\langle {F}_{\mathrm{bound}}\rangle$	$\left\langle \tfrac{{t}_{\mathrm{dis}}}{{10}^{9}\ \mathrm{yr}}\right\rangle$
(1)	(2)	(3)	(4)	(5)	(6)
5225	3000	0.15	21 (21)	0.07 ± 0.05	0.19 ± 0.11
5225	3000	0.17	1 (1)	0.21	0.52
5225	3000	0.20	1 (1)	0.29	0.71
5225	3000	0.25	1 (1)	0.53	1.07
10455	6000	0.15	26 (26)	0.06 ± 0.04	0.23 ± 0.16
10455	6000	0.17	1 (1)	0.21	0.66
10455	6000	0.20	3 (1)	0.31 ± 0.02	1.29
10455	6000	0.25	3 (1)	0.50 ± 0.03	1.56
17425	10000	0.15	26 (26)	0.06 ± 0.03	0.33 ± 0.15
17425	10000	0.20	1 (1)	0.38	2.11
17425	10000	0.25	1 (1)	0.51	2.55
26138	15000	0.15	22 (22)	0.08 ± 0.04	0.53 ± 0.24
26138	15000	0.17	1 (1)	0.18	1.02
26138	15000	0.20	3 (1)	0.33 ± 0.02	1.84
26138	15000	0.25	3 (1)	0.51 ± 0.01	2.86
52277	30000	0.15	16 (16)	0.08 ± 0.03	0.64 ± 0.32
52277	30000	0.17	1 (1)	0.16	1.37
52277	30000	0.20	1 (1)	0.34	3.58
52277	30000	0.25	1 (1)	0.56	5.49
104554	60000	0.15	15 (9)	0.08 ± 0.02	0.75 ± 0.35
174257	100,000	0.15	3 (3)	0.08 ± 0.01	0.96 ± 0.34

Note. Columns are as follows: (1) total number of stars, (2) birth mass, (3) global SFE, (4) number of random realizations per model (number of those calculated up to cluster dissolution), (5) mean bound mass fraction at the end of violent relaxation and its standard deviation when there is more than one random seed, (6) mean dissolution time.

Download table as:  ASCIITypeset image

In Shukirgaliyev et al. (2017), we mentioned that the bound mass fraction at the end of violent relaxation, F_bound (i.e., the stellar mass fraction remaining bound to the cluster at an age of 20–30 Myr), can vary by about 6%–10% of the birth mass for different random realizations of a given model. Such variations are almost negligible for model clusters with a high SFE_gl > 0.15, as they survive with more than 20% of the birth mass. However, the situation changes for our SFE_gl = 0.15 models, as they survive with the lowest bound mass fraction (about or even below 0.10), being the transition models from full destruction following gas expulsion to survival in the parameter space of SFE_gl. For such a low bound mass fraction, the variations mentioned above are significant. Therefore, we performed additional random realizations of the SFE_gl = 0.15 models to have better statistics for our study. In fourth column of Table 1 we give the numbers of random realizations per model and in parentheses the number of runs completed to cluster dissolution. Column 5 shows the mean bound mass fractions and their standard deviations measured at t = 30 Myr. The discussions and results on random realizations are presented further in Section 3.

Estimating cluster masses is a complex issue, from both an observational and theoretical point of view. Even the very definition of a star cluster varies through the literature (e.g., see review in Renaud 2018). Estimating the luminous mass of an observed star cluster is not straightforward, due to incompleteness issues and field star contamination. Estimating the dynamical mass of an observed cluster is hindered by the contribution of binaries to the cluster overall velocity dispersion. From a theoretical point of view, the mass determination is not straightforward, due to the unknown second integral of motion.

In Shukirgaliyev et al. (2017), we refer to the cluster mass as the stellar mass within one Jacobi radius, R_J, which is also supposed to be the bound mass. To calculate the Jacobi radius of a cluster at a given age, we start with its value at the time of gas expulsion and with the stellar mass it contains using Equation (13) from Just et al. (2009). Then we recalculate the Jacobi radius using the mass within the previously defined Jacobi radius. We iterate until the Jacobi radius converges. Our Jacobi radius calculation method works only if we define the cluster density center correctly, which can prove an issue for the following reasons. In our N-body simulations, we keep track of all stars, even those that have definitively escaped the cluster in which they initially formed. Therefore, escaped stars live in the tidal tails of our model clusters, which can extend so far as to wrap around the Galaxy, making it difficult to define the cluster center as a center of mass of all stars. That the tidal tails contain epicyclic overdensities also creates difficulties in finding the exact cluster center. That is, we have struggled to find the correct cluster center in the late stages of cluster evolution, as clusters are becoming low-mass and diffuse objects. Their low densities can then be comparable to that of the surrounding field, or to the epicyclic overdensities of the tidal tails. Our algorithm can thus incorrectly identify the center of a tidal tail overdensity as the cluster center. In order to prevent this, to calculate the cluster center, we use only those stars that were within 2R_J of the center in the previous N-body simulation snapshot. This allowed us to steer clear of the epicyclic overdensities of the tidal tails. Yet, it remains difficult to correctly define the cluster center when it consists of few stars, has an extended core, a low volume density, or presents substructures. Usually, this happens in the last stages of evolution of star clusters, before they are fully destroyed by the Galactic tidal field. In such clusters, relatively small shifts of the assumed cluster center can lead to significantly wrong mass estimates because of a wrong Jacobi radius, R_J.

In order to avoid this, we assume that clusters are dissolved if they have less than 100 ${M}_{\odot }$ left within two Jacobi radii:

$$\begin{eqnarray}&&{t}_{\mathrm{dis}}=t({M}_{2J}\lt 100{M}_{\odot }),\end{eqnarray} \tag{ 1 }$$

where M_2J is the stellar mass enclosed within two Jacobi radii.

We have checked how long some of our model clusters can live beyond the dissolution time, t_dis, defined by Equation (1). Typically the difference is about a few tens of megayears, which is not significant, especially for clusters whose lifetime scales up to a gigayear.

In this contribution, we consider two types of cluster mass estimates. One is the Jacobi mass (or "bound mass"), M_J, which is the stellar mass within one Jacobi radius, R_J. The other one, which we refer to as the "extended mass" M_2J, is the stellar mass within 2R_J. The second mass estimate is important for young clusters as they are surrounded by an envelope of unbound stars (Elson et al. 1987), most of them located beyond one Jacobi radius but still within 2R_J. Such envelopes persist for many megayears, as Figure 1 will show.

Zoom In Zoom Out Reset image size

These envelope stars can count toward the mass of clusters in extragalactic studies, where a membership analysis is impossible, and the cluster mass is estimated by fitting the cluster surface brightness profile.

Figure 1 visualizes the evolution of star clusters with a birth mass of M_⋆ = 15k ${M}_{\odot }$ in the form of volume density maps. Each point in these plots corresponds to the position of one star projected onto the Galactic disk plane. The colors correspond to the local volume density obtained by a nearest-neighbor scheme with N_nb = 50 neighbor stars. The coordinate system is centered on the density center of the clusters, and the two circles correspond to 1R_J and 2R_J. The left and middle panels depict two random realizations of the model with SFE_gl = 0.15, and the right column corresponds to SFE_gl = 0.25. In each column, different snapshots are presented (t = 0, 5, 10, 30, 70 Myr).

In Figure 2 we show the radial volume density profiles of the same three model clusters as in Figure 1 calculated at the end of violent relaxation (t = 30 Myr, top panels) and at a later time when clusters are almost cleared of their envelope stars (t = 70 Myr, bottom panels). Each point represents the density at the location of one star, calculated using a 50-nearest-neighbor scheme. In each panel, the vertical dashed lines correspond to one and two Jacobi radii, and the horizontal dotted line corresponds to the mean stellar density within one Jacobi radius, that is, $\langle {\rho }_{J}\rangle ={M}_{J}/(4/3\pi {R}_{J}^{3})$. According to Equation (13) of Just et al. (2009), for a given environment, that is, for a fixed circular orbit in the Galactic disk plane, the mean density within one Jacobi radius, $\langle {\rho }_{J}\rangle$, is constant and independent of cluster parameters, since ${R}_{J}\propto {M}_{J}^{1/3}$. In our case, for a circular orbit with R_G = 8 kpc, the mean density within one Jacobi radius is $\langle {\rho }_{J}\rangle$ ≈ 0.1 ${M}_{\odot }$ pc⁻³.

Zoom In Zoom Out Reset image size

In the left panels of Figure 2, the region between 1R_J and 2R_J is well populated by stars such that the distant observer, when calculating the cluster mass by fitting its projected density profile, will use the envelope stars too. For the middle and right columns, however, the envelope stars do not contribute much to the cluster mass, as evidenced by the density contrast between the central part and the outskirts of the clusters (see Figures 1 and 2). These density contrasts are about one to two orders of magnitude in the left panels, about three orders of magnitude in the middle panels, and about four to five orders of magnitude in the right panels.

In this section, we study the relation between cluster dissolution time and mass. From the random realizations of model clusters with SFE_gl = 0.15, we find that the differences in bound mass fraction at the end of violent relaxation also result in different star cluster lifetimes. Figure 4 presents the dissolution time t_dis of star clusters as a function of the bound mass fraction F_bound at the end of violent relaxation. The cluster dissolution time is defined according to Equation (1). The color coding and symbol coding correspond to cluster birth mass and global SFE, respectively (see the key). The general trend is that the higher the bound mass fraction, the longer a cluster lives.

Zoom In Zoom Out Reset image size

However, the large scatter and the fact that Figure 4 considers various birth masses do not give us much more information about the relation between cluster dissolution time and its mass.

The cluster dissolution time, t_dis, as a function of cluster "initial" mass is presented in Figure 5.

Zoom In Zoom Out Reset image size

By "initial" mass, M_init, we mean here the cluster mass once violent relaxation is over, that is, when the long-term secular evolution starts. Our cluster "initial" masses are therefore measured at t = 30 Myr (top panels) and t = 70 Myr (bottom panels) and are lower than the birth masses given in Table 1 (the ratio between Jacobi and birth masses is the bound fraction F_bound). We have done so to be consistent with studies that ignore the violent relaxation phase of cluster evolution when inferring the cluster dissolution time as a function of cluster mass (for instance, by considering only clusters older than 10–30 Myr).

When plotting the cluster dissolution time as a function of cluster initial mass, not only do we measure the mass at two different cluster ages, but we also consider two definitions of cluster masses in terms of cluster spatial coverage. Specifically, the cluster initial mass is defined either as the Jacobi mass M_J (left panels) or as the "extended mass" M_2J (right panels). We therefore investigate four different cases in total. The color coding shows the birth mass of the model clusters, while the different symbols correspond to different SFE_gl (see the key). The scatter arising from the random realizations of a given model (given birth mass and given SFE_gl) is therefore illustrated by symbols of a given color and of a given type.

The solid line, with the shaded area accounting for the error bars, corresponds to the MDD relation of Boutloukos & Lamers (2003):

$$\begin{eqnarray}&&{t}_{\mathrm{dis}}={t}_{4}^{\mathrm{dis}}{\left(\displaystyle \frac{{M}_{\mathrm{init}}}{{10}^{4}{M}_{\odot }}\right)}^{\gamma },\end{eqnarray} \tag{ 2 }$$

where M_init is the cluster "initial" mass, and t₄^dis is the dissolution time for a cluster with "initial" mass of M_init = 10⁴${M}_{\odot }$. Here we recall that M_init is the equivalent of the Jacobi mass, M_J, or the extended mass, M_2J, at t = 30 or 70 Myr, but not of the birth mass, M_⋆. The values of ${t}_{4}^{\mathrm{dis}}=(1.3\pm 0.5)\cdot {10}^{9}\ \mathrm{yr}$ and γ = 0.6 are taken from Lamers et al. (2005a) for the solar neighborhood. The dashed and dash-dotted lines show the best fits to our high-SFE (SFE_gl ≥ 0.20) and low-SFE (SFE_gl = 0.15) model clusters, respectively. The bold red curve in the lower right panel connects, for each cluster birth mass (i.e., for each symbol color), the medians of SFE_gl = 0.15 model random realizations.

The overall impression from Figure 5 is that star clusters dissolve in agreement with MDD, although with some significant scatter. Especially, model clusters formed with a relatively high global SFE (SFE_gl ≥ 0.20, open squares and triangles in Figure 5) nicely follow an MDD relation (dashed line), with a slope a bit steeper (γ ∼ 0.7) than that given by Boutloukos & Lamers (2003) and a dissolution time longer (${t}_{4}^{\mathrm{dis}}\sim 3.8\,\mathrm{Gyr}$, lower left panel of Figure 5) than the observational estimate of Lamers et al. (2005a) for the solar neighborhood (solid line). Our result is consistent with other theoretical works, that Roche-volume-filling or underfilling clusters in virial equilibrium dissolve in a mass-dependent way. However, the cluster dissolution time obtained from our high-SFE models (${t}_{4}^{\mathrm{dis}}\sim 3.8\ \mathrm{Gyr}$) is almost a factor of two shorter than that predicted by Baumgardt & Makino (2003), who found ${t}_{4}^{\mathrm{dis}}=6.9\ \mathrm{Gyr}$. This can result from different models for the Galactic gravitational potential: while we consider an axisymmetric three-component model (Just et al. 2009), Baumgardt & Makino (2003) consider a spherically symmetric logarithmic potential. Another and probably more crucial reason is the "initial" density profile of the star clusters: our density profiles at t = 30 Myr differ from the virialized King models (with W₀ = 5 and W₀ = 7) used by Baumgardt & Makino (2003).

We do not notice any significant difference for the high-SFE models between the four panels. Therefore, for the high-SFE models, how the cluster initial mass is defined (at $t=30\ \mathrm{or}t=70\ \mathrm{Myr};$ bound mass, M_J, or extended mass, M_2J) hardly influences the corresponding predicted MDD relation.

In contrast to high-SFE ones, low-SFE models show broad scatter and differ from panel to panel. The scatter is maximum for extended mass at t = 30 Myr (upper right panel) and minimum for Jacobi mass at t = 70 Myr (lower left panel). Since the difference between the Jacobi mass and the extended mass is larger for low-SFE clusters and negligible for high-SFE clusters, the initial masses are characterized by a broader scatter in the right panels than in the left panels of Figure 5. When the clusters have evolved a bit more in time, the envelope gets cleaned up by an age of t = 70 Myr, and the extended mass of our model clusters then becomes comparable to their Jacobi mass (bottom panels of Figure 5), although the scatter persists. Our best fits to low-SFE models provide a dissolution time in agreement with Lamers et al. (2005a), ${t}_{4}^{\mathrm{dis}}\sim 1.3\ \mathrm{Gyr}$, but a slightly shallower slope, γ ∼ 0.5 (lower left panel of Figure 5). For the SFE_gl = 0.15 models, we can see that the scatter in bound mass fractions can yield significant differences in the cluster lifetime (open circles of a given color). As shown in Figures 1 and 2 (left and middle panels), this is the result of the development with time of markedly different density profiles. Additionally, if we consider a vertical bin embracing cluster masses from 2 × 10³ to 4 × 10³${M}_{\odot }$ in the right panels of Figure 5, we can see that such clusters can dissolve as fast as within 100–200 Myr or can live as long as about a gigayear.

For the sake of clarity, Figure 6 shows each low-SFE model (i.e., given birth mass and global SFE) represented by its mean and standard deviation.

Zoom In Zoom Out Reset image size

The initial mass is defined as the extended mass, M_2J, and the age is t = 70 Myr. Figure 6 is thus equivalent to the lower right panel of Figure 5, apart from the low-SFE models being represented by mean values. Now the results for the SFE_gl = 0.17 models can be seen clearly as open diamonds with the corresponding fit having ${t}_{4}^{\mathrm{dis}}=1.8\ \mathrm{Gyr}$ and γ ∼ 0.5 (dotted line). The other lines are as in the lower right panel of Figure 5, namely, the MDD relation of Boutloukos & Lamers (2003) and the best fits to the SFE_gl ≥ 0.20 and SFE_gl = 0.15 models. The combination of these three groups of star clusters can yield a relation close to that of Boutloukos & Lamers (2003), especially if low-SFE clusters dominate the cluster population. This is a real possibility for the solar neighborhood, since observations of nearby gas-embedded clusters tell us that low-SFE embedded clusters are more common than high-SFE ones (Lada & Lada 2003; Evans et al. 2009; Peterson et al. 2011; Kennicutt & Evans 2012; Kainulainen et al. 2014).

Gieles et al. (2006) noted that the dissolution time of a 10⁴${M}_{\odot }$ cluster in the solar neighborhood as inferred from observations by Lamers et al. (2005a; ${t}_{4}^{\mathrm{dis}}=1.3\pm 0.5\ \mathrm{Gyr}$) differs by about a factor of five from what was predicted by Baumgardt & Makino (2003; ${t}_{4}^{\mathrm{dis}}=6.9\ \mathrm{Gyr}$). They suggested that the discrepancy can be eliminated by accounting for the influence on star clusters of encounters with giant molecular clouds. Here, we propose yet another channel to explain the shorter dissolution time of Lamers et al. (2005a), namely, that the star clusters of the solar neighborhood are predominantly the survivors of embedded clusters formed with a global SFE of SFE_gl ≈ 0.15. We expect that most of these clusters will dissolve by the time they reach an age of 1 Gyr. In contrast to other theoretical studies, in which compact clusters in virial equilibrium are considered as initial conditions, our model clusters have experienced violent relaxation, which is a natural process affecting the evolution of young clusters. With our approach, we are thus able to simulate the evolution of clusters that have survived violent relaxation as bound, but diffuse, objects. They dissolve faster than their compact counterparts, even for otherwise equal "initial" masses, and we probably observe them as open clusters in the solar neighborhood.

If the low-SFE clusters dominate the cluster census of the Galactic disk, then not many clusters are able to live beyond 1 Gyr, as we can see from Figure 7.

Zoom In Zoom Out Reset image size

Figure 7 shows the histogram of cluster dissolution times for low-SFE (SFE_gl = 0.15, blue) and higher-SFE (SFE_gl ≥ 0.17, orange) clusters. Each histogram is normalized such that the sum of all bins equals unity. The combination of both distributions, each with its own contribution, will provide the distribution of dissolution times of a cluster population. The latter can give us a hint about the shape of the corresponding cluster age distribution. A cluster population dominated by low-SFE clusters should feature a peak in cluster logarithmic age distribution (dN/d log t) earlier than 1 Gyr (if the cluster formation rate is constant), since most low-SFE clusters die before 1 Gyr. We coin this process the teenage mortality of star clusters. The peak at about a few hundred megayears in the age distribution of solar-neighborhood star clusters has been discussed in Lamers et al. (2005a) and Piskunov et al. (2018).

An interesting question is whether the dissolution time distribution of low-SFE clusters extends significantly beyond 1 Gyr, if our simulations were to include clusters more massive than those modeled so far. Here we remind the reader that the initial cluster mass function is a power law with a slope of −2, that is, $\tfrac{{dN}}{{dm}}\propto {m}^{-2}$, equivalent to $\tfrac{{dN}}{d\mathrm{log}m}\propto {m}^{-1}$. This implies that for a mass spectrum with a lower mass limit of 100 ${M}_{\odot }$, star clusters more massive than 10⁴${M}_{\odot }$ (10⁵${M}_{\odot }$) represent 1% (0.1%) of the cluster census.

In addition, to produce a low-SFE cluster with M_init = 10⁵${M}_{\odot }$ (10⁶${M}_{\odot }$), we need a molecular clump with a gas mass of 10⁷${M}_{\odot }$ (10⁸${M}_{\odot }$). This is because the "initial" mass of such clusters is only about 1% (F_bound × SFE_gl = 0.07 × 0.15) of the total mass of the star-forming clump. In our Galaxy, very massive low-SFE clusters should thus be very rare, if they exist at all, because the mass function of molecular clumps is also a decreasing power law with an index of −1.7 to −2.

Finally, according to Rahner et al. (2017), the higher the mass of a star-forming clump, the higher the SFE it requires to be destroyed by its newly formed star cluster. Should this not be the case, the recollapse of the stellar-feedback-driven gaseous shell may happen and produce more stars, enhancing the SFE (Rahner et al. 2018). Therefore, the distribution presented in Figure 7 is robust and will not change significantly if more massive low-SFE clusters were included.

An interesting trend emerges when we consider SFE_gl = 0.15 models in Figure 6. If we ignore the low-mass end (≤10³${M}_{\odot }$, i.e., the mass range for which the cluster sample is incomplete in most extragalactic studies) and if we take into account that, in extragalactic observations, the masses of young clusters can be overestimated because of the contribution of an envelope of unbound stars (i.e., the cluster "initial" mass has to be defined as our extended mass, M_2J), an apparent MID mode may be emerging. For initial mass M_2J between 10³ and 10⁴${M}_{\odot }$, low-SFE clusters actually show similar mean (and median) dissolution times, with even increasing standard deviations. On top of that, if we consider model clusters formed with very low SFE_gl < 0.15, which do not survive instantaneous gas expulsion and the resulting violent relaxation (t_dis < 30 Myr), but which are still observable as young clusters, then the effect of an apparent MID can be further strengthened.

However, since we do not cover that large a range of cluster masses and do not have that high a number of random realizations for good statistics (especially at high mass), we cannot argue firmly in favor of a mass-independent dissolution relation for low-SFE star clusters. Here, we stress that with a mean bound fraction at the end of violent relaxation of F_bound = 0.07, modeling a star cluster with an "initial" mass >10⁴${M}_{\odot }$ and that formed with SFE_gl = 0.15 requires a birth mass of at least 10⁵${M}_{\odot }$. We are currently expanding our sample for M_⋆ = 10⁵${M}_{\odot }$. The modeling of even more massive clusters is highly desirable for comparison with the results of observational extragalactic studies. Nevertheless, we can firmly say that we have found a large scatter in the relation between the cluster "initial" masses and the cluster dissolution times for low-SFE clusters and for a given birth mass. This scatter results from the massive-star-driven stochastic effects taking place during violent relaxation. Such effects yield, for a given birth mass and given SFE_gl = 0.15, different density profiles (and hence different degrees of cluster compactness) and different bound fractions at the end of violent relaxation (and hence different cluster "initial" masses).

In summary, so far we have found that clusters formed with a high SFE_gl (>0.15) dissolve in an MDD regime. Clusters formed with a low SFE_gl (=0.15) also dissolve in an MDD regime, albeit with a significant scatter. Their dissolution time is comparable to that observationally inferred by Lamers et al. (2005a). There is a strong mass-dependent upper limit to the cluster dissolution time, which means that, for a given environment, low-mass clusters cannot live as long as their high-mass compact counterparts. However, some high-mass clusters can dissolve as quickly as low-mass ones in the same environment.

Investigating further the parameters of our model clusters, we have found a correlation between the cluster dissolution time, t_dis, and the Roche volume filling factor measured at the end of violent relaxation, although with significant scatter (Figure 8).

Zoom In Zoom Out Reset image size

The Roche volume filling factor is defined as the ratio of the half-mass radius to the Jacobi radius, and the half-mass radius refers to the half-Jacobi mass radius or the half-extended mass radius (left and right panels of Figure 8, respectively). The filling factors are calculated at t = 30 and 70 Myr in the top and bottom panels, as in Figure 5. The color coding depicts the Jacobi mass, M_J (left panels), and the extended mass, M_2J (right panels), indicated by the common color bar at the right. Note that the color coding is different from Figure 5, where it refers to the birth mass of clusters. The shapes of markers still show the global SFEs. With gray arrows, we indicate those M_⋆ = 15,000 ${M}_{\odot }$ clusters presented in Figures 1 and 2, namely, from right to left, two random realizations of the SFE_gl = 0.15 model (open circles) and one random realization of the SFE_gl = 0.25 model (open triangle). The indicated model clusters here are in reverse order to that of the panels in Figure 1.

The correlation is such that for a given "initial" mass, the higher the filling factor, the shorter the dissolution time, as found in Ernst et al. (2015). However, a comparison between our results and those of Ernst et al. (2015) is not fully self-consistent. First, our model clusters are still expanding and have not yet returned to virial equilibrium when they overfill their Roche volume. In contrast, model clusters of Ernst et al. (2015) are initially in virial equilibrium while overfilling their Roche volume. Second, our star clusters can present different density profiles (shallow or steep, with an extended or compact core) after violent relaxation, so the Roche volume filling factor, as defined here, cannot characterize them universally.

We also find a correlation between the cluster life expectancy, t_dis, and the volume density contrast between cluster center and outskirts, defined as the ratio of the central volume density, ρ_c, to the volume density at the Jacobi radius, ρ_J:

$$\begin{eqnarray}&&{\rho }_{c}/{\rho }_{J}=\displaystyle \frac{\rho (r=0)}{\rho (r={R}_{J})}.\end{eqnarray} \tag{ 3 }$$

This is shown in the upper left panel of Figure 9 as measured at the end of violent relaxation (t = 30 Myr). The smaller the density contrast, the faster the cluster dissolves and vice versa. The correlation is the tightest if we consider the density contrast of slightly more evolved clusters (t = 70 Myr, top right panel of Figure 9), when the envelope at R_J < r < 2R_J has been cleaned up of most of its stars.

Zoom In Zoom Out Reset image size

The color coding corresponds to the Jacobi mass at the corresponding ages. The best fits, in the form of

$$\begin{eqnarray}&&{t}_{\mathrm{dis}}=3.4\left[{\left(\displaystyle \frac{{\rho }_{c}}{{\rho }_{J}}\right)}^{0.25}-1\right]\times {10}^{8}\ \mathrm{yr}\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&{t}_{\mathrm{dis}}=1.1\left[{\left(\displaystyle \frac{{\rho }_{c}}{{\rho }_{J}}\right)}^{0.37}-1\right]\times {10}^{8}\ \mathrm{yr},\end{eqnarray} \tag{ 5 }$$

are shown with black curves in the top left and top right panels, respectively. For high density contrast (ρ_c/ρ_J > 10³) clusters, the top panels of Figure 9 show that their dissolution times are mass dependent, which is again consistent with earlier studies of dissolution of Roche volume filling or underfilling clusters.

When a cluster has a shallow density profile, the shrinking of its Jacobi radius due to cluster mass loss (e.g., stellar evolution, tidal stripping) will leave outside the new Jacobi radius a number of stars higher than in the case of a steeper density profile. Such clusters therefore dissolve more quickly than those with a steep density profile. Alternatively, if a cluster has developed a large core, comparable to the Roche volume, this also leads to the total destruction of the cluster. Such clusters are located at the very left of the top panels in Figure 9, with low density contrasts and therefore short lifetimes.

In order to see better the correlation between these three parameters, ρ_c/ρ_J, t_dis, and M_J, we swapped the dissolution time for the Jacobi mass in the bottom panels. The bottom panels of Figure 9 show the Jacobi mass, M_J, as a function of the density contrast, ρ_c/ρ_J, color coded by the dissolution time, t_dis. As now seen clearly here with a more or less nice pattern, the dissolution time depends on both "initial" mass and density contrast, wherein short dissolution times of massive clusters are explained by their low density contrast (low concentration) and vice versa. This is also consistent with other theoretical works, where the evolution of globular clusters with initially low concentration has been discussed (see, e.g., Fukushige & Heggie 1995; Takahashi & Portegies Zwart 2000; and Vesperini & Zepf 2003)

As we see from the foregoing results, when studying the long-term secular evolution of star clusters, it is worth estimating their "initial" mass at t = 70 Myr after gas expulsion in order not to be biased by the processes taking place during violent relaxation. The latter is very different from the subsequent long-term evolution of star clusters. As we showed in Shukirgaliyev et al. (2017), clusters formed with very low global SFE (<0.15) dissolve during violent relaxation, and their dissolution time is independent of their birth mass.

We remind the reader that we have modeled the clusters in a given tidal-field environment (i.e., star clusters have the same mean density at the time of instantaneous gas expulsion), and that they have the same stellar density profile at the time of instantaneous gas expulsion. All relations we have found could thus be affected if birth conditions are different. Star clusters that at the time of gas expulsion are more compact or more diffuse than our models may have evolutionary tracks different from those analyzed here, forming denser or more diffuse bound clusters after violent relaxation.