© The Authors 2025

1 Introduction

Open clusters (OCs) are fundamental building blocks of galaxies and serve as important tracers of galactic structure and evolution (Lada & Lada 2003; Portegies Zwart et al. 2010; Feigelson et al. 2013). Understanding their formation, dynamical evolution, and survival rates is crucial for interpreting the star formation history and chemical enrichment of the Milky Way (Krumholz et al. 2019). However, the processes governing the evolution of OCs, particularly during their early stages, are still not understood.

The early evolution of OCs is complex due to the interplay of various physical processes. OCs form within giant molecular clouds through the gravitational collapse and fragmentation of gas (McKee & Ostriker 2007). During this process, only a fraction of the gas is converted into stars, while the remainder persists as residual gas in the cluster environment. The star formation efficiency (SFE) – the fraction of gas converted into stars – is a critical parameter that significantly influences the ability of clusters to remain bound after the expulsion of the residual gas.

The gas expulsion process, driven by feedback from massive stars through stellar winds, ionizing radiation, and radiation pressure, occurs rapidly compared to the dynamical timescale of the cluster, well before supernova events take place (Krause et al. 2020; Lewis et al. 2023). This rapid loss of binding mass can lead to cluster expansion and potential dissolution via violent relaxation. Accurately modelling this process is challenging due to the complex interplay between stellar dynamics, gas hydrodynamics, and stellar feedback (Krumholz et al. 2019; Rahner et al. 2017, 2019; Pellegrini et al. 2020a,b).

Several studies have investigated this phenomenon using combined hydrodynamics and N-body simulations (Wall et al. 2019; Li et al. 2019; Fukushima & Yajima 2021; Fujii et al. 2021a,b; Guszejnov et al. 2022; Rieder et al. 2022; Rantala et al. 2023, and many others). However, these methods are computationally expensive, limiting their applicability to modelling a large number of clusters with diverse initial conditions. Interestingly, Fujii et al. (2021a) noted that simplified initial conditions with instantaneous gas expulsion can be a reasonable approximation when focusing on the long-term evolution of clusters.

Various approaches employing N-body simulations with the assumption of instantaneous gas expulsion at the outset have been developed (e.g. Baumgardt & Kroupa 2007; Goodwin & Bastian 2006; Smith et al. 2011; Shukirgaliyev et al. 2017, 2021). These studies have demonstrated that clusters with low SFEs are more likely to dissolve rapidly, while those with higher SFEs have a greater chance of survival. Nevertheless, the precise relationship between SFE, gas expulsion timescales, and cluster survival remains an active area of research (Lee & Goodwin 2016; Farias et al. 2018; Li et al. 2019; Shukirgaliyev et al. 2021).

Of particular interest is the method introduced by Shukirgaliyev et al. (2017), which employs more physically motivated initial conditions for the density distribution of stars and gas prior to instantaneous gas expulsion. This approach utilises the local-density-driven clustered star formation model of Parmentier & Pfalzner (2013) to reconstruct the density profile of the residual gas while the cluster is still embedded within it. The model assumes that star formation occurs within a centrally concentrated gas clump, with a constant efficiency per free-fall time, resulting in a star cluster (SC) with a steeper, more centrally concentrated density profile than the gas. This configuration allows even low-SFE clusters to survive the aftermath of gas expulsion by reducing their susceptibility to the gas potential. In particular, Shukirgaliyev et al. (2017) demonstrated that model star clusters with a Plummer profile could survive gas expulsion with a global SFE as low as 15.

In this work, we investigated the long-term evolution of SC with peaked star formation in the centre, following instantaneous gas expulsion under various Galactic orbital conditions. We adopted the initial conditions for cluster formation based on Shukirgaliyev et al. (2017), extending our study to elliptical and non-planar orbits. This extension enables us to explore a broader range of realistic scenarios, providing a more comprehensive understanding of cluster evolution within the Milky Way potential.

Section 2 details our methods and simulations, including the Milky Way potential employed and the orbital configurations of the model clusters. We also describe the initial conditions and the simulations conducted. In Sect. 3, we present our findings on mass loss in the model clusters. Additionally, we develop an analytical framework to characterise the rate and timescale of gas expulsion in Sect. 4. Section 5 compares our results with observations, particularly in terms of the Roche-lobe filling factors of Galactic OCs. Finally, we summarise our results and discuss their broader implications in Sect. 6.

2 Methods and models

For the generation of our SC initial conditions, we employed a multi-step approach comprising several key elements: the Milky Way (MW) rotation curve and tidal field; the calculation of the orbital motion of OCs within the Galaxy’s external tidal field; the generation of non-virial equilibrium initial OC models with varying global star formation efficiencies SFEs; and the combination of OC initial data with a chosen initial mass function (IMF) and physical normalisation of the initial conditions. These steps produce a set of initial OC conditions suitable for direct N -body modelling of our stellar systems.

2.1 Galactic tidal field

For this study, we adopted a static composite Galactic tidal field to model the external influences on the OCs. Specifically, we used the well-known Miyamoto-Nagai potential (Miyamoto & Nagai 1975), which combines a spherically symmetric Plummer model (Plummer 1911) with an axisymmetric Kuzmin disc (Kuzmin 1955). The complex MW potential is represented as a combination of three independent components: the bulge, disc, and halo. This is written as follows: $${\text{Φ}}_{tot}\equiv {\text{Φ}}_B+{\text{Φ}}_D+{\text{Φ}}_H.$$(1)

Each of these components can be described by the equation $${\text{Φ}}_{B,D,H}(R,Z)=-\frac{G\cdot M_{B,D,H}}{\sqrt{R^2+{\left(a_{B,D,H}+\sqrt{b_{B,D,H}^2+Z^2}\right)}^2}},$$(2)

where a and b are the scale radii of the flattened models for each component. Here, R represents the cylindrical radial distance and Z denotes the axial coordinate.

A similar axisymmetric Galactic potential was used in our previous simulations of the dynamical evolution of MW OCs (Just et al. 2009; Kharchenko et al. 2009; Ernst et al. 2010). Compared to these earlier works, the mass of the halo has been slightly adjusted to achieve a better fit to the circular velocity at the updated location of the Sun in the Galaxy’s rotation curve.

We present our parameters, including the bulge, disc, and halo masses, along with their characteristic scales for the static potential marked as FIX. We compare them with other proposed Milky Way potentials from Bennett et al. (2022) and Bland-Hawthorn & Gerhard (2016) in Table 1.

2.2 Orbital motion of OCs

For our set of dynamical simulations, we selected two families of galactocentric orbits with apocentres (r_apo) of 10 and 3 kpc. For these families, we used orbital eccentricities of e = 0 (circular orbit), 1/3, and 1/2. We also created models with the same orbital parameters but with an off-plane shift in the initial position, above the Galactic disc plane at Z = 200 pc and Z = 60 pc (corresponding to their respective r_apo values). As an illustration, in Fig. 1, we show two elliptical orbits with r_apo = 10 kpc and 3 kpc for e = 1/3. By using this wide range of orbital parameters, we aim to reproduce a mixture of possible OC orbital configurations in the Galactic disc near the Sun. The set of models with r_apo = 3 kpc is added specifically to investigate the potential differences between OCs near the Galactic centre and those in the solar vicinity.

Fig. 1 Orbital evolution of OC with e = 1/3 is shown in the X-Y and R- Z planes (where R is the planar galactocentric radius). The integration time is 1 Gyr, as indicated by the colour bar. Upper panel: orbit with rapo = 10 kpc (Z = 200 pc). Lower panel: orbit with rapo = 3 kpc (Z = 60 pc).

2.3 Generation of OCs’ non-virial initial conditions

The non-virial initial conditions for open clusters (OCs) prior to gas expulsion, as proposed by Shukirgaliyev et al. (2017), have been successfully implemented in a series of studies (Shukirgaliyev et al. 2018, 2019, 2021; Bissekenov et al. 2024). The basic idea of the proposed approach is straightforward: we assume dynamical equilibrium for the self-gravitating two-component system, consisting of a mixture of gas and stars. By assuming a theoretical stellar density distribution and using the approximation of local dynamical equilibrium between self-gravity and the velocity dispersion of the gas-star mixture in each shell, we derive the corresponding equilibrium gas distribution. A key physical assumption in this framework is that star formation occurs locally within each shell on the local free-fall timescale. Using this physical premise, Parmentier & Pfalzner (2013) developed a semi-analytical model for open cluster formation from centrally concentrated, spherically symmetric gas clumps with a constant SFE ϵ_ff per free-fall time t_SF.

For the stellar-mass density distribution, we used the well- known Plummer model (Plummer 1911): $${\rho }_{\star }(r)=\frac{3\cdot M_{\star }}{4\pi a_{\text{p}}^3}{\left(1+\frac{r^2}{a_{\text{p}}^2}\right)}^{-5/2},$$(3)

where M_⋆ is the cluster’s total stellar mass, and a_p is the so- called Plummer radius, which corresponds to the projected halfmass radius of the star cluster.

By varying the SFE per free-fall time and the star formation duration, we created a set of simulations with different global SFEs, specifically 0.15, 0.17, 0.20, 0.25, 0.30, and 0.35. In Shukirgaliyev et al. (2017), it was shown that a minimal global SFE of 0.15 is necessary to form a bound stellar system after instantaneous gas expulsion from the cluster. Our range of modeled global SFEs is in good agreement with observations of open clusters in the Solar vicinity (Higuchi et al. 2009; Murray 2011).

As an additional set of runs for comparison with previous classical works using initial conditions in virial equilibrium for stellar systems, we also generated a special configuration with a maximum SFE of 0.99. These initial conditions mimic the formation of purely stellar, gas-free OCs (see Fig. 2).

Fig. 2 Concentration profiles of cluster’s stellar (n⋆) and gaseous (n𝑔as) material are shown for different SFEs, ranging from 0.15 to 0.35 and 0.99. The dotted and dash-dotted black vertical lines represent the radii enclosing 50% and 90% of the cluster’s stellar mass, respectively. The dashed black vertical line marks the simulated cluster boundary, rcluster, which is defined as 20 ⋅ ap. For the rapo = 10 kpc and λ0 = 0.09 models (see Sect. 2.4), this corresponds to ~19.5 pc.

2.4 Cluster IMF and physical normalisation

The stellar masses are randomly drawn from an adapted Kroupa IMF (Kroupa 2001), with lower and upper mass limits of 0.08 – 100 M_⊙. For the stellar metallicity, we adopted the classical solar metallicity value of Z_⊙ = 0.02 (Grevesse & Sauval 1998). The total stellar mass of our theoretical OC is set to 6000 M_⊙, representing a typical mass for observed young clusters (Liu & Pang 2019; Hunt & Reffert 2023).

The initial Roche volume filling factor, λ₀, is defined as the ratio of the half-mass radius to the Jacobi radius (see definitions of Jacobi mass and Jacobi radius in Just et al. 2009), where the half-mass radius corresponds to the radius enclosing half the Jacobi mass. Based on earlier studies of OCs (Ernst et al. 2010, 2015), we adopt a standard value of λ₀ = 0.09 and, in some cases, a value two times smaller, λ₀ = 0.045. Both values correspond to tidally under-filled clusters.

Using the theoretical OC initial mass and distance from the Galactic centre, we derived the half-mass radius of clusters for the selected λ₀. For instance, for r_apo = 10 kpc and a cluster mass of 6000 M_⊙, the Plummer scale radius is a_p = 1.95 pc, yielding r_hm = 2.54 pc for λ₀ = 0.09 and a_p = 0.97 pc with r_hm = 1.27 pc for λ₀ = 0.045. Figure 3 illustrates the initial OC density distribution for λ₀ = 0.09.

As a final step, we placed our theoretical OC models in an external Galactic field at positions of r_apo = 10 and 3 kpc from the Galactic centre. In addition to the planar models, we ran the full set of models with initial positions at Z = 200 and 60 pc above the Galactic plane (see Sect. 2.2). The complete set of models and values of our SFEs is presented in Table 2.

For numerical integration, we scaled our dynamical N-body modelling using the so-called NB units (also known as Hénon units Hénon 1971; Heggie & Mathieu 1986). The chosen NB units are as follows: $$\begin{array}{l}{M}_{\text{NB}}=222.2M_{\odot },\hfill \\ V_{\text{NB}}=0.977\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1},\hfill \end{array}\begin{array}{l}{R}_{\text{NB}}=1.0\text{\hspace{0.17em}pc},\hfill \\ T_{\text{NB}}=1.0\text{\hspace{0.17em}Myr}.\hfill \end{array}$$

Fig. 3 Initial density distribution of OC for λ0 = 0.09. The blue circle indicates the half-mass radius of the cluster.

2.5 Integration procedure

For the dynamical orbital integration of OCs, including the effects of stellar evolution, we employed the high-order parallel N-body code φ-GPU¹, which is based on a fourth-order Hermite integration scheme with hierarchical individual block time steps (for details, see Berczik et al. 2011 and Berczik et al. 2013). This code is well tested and has been successfully used to obtain significant results in our prior star-cluster simulations (Shukirgaliyev et al. 2017, 2018, 2021; Ishchenko et al. 2024).

We integrated our OCs up to their dissolution time, with a maximum integration time of 4 Gyr. We define the dissolution of a cluster as the point where the remaining tidal mass decreases to less than 10% of the initial mass. This criterion is used in conjunction with the central density of the cluster, with dissolution assumed when the central density falls below 0.06 M_⊙/pc³ (Gregersen 2010).

3 Results

3.1 Mass-loss evolution

We followed the detailed dynamical and stellar evolution of mass loss in our set of OCs. In Fig. 4, we showed the evolution of the OC tidal mass M_tid as a function of time. This figure presents the two main orbital families with r_apo = 10 and 3 kpc. For each apocentre value, we generated the initial conditions with three different r_per. In Panels 1 and 4, we show the results for circular orbits with e = 0. In Panels 2, 3, 5, and 6, we show the OC models with elliptical orbits of e = 0.33 and 0.5. The colour palette, from black to light purple, represents the initial λ₀ = 0.09 values for different SFEs. The solid lines represent the OC models with initial positions exactly in the X–Y Galactic plane. Dashed lines represent the models with initial Z values of 0.2 kpc (Panels 1, 2, 3) and 0.06 kpc (Panels 4, 5, and 6) above the Galactic X – Y plane.

For easier comparison and representation of the OCs with different types of orbits, we present time on the X-axis in units of orbital periods. For example, for circular orbits with r_apo = 10 kpc, the full period is ~270 Myr, but for orbits with r_apo = 3 kpc, it is only ~80 Myr. Thus, the real-time evolution of the clusters’ mass losses differs significantly. However, by presenting time in orbital periods, we can easily compare the mass losses of clusters on significantly different orbits. The real evolutionary time of the clusters is also shown at the top of each panel.

Comparing all six plots, we can conclude that the models that come close to the Galactic centre have shorter dissolution times. This process is almost independent of their initial SFE values. This can be explained by the influence of the relatively high mass bulge in our FIX potential model. OCs with r_apo = 3 kpc and e = 0.33 and 0.5 are disintegrated after ∼3–6 rotation periods, as shown in Panels 5 and 6, depending on their initial SFE values. As we can see from Fig. 4, models with low SFE (0.15), regardless of the r_apo and e values, are destroyed much more quickly – only after several million years of evolution (less than 0.1–0.2 periods) – as shown in all panels.

The longest lived models are those with orbits of r_apo = 10 kpc and with e = 0 and 0.33, where we see smooth mass loss for almost all SFEs; this is shown in Panels 1 and 2. We also performed additional simulations with a higher SFE = 0.35, which makes the cluster more concentrated, as shown by the black lines in Panels 1, 2, 3. As expected, in this case the cluster survival rate increases for all orbital eccentricities.

We also ran a few simulations of our basic models (Panel 1) with the reduced initial λ₀ = 0.045 parameter. These models are shown in cyan. As we can see, for SFE = 0.20, there is no difference between λ₀ = 0.045 and 0.09. However, for SFE = 0.25 and 0.35, the OCs with smaller λ0 demonstrate significantly longer lifetimes.

In Panels 3 and 6, we observe more significant peaks in M_tid values for all different SFEs during the first few periods. This behaviour can be understood by considering the tidal radius evolution during the OC’s eccentric orbit. As we know, the OC tidal radius, in a first approximation, is linearly proportional to the OC’s distance from the Galactic centre. So, during the closest approach to the Galactic centre, the OC tidal radius shrinks significantly, and the tidal mass can even approach zero. As the OC flies by the Galactic centre, the tidal radius grows linearly again, and the corresponding tidal mass increases significantly.

In addition to the simple X−Y Galactic planar models, we also ran a set of simulations with the initial position of the cluster above the plane at Z = 0.2 kpc for r_apo = 10 kpc and at 0.06 kpc for r_apo = 3 kpc, including all variations of e. These models are represented by dashed lines in Fig. 4. As we can see in the plot, there are no significant differences between these two sets of runs.

In Table 3, we summarise the lifetime orbital periods of our OCs. As a formal dissolution condition, we set a 10% limit for M_tid in units of the cluster’s initial mass. From the table, we can clearly see that OCs with circular orbits at r_apo = 10 and 3 kpc are the most stable systems. For example, the lifetime of the models with SFE = 0.30 in the first case is ∼2.5 Gyr, and in the second case it is ∼0.6 Gyr. The special case of OC evolution with an artificially maximum SFE = 0.99 is presented in Fig. 9 and discussed in Sect. 5.

In Fig. 5, we present a visualisation of the data from Table 3 using a simple 3D surface. We observe a clear increasing trend in the cluster dissolution period as a function of SFE. A larger SFE (indicating a smaller remaining gas fraction in the cluster) results in a longer dissolution period. In the figure, the dependence on orbital eccentricity (e) is much less pronounced. The two sets of models with different r_apo exhibit similar behaviour, with the two surfaces nearly overlapping. In this context, our work complements the study by Cai et al. (2016), where the authors investigate the evolution of open clusters on Galactic orbits with varying orbital eccentricities. Our results align closely with their conclusions.

Fig. 4 Evolution of mass loss in OCs for different SFEs with λ0 = 0.09 (solid lines). Models with rapo = 10 kpc and e = 0, 0.33, and 0.5 are shown in the upper panels. Models with rapo = 3 kpc are shown in the bottom panels. Dashed lines represent orbits with initial Z values of 0.2 kpc upper panels) and 0.06 kpc (bottom panels) above the Galactic X − Y plane. Cyan lines represent models with λ0 = 0.045 for OCs with rapo = 10 kpc.

Fig. 5 Cluster dissolution time in orbital periods as a function of SFE and orbital eccentricity.

3.2 Evolution of the Roche volume filling factor λ

In Fig. 6, we show the evolution of the Roche volume filling factor for all our sets of models. Models with circular orbits (Panels 1 and 4) show smoother behaviour in a λ range from 0.2 to 0.4 for all SFE values, except for SFE = 0.15. Models with a minimum SFE = 0.15 fall from 0.7 to 0.45 at the beginning of the evolution before the cluster undergoes rapid final dissolution, all panels. Another prominent feature in the panels is the significant oscillation of λ as a function of time. The oscillation peaks obviously correspond to the maximum in the r_apo of the OC positions (e.g. Panel 3 in Fig. 4). The oscillations exhibit a stable period and amplitude.

In Fig. 6 (all panels), we observe no significant differences between planar (Z = 0) and non-planar models. However, we do see notable differences between the two initial λ0 sequences for all different SFEs (see cyan colour in Panel 1).

In Figs. 7 and A.1, we show the density distribution in the OCs for several models, namely for r_apo = 10 and 3 kpc with e = 0.33. For each of these models, we also present the particle distributions for different SFEs of 0.15, 0.20, 0.25, and 0.30. The evolution of the half-mass radii (marked as a blue circle) is shown for several time intervals of the OC evolution: 5, 10, 25, 50, 100 Myr, and 500 Myr of the integration. As we can see, the cores for models with SFE = 0.15 and 0.20 were almost dissolved at T = 25 and 50 Myr, respectively.

4 Cluster residual star-forming gas expulsion

To estimate the gas expulsion timescale and the residual gas after the supernova explosions in the cluster volume, we studied one example of an open cluster located inside a typical Galactic giant, cold, molecular cloud with a concentration of n_cloud ∼ 10–100 hydrogen atoms per cm³. The cluster stellar matter concentration, n_⋆ , is described by the spherically symmetric Plummer model (Plummer 1911), with a total cluster stellar mass of M_⋆ = 6000 M_⊙ and a Plummer radius of a_p = 1.95 pc. The residual star-forming gas concentration, n_gas, corresponding to a given global star formation efficiency (SFE), has been recovered using the spherically symmetric semi-analytical model from Parmentier & Pfalzner (2013), mentioned in Sect. 2. The calculated cluster stellar matter n_⋆ concentration profiles for two cases of gaseous external medium, n_cloud, are presented in Panel 1 of Fig. 8. Here, we show the distributions for two cases: SFE = 0.15 with n_cloud = 100 cm⁻³ and SFE = 0.35 with n_cloud = 10 cm⁻³.

The young hot stars of the cluster disturb the surrounding gaseous medium by wind-blowing bubbles and by creating extended regions of strong dissociation and ionisation in the remaining gas after star formation (Rahner et al. 2017, 2019; Pellegrini et al. 2020a,b). For simplicity in further calculations, we do not take into account the existence of wind-blowing bubbles around hot stars. This will result in a certain underestimation of the speed of the cluster gas expulsion process. Instead, we evaluate the phenomenon of hot stars creating the extended HII ionised hydrogen region in the cluster medium. That is, we calculate the radius of the so-called Strömgren sphere. For this purpose, we use the equation from Table 5 of the original Strömgren paper (Strömgren 1939) (for more details and discussion, see also Ritzerveld 2005). Using our standard Kroupa IMF (Kroupa 2001), we can simply estimate the number of high-mass (m ≳ 8 M_⊙) and high-temperature (T ≳ 30 000 K) main-sequence stars with O/B spectral type for our model clusters (with total stellar mass M_⋆ = 6000 M_⊙) as N_O/B ≈ 70. Using this approximate number, we can calculate the ionisation inside the gaseous medium of the cluster. Accordingly, we obtain the Strömgren sphere radius equal to r_HII = 0.9 pc for the case with SFE = 0.15 and r_HII = 2.3 pc for the case with SFE = 0.35.

As we can see, a large fraction of the gaseous medium in our model clusters mainly contains non-ionised molecular gas. The strong shockwave from the first supernova stars in the stellar system propagating throughout such an environment produces full dissociation and ionisation of the gas. This results in a noticeably higher compression ratio, more than four times that of a strong shockwave propagating through plasma, at the front of the shockwave moving with a speed of ∼100 km/s or less. The shockwave, generated by supernova explosions in the cluster, will greatly slow down when reaching the outside transition region between the gaseous medium of the cluster and the surrounding giant, cold, molecular clouds. After that, its speed may be comparable to ∼100 km/s.

Using the stellar evolution library (Kamlah et al. 2022), we calculate that (according to current stellar evolution models) out of the total number of stars (N_⋆ = 10 000 stars) in our model clusters, approximately N_SNe ≈ 90 stars will explode as supernovae of various types. The start of this process occurs at time T_SNe ≈ 4.0 Myr after the cluster formation time and stops around 60 Myr. The corresponding chronology of energy transfer by supernova explosions into the gaseous medium of the cluster, calculated under the assumption that each supernova has the same explosion energy Q_SN equal to 10⁵¹ erg, is presented in Panel 2 of Fig. 8.

The envelope ejected by the supernova explosion, moving through the supernova progenitor surrounding the gaseous medium, generates a strong shockwave. The time evolution of the shockwave’s physical parameters, such as shock-front speed, the velocity of gas behind the shock, and the spatial distribution of density and pressure, are described by the well-known Euler equations in fluid dynamics, which are a set of partial differential equations, together with the Rankine–Hugoniot conditions. Each supernova in our model, at the moment of its explosion, is located at some distance from the cluster centre. As a simplifying approximation, we neglect this fact and assume that all supernovae in our model explode near the centre of the cluster. This allows us to simplify the problem from a 3D hydrodynamic problem of a strong explosion in a gaseous medium with non-uniform density distributions to a one-dimensional problem of a strong, spherically symmetric explosion at the centre of a medium with a spherically symmetric non-homogeneous density distribution. Even such a highly idealised problem does not have an exact, complete analytical solution for the general case. However, this problem does have partial self-similar analytical solutions for the case of a gaseous medium with a power-law density distribution of ρ_gas(r) ∝ r^k, as described, for example, in Chapter 4 of Sedov (1993). This solution is an extension of the well-known Taylor–von Neumann–Sedov blast wave, which represents a strong explosion inside a uniform gaseous medium.

We adapted the mathematical expression for the cluster residual star-forming gaseous medium density distribution ρ_gas(r) to a quasi-power-law form as follows (for an arbitrarily chosen r₀ and ${\left(r-r_0\right)}^2\ll r_0^2)$: $${\rho }_{𝑔as}(r)\simeq {\rho }_{𝑔as}\left(r_0\right){\left(\frac{r}{r_0}\right)}^{-k\left(r_0\right)}.$$(4)

In the cases of gaseous medium concentration profiles shown in Panel 1, Fig. 8, the exponent k(r) increases from a value of zero at the cluster centre to approximately 2.5 at a distance of approximately 5 pc from the cluster centre. It then decreases to approximately 0.5 and lower values, reaching zero at the cluster boundary, which is approximately 20 pc.

According to Eq. (14.4) from paragraph 14 in Sedov (1993), and taking into account that the cluster gaseous medium density ρ_gas(r) ∝ r^k, the front velocity of a strong spherically symmetric shockwave formed by serial explosions of supernovae at the centre of the cluster can be written as $$V_{shock}(r)\simeq \frac{2}{5-k}\cdot {\left[Q_{{\displaystyle \sum \text{SNe}}}(r)/\left({\rho }_{gas}(r)\cdot r^3\right)\right]}^{1/2},$$(5)

where r is the radius of the shockwave front at this moment, ρ_gas(r) is the density of the gaseous medium before the front of the shockwave, and Q_ΣSNe(r) is the total energy currently deposited by supernovae explosions to the shockwave.

The calculated travel time from the centre of the cluster and the shockwave velocity are presented in Panel 4, Fig. 8. In this figure, the three horizontal coloured lines show the moments when the shockwave, formed by the explosion of the first supernova inside the cluster, is joined by the explosions of the second, third, and fourth supernovae, respectively. These later explosions also begin to deposit energy into the cluster medium.

As seen, at the moment of the explosion of the fourth supernova (~θ.9 Myr after the start of the supernovae explosions), the radius of the front of the total shock wave from the explosions of supernovae significantly exceeds the cluster boundary of ~20 pc.

Since, at this time, the shockwave is moving through the environment of the surrounding molecular cloud, the density distribution behind the front of the shockwave is close to the distribution described by the classical Taylor–von Neumann–Sedov blast-wave solution. Accordingly, almost all the gaseous mass inside the shockwave front is concentrated at a relatively small distance behind the front, towards the centre of the shockwave. Therefore, within the radius that encloses 90% of the cluster stellar mass, which is almost an order of magnitude smaller than the radius of the shockwave front, only a few percent of the total initial gaseous mass of the cluster remains inside it. In our further dynamical modelling, we neglect the gravitational influence of this remaining gaseous mass on the cluster stars.

In general, our calculations show that the first few explosions of the cluster supernovae are sufficient to produce gas expulsion inside the whole cluster. We can conclude that after ~0.9 Myr from the start of the supernovae explosions, only a few percent of the total initial cluster gaseous mass remains inside the volume.

Fig. 6 Same as in Fig. 4, but for the λ evolution.

5 Comparison with the observational data

To compare our theoretical OC models with the observed clusters in the solar neighbourhood, we used observations from two different catalogues. The first catalogue, MWSC (Kharchenko et al. 2012, 2013), contains the full set of OCs within an ~1.6 kpc radius around the Sun. It is a heliocentric sample of Galactic clusters constructed using astrometry from the Catalogue of Positions and Proper Motions on the International Celestial Reference Frame (Roeser et al. 2010), in combination with near-infrared photometry from the Two Micron All Sky Survey (2MASS Skrutskie et al. 2006). The MWSC catalogue contains basic parameters for 3061 open clusters (Kharchenko et al. 2012), 202 of which were newly discovered in the MWSC survey (Schmeja et al. 2014; Scholz et al. 2015).

The second, more recent catalogue (Hunt & Reffert 2021, 2023, 2024) contains improved data for all-sky MW OCs based on Gaia DR3 observations of more than 729 million sources. The new catalogue contains more than 5647 clusters, with the completeness limit at approximately 1.8 kpc around the Sun. The most important differences between the two catalogues are that the Gaia catalogue has a more strongly cut high-quality sample of the 3560 OC objects. These differences are also clearly visible in the range of λ parameters of OCs.

In Fig. 9, we show a comparison of the model λ for different e values and at different moments of time (ranging from 0.008 up to 0.7). This shape parameter λ increases during the evolution, similarly to the result of Ernst et al. (2015), and independently of the type of orbit. We observe a clear effect for models with different SFEs. The lower SFE values (0.15 and 0.17) demonstrate faster evolution, quickly growing up to 0.7, which is the upper limit for the observed clusters (Gaia DR3). The models with higher SFEs show moderate growth up to 0.4–0.6. Generally, from our set of SC models, we see good coverage of the observed λ range in both catalogues. It should be noted that only the recent Gaia DR3 catalogue presents OCs with λ > 0.4. These observations are in very good agreement with our eccentric orbital models. Models with circular orbits can also reach such high values of λ, but only in the case of low SFE (<0.17).

We also carried out a special set of simulations with the theoretically maximum SFE = 0.99. In this set of simulations, we mimicked the previous (older) gas-free cluster initial conditions. As we can see, these models are compatible with the observations if we assume Galactic eccentric orbits for our clusters.

We also checked the influence of the initial λ₀ = 0.045 on the λ time evolution of the models with circular orbits for r_apo = 10 kpc (black symbols). The initial λ₀ has a strong influence on the cluster’s filling factor and, as a consequence, on the OC’s survival time.

In Fig. 10, we present the cluster filling factor distribution for a different set of SFE models, including the special models with SFE = 0.35. Generally, we can conclude that models with lower SFE obviously cannot reproduce the observed cluster distributions. With higher SFE, we managed to reproduce the general trends in both catalogues’ distributions. Especially good coverage was achieved for the models with SFE = 0.30 and 0.35.

Our results confirm the general idea that different clusters, at least in the solar vicinity, formed within a given SFE range. Only in this case can we reproduce both the lower λ range of the distribution and the higher λ range. Small SFE systems can reproduce the higher λ range (0.6–0.8), while larger SFEs can reproduce the lower λ range (0.01–0.2). Models with non-planar orbits show additional features compared to the planar models. In the end, we observe that the combination of our different model sets can reproduce the wide range of observed distributions. Additionally, it should be mentioned that the number of model clusters in the theoretical curves on our histograms are used directly, without any extra re-normalisation. The good correspondence between observed and theoretical numbers can be interpreted as additional proof of our concept of different SFE formation for OCs in the solar vicinity.

In Fig. 11, we present the λ evolution of the OCs as a function of the remaining cluster stellar mass. As we already observed from the previous plots, the lower SFE models describe the whole range of observed Roche-volume filling factors in both catalogues better. The higher SFEs essentially reproduce only the MWSC observed values. We clearly see that the r_apo = 3 kpc and r_apo = 10 kpc models are well over-plotted, essentially covering the same λ range. As mentioned in the previous paragraph, to describe all the different observed filling factors, we probably need to use a weighted set of different models with varying SFEs (and possibly different initial masses), but such a comprehensive study is well beyond the scope of our current paper.

Fig. 7 Evolution of density distribution in the centre of the OC for an elliptical orbit with rapo = 10 kpc and e = 0.33. Snapshots for 5, 10, 25, 50, 100, and 500 Myr are shown. The blue circle represents the half-mass radius.

Fig. 8 Characteristic profiles of gas expulsion processes inside the star cluster (from left to right, top to bottom). Calculated concentration profiles of the cluster stellar n★ and gaseous nɡas material for the two cases: SFE = 0.15 with ncload = 100 cm−3 and SFE = 0.35 with ncload = 10 cm−3. The supernova explosion energy rate inside the cluster as a function of time is shown. The shockwave velocity and the shockwave travel time from the centre of the cluster as a function of the front radius. The dotted and dashed-dotted black vertical lines represent the radii within which 50% and 90% of the cluster stellar mass are enclosed, respectively. The dashed black vertical line represents the simulated cluster boundary, rcluster, which we assume to be equal to 20 ap (for the example model, it is ~19.5 pc).

Fig. 9 Evolution (vertical points) of clusters’ Roche-volume filling factors − λ as a function of orbital eccentricities e = 0, 0.33, 0.5 for the two values of orbital apocentres rapo = 10 and 3 kpc. With different coloured circles we represent the different values of global SFE. All the basic runs start with the default λ0 = 0.09. Additionally, the set of three runs with the two-times smaller λ0 = 0.045 presented with the crosses. Filled circles represent the planar orbits (Z = 0), while non-filled circles represent the non-planar orbits. The range of the observed λ based on the OC data from Gaia DR3 is shown in light red, and the MWSC catalogue is shown in light blue.

Fig. 10 Number of clusters (inside the catalogue’s completeness limits around the Sun) as a function of the observed cluster’s Roche-volume filling factor − λ. The distribution of the observational OC data from Gaia DR3 is shown in light red, and the MWSC catalogue is shown in light blue. The different panels show the theoretical distributions for different SFE models (with rapo = 10 kpc). The lines with different colours show models with different orbital eccentricities. The dotted lines present the results for the non-planar orbits (Z0 = 200 pc). The cyan lines show the results for the three special models with the half-reduced initial λ0 = 0.045.

Fig. 11 Evolution of half-mass Roche-volume filling factor λ as a function of the bound mass M/M⊙ for SFE = 0.15, 0.17, 0.20, 0.25,0.30, and 0.35 with λ0 = 0.09. The blue coloured palette curves represent rapo = 10 kpc with e = 0, 0.5, and 0.33; while the red coloured palette curves represent the same, but for rapo = 3 kpc. Solid curves represent the data for Z = 0, and dashed curves represent Z = 0.2 kpc for rapo = 10 kpc and Z = 0.06 kpc for rapo = 3 kpc. The cyan curves represent data for λ0 = 0.045 for rapo = 10 and e = 0. The distribution of the observational OC data from Gaia DR3 is shown in light red, and that from MWSC is shown in light blue.

6 Conclusions and discussions

In this study, we considered the survivability of open clusters with different global star formation efficiency models and Plummer density profiles after instantaneous gas expulsion (Shukirgaliyev et al. 2017, 2018, 2019). Specifically, we considered model clusters with a fixed initial stellar mass of 6000 M_⊙ – corresponding to the Plummer model with λ₀ = 0.09 – with an additional 0.045 in virial equilibrium within the residual gas immediately before the instantaneous gas expulsion.

Our calculations show that a few initial supernova explosions in the cluster are sufficient to expel gas far beyond the cluster boundary (≈20 pc). Therefore, within the radius that encloses 90% of the cluster stellar mass – almost an order of magnitude smaller than the radius of the shockwave front at ≈0.9 Myr after the onset of supernova explosions - only a few percent of the total initial gaseous mass remains inside the system.

The density profiles of the residual star-forming gas corresponding to a given global SFE were recovered assuming that stars formed with a constant efficiency per free-fall time (Parmentier & Pfalzner 2013; Shukirgaliyev et al. 2017). We then performed direct N-body simulations of cluster evolution after instantaneous gas expulsion, assuming they are orbiting in circular and eccentric trajectories in the Galactic disc, with two fixed apocentres at 10 and 3 kpc.

The theoretical models for the Roche-volume filling factor distribution and time evolution coincide well with the observed Gaia DR3 and MWSC catalogues. Even our restricted set of models shows good coverage of the clusters’ filling parameters. The set of non-planar, orbiting, open-cluster evolution shows a very similar mass loss and filling parameter evolution compared to the planar orbiting star clusters. Compared to previous models (Ernst et al. 2010, 2015), where the authors used 100% star-forming clusters, our different SFE initial condition models can easily cover the observed cluster filling factor ranges, even up to λ = 0.7. To describe all the observed filling factors, we definitely need to use some form of weighted set of theoretical models with different SFEs and initial masses.

As a general conclusion from our set of different SFE models, we observe a relatively short lifespan for typical MW disk OCs that is usually around ten or fewer orbital periods. We note that OCs on eccentric orbits have only slightly shorter lifespans. Our results for the cluster dissolution time from the eccentric orbital simulations show the same behaviour as the models presented in Cai et al. (2016).

The MW OC dissolution models presented in this paper can serve as a basis for future open-cluster population modelling in the Galaxy. The main feature of our models is the wide range of different star formation efficiency prescriptions. This SFE parameter highlights the strong importance of gas dynamical processes and gas expulsion efficiency for young OCs.

Acknowledgements

We thank the referee for the helpful comments. This research has been funded by the Committee of Science of the Ministry of Science and Higher Education of the Republic of Kazakhstan (Grant Nos. BR24992759, AP19677351). MI and PB thanks the support from the special program of the Polish Academy of Sciences and the U.S. National Academy of Sciences under the Long-term program to support Ukrainian research teams grant No. PAN.BFB.S.BWZ.329.022.2023. PB acknowledges the support by the National Science Foundation of China under grant NSFC No. 12473017. BS acknowledges the support provided through the Nazarbayev University Faculty Development Competitive Research Grant Program, with Grant No. 11022021FD2912. BS acknowledges the use of free version of ChatGPT (OpenAI) for language polishing and manuscript text refinement.

Appendix A Evolution of the density distribution in the centre of the OC for model with r_apo = 3 kpc and e = 0.33

Fig. A.1 Same as Fig. 7, but for model with rapo = 3 kpc and e = 0.33.

References

All Tables

All Figures

	Fig. 1 Orbital evolution of OC with e = 1/3 is shown in the X-Y and R- Z planes (where R is the planar galactocentric radius). The integration time is 1 Gyr, as indicated by the colour bar. Upper panel: orbit with rapo = 10 kpc (Z = 200 pc). Lower panel: orbit with rapo = 3 kpc (Z = 60 pc).
In the text

Fig. 2 Concentration profiles of cluster’s stellar (n⋆) and gaseous (n𝑔as) material are shown for different SFEs, ranging from 0.15 to 0.35 and 0.99. The dotted and dash-dotted black vertical lines represent the radii enclosing 50% and 90% of the cluster’s stellar mass, respectively. The dashed black vertical line marks the simulated cluster boundary, rcluster, which is defined as 20 ⋅ ap. For the rapo = 10 kpc and λ0 = 0.09 models (see Sect. 2.4), this corresponds to ~19.5 pc.

In the text

	Fig. 3 Initial density distribution of OC for λ0 = 0.09. The blue circle indicates the half-mass radius of the cluster.
In the text

Fig. 4 Evolution of mass loss in OCs for different SFEs with λ0 = 0.09 (solid lines). Models with rapo = 10 kpc and e = 0, 0.33, and 0.5 are shown in the upper panels. Models with rapo = 3 kpc are shown in the bottom panels. Dashed lines represent orbits with initial Z values of 0.2 kpc upper panels) and 0.06 kpc (bottom panels) above the Galactic X − Y plane. Cyan lines represent models with λ0 = 0.045 for OCs with rapo = 10 kpc.

In the text

	Fig. 5 Cluster dissolution time in orbital periods as a function of SFE and orbital eccentricity.
In the text

	Fig. 6 Same as in Fig. 4, but for the λ evolution.
In the text

	Fig. 7 Evolution of density distribution in the centre of the OC for an elliptical orbit with rapo = 10 kpc and e = 0.33. Snapshots for 5, 10, 25, 50, 100, and 500 Myr are shown. The blue circle represents the half-mass radius.
In the text

Fig. 8 Characteristic profiles of gas expulsion processes inside the star cluster (from left to right, top to bottom). Calculated concentration profiles of the cluster stellar n★ and gaseous nɡas material for the two cases: SFE = 0.15 with ncload = 100 cm−3 and SFE = 0.35 with ncload = 10 cm−3. The supernova explosion energy rate inside the cluster as a function of time is shown. The shockwave velocity and the shockwave travel time from the centre of the cluster as a function of the front radius. The dotted and dashed-dotted black vertical lines represent the radii within which 50% and 90% of the cluster stellar mass are enclosed, respectively. The dashed black vertical line represents the simulated cluster boundary, rcluster, which we assume to be equal to 20 ap (for the example model, it is ~19.5 pc).

In the text

Fig. 9 Evolution (vertical points) of clusters’ Roche-volume filling factors − λ as a function of orbital eccentricities e = 0, 0.33, 0.5 for the two values of orbital apocentres rapo = 10 and 3 kpc. With different coloured circles we represent the different values of global SFE. All the basic runs start with the default λ0 = 0.09. Additionally, the set of three runs with the two-times smaller λ0 = 0.045 presented with the crosses. Filled circles represent the planar orbits (Z = 0), while non-filled circles represent the non-planar orbits. The range of the observed λ based on the OC data from Gaia DR3 is shown in light red, and the MWSC catalogue is shown in light blue.

In the text

Fig. 10 Number of clusters (inside the catalogue’s completeness limits around the Sun) as a function of the observed cluster’s Roche-volume filling factor − λ. The distribution of the observational OC data from Gaia DR3 is shown in light red, and the MWSC catalogue is shown in light blue. The different panels show the theoretical distributions for different SFE models (with rapo = 10 kpc). The lines with different colours show models with different orbital eccentricities. The dotted lines present the results for the non-planar orbits (Z0 = 200 pc). The cyan lines show the results for the three special models with the half-reduced initial λ0 = 0.045.

In the text

Fig. 11 Evolution of half-mass Roche-volume filling factor λ as a function of the bound mass M/M⊙ for SFE = 0.15, 0.17, 0.20, 0.25,0.30, and 0.35 with λ0 = 0.09. The blue coloured palette curves represent rapo = 10 kpc with e = 0, 0.5, and 0.33; while the red coloured palette curves represent the same, but for rapo = 3 kpc. Solid curves represent the data for Z = 0, and dashed curves represent Z = 0.2 kpc for rapo = 10 kpc and Z = 0.06 kpc for rapo = 3 kpc. The cyan curves represent data for λ0 = 0.045 for rapo = 10 and e = 0. The distribution of the observational OC data from Gaia DR3 is shown in light red, and that from MWSC is shown in light blue.

In the text

	Fig. A.1 Same as Fig. 7, but for model with rapo = 3 kpc and e = 0.33.
In the text