© ESO 2022

1. Introduction

Classical Cepheids occurring in star clusters (henceforth clustered Cepheids) are doubly important objects for stellar physics and cosmology (via the distance scale) because both their pulsations and their host cluster’s member stars provide stringent tests of stellar models and allow us to calibrate the Leavitt law (Leavitt & Pickering 1912, LL, period-luminosity relation). The identification of Cepheids in clusters was initially a one-by-one procedure (Irwin 1955; Kholopov 1956; Sandage 1958) and required painstaking reddening corrections (e.g. Turner 1976; Turner & Burke 2002, and references therein) before large amounts of high-quality astrometric data became available. Anderson et al. (2013) adopted a new strategy for systematically and quantitatively establishing the membership of Cepheids in Milky Way (MW) open star clusters using eight-dimensional phase-space information with an eye on the Gaia mission launched in December 2013 (Gaia Collaboration 2016). Gaia’s unprecedented astrometric census greatly facilitates the identification of MW star clusters based chiefly on proper motion clustering (e.g. Cantat-Gaudin et al. 2018; Castro-Ginard et al. 2020), while ground- and space-based time-domain surveys (e.g. Pawlak et al. 2019; Soszyński et al. 2020; Holl et al. 2018; Ripepi et al. 2019) considerably increased the number of MW Cepheid candidates. This combination has led to reassessments of cluster membership as well as several new interesting cluster Cepheid candidates (Medina et al. 2021; Zhou & Chen 2021). However, a main limitation of such modern astro-statistical approaches is the formulation of prior probabilities based on (2D) on-sky locations. Understanding how and how frequently Cepheids occur in star clusters or their halos would significantly help to improve the reliability of such priors.

The overall empirical result of such studies has been that (a) the vast majority of Cepheids do not reside in clusters and the clustered Cepheid fraction f_CC, true ≲ 10% (Anderson & Riess 2018; b) MW clusters that do host Cepheids generally contain a single, occasionally two (e.g. Dékány et al. 2015), and at the very most three Cepheids (Majaess et al. 2013; Anderson et al. 2013; Medina et al. 2021, NGC 7790); (c) a sizeable fraction of reported MW cluster Cepheids are located in the distant outskirts (halos) of the putative host clusters (Anderson et al. 2013; Zhou & Chen 2021; Medina et al. 2021); and (d) f_CC, true is generally higher for the younger long-period Cepheids (Anderson & Riess 2018). The most significant exceptions to this rule are very massive (≈5 × 10⁴ M_⊙) clusters in the Large Magellanic Cloud (LMC), such as NGC 1866 or NGC 2031, that each contain up to 24 Cepheids, most of them with a fairly short period (Testa et al. 2007; Musella et al. 2016).

The identification of clustered Cepheids becomes increasingly difficult at distances beyond 100 kpc both because the detection of the main sequence in the clusters becomes more difficult and because the angular resolution element of a detector corresponds to increasingly larger physical scales. Blending Cepheids with cluster member stars depresses observed photometric amplitudes and alters light curve shapes, thus precluding the detection of Cepheids occurring in massive clusters (Ferrarese et al. 2000; Riess et al. 2020). This results in a lower observed effective clustered extragalactic Cepheid fraction, f_CC, eff, compared to the true number of clustered extragalactic Cepheids, and any uncertainties in f_CC, eff affect stellar association bias corrections for the distance ladder linearly (Anderson & Riess 2018). Understanding the difference between f_CC, eff and f_CC, true and possible variations in f_CC among different galaxies will be valuable for improving distance measurements based on Cepheids (Anderson & Riess 2018). Providing a better astrophysical understanding of how Cepheids relate to their (former) host clusters further aids the calibration of the LL, which can be significantly improved by adopting precise average cluster parallaxes for LL calibration (Breuval et al. 2020).

The occurrence of clustered Cepheids is intimately related to star cluster dynamics because most stars are thought to be born in young embedded clusters (Lada et al. 1991; Kroupa 1995a; Lada & Lada 2003; Porras et al. 2003; Megeath et al. 2016) that evolve dynamically over time, maintaining some imprints of the initial dynamical interactions within their dense birth environments. Following star cluster formation, these objects relax and mass segregate (e.g. Spitzer 1969; Bonnell & Davies 1998; Mouri & Taniguchi 2002; McMillan et al. 2007; Spera et al. 2016), eject stars either in close interactions of three or more stars (e.g. Aarseth 1971; Heggie 1975; Perets & Šubr 2012; Tanikawa et al. 2012; Oh et al. 2015) or as a result of a decay of primordial triples and quadruples (e.g. Sterzik & Durisen 1998; Goodwin & Kroupa 2005; Pflamm-Altenburg & Kroupa 2006; Wang et al. 2019), and evaporate stars as the result of weak encounters between stars (e.g. Küpper et al. 2008, 2010). The loss of stars from clusters is enhanced by the tidal field of the cluster hosting galaxy (Baumgardt et al. 2002; Baumgardt & Makino 2003; Sollima 2020) and by removal of the residual gas which was not consumed in the star forming process (e.g. Lada et al. 1984; Kroupa et al. 2001; Baumgardt & Kroupa 2007). The remaining fraction of stars that resides within their birth clusters then decreases as a function of age (assuming a constant star formation rate), as is observed for clustered Cepheids in the MW, LMC, Small Magellanic Cloud (SMC), and M 31 (Anderson & Riess 2018, their Fig. 13). Such a dynamical picture also explains the occurrence and aids the discovery of extended tidal tails of clusters dissolving into the MW field population (e.g. Röser et al. 2019; Meingast & Alves 2019; Jerabkova et al. 2021). It may also help to confront predicted period-age relations (Bono et al. 2005; Anderson et al. 2016) to observations (Senchyna et al. 2015; Medina et al. 2021) because the dynamical timescale is independent of the assumptions or input physics of stellar models. This is particularly interesting, since Cepheid period-age relations are important to trace Galactic evolution even though the ages of Cepheids depend by nearly a factor of two on the effects of main-sequence rotation (Anderson et al. 2016).

To address the above issues, we computed an extensive library of N-body cluster simulations to investigate cluster dispersal over a timescale of up to 300 Myr. In this first paper of a series exploiting these simulations, we present the computed library and its predictions concerning the occurrence of Cepheids in star clusters for different cluster masses, chemical composition, and galactocentric radius. Of course, cluster dynamics also impacts the orbital parameters of primordial binary stars (Kroupa 1995b; Oh & Kroupa 2016), and Cepheids are particularly likely to have experienced binary interactions, since the close orbital configurations observed among their progenitors (B-type main sequence stars) must adapt to the significant evolutionary increase in radius (typical Cepheid radii 20 − 200 R_⊙). Elucidating the nature and frequency of these binary interactions is crucial for the forward modelling of populations of Cepheids and for understanding the intriguing hierarchy of multiple star systems hosting Cepheids (Evans et al. 2020). Paper II of this series will focus on such effects (Dinnbier et al., in prep.).

This article is structured as follows. Section 2 describes the set up of the simulations, initial conditions, and other choices made to create the cluster library. Section 3 presents how Cepheids are identified and traced in the NBODY6 simulations. Section 4 investigates the evolution and occurrence of Cepheids in dynamically evolving star clusters as a function of the cluster mass, metallicity, and galactocentric radius. Section 5 seeks to estimate the number of Cepheids occurring in clusters, f_CC, true, and establishes upper and lower limits based on the N-body simulations. The discussion in Sect. 6 focuses on the impact of gas expulsion, initial conditions, and the role of molecular clouds. The final section, Sect. 7 summarises our results and conclusions.

2. Numerical methods and initial conditions

2.1. Numerics

The cluster simulations presented in this work are calculated by use of the Aarseth code NBODY6 (Aarseth 1999, 2003). The code integrates stellar trajectories with a fourth-order Hermite predictor-corrector scheme (Makino 1991) with individual block time steps coupled to the Ahmad-Cohen method (Ahmad & Cohen 1973; Makino & Aarseth 1992). Stellar evolution is treated by stellar evolutionary tracks due to Tout et al. (1996) and Hurley et al. (2000). To tackle the enormous range of relevant timescales, NBODY6 employs various sophisticated numerical tools such as two-body (Kustaanheimo & Stiefel 1965), three-body (Aarseth & Zare 1974), and multiple-body (Mikkola & Aarseth 1990) regularisation, which include the dynamical formation of these systems, exchanges, and ionisations. The regularisation techniques typically transform the coordinates to a four-dimensional space, and they utilise the properties of the Levi-Civita matrix (Kustaanheimo & Stiefel 1965; Bettis & Szebehely 1971; Stiefel & Scheifele 1975). The code initiates or terminates a particular regularisation technique according to given criteria. Not all compact systems are necessarily regularised. For example, unperturbed hierarchical systems are evolved as single bodies until an external perturbation or stellar evolution within the hierarchy necessitates an update of the system internal conditions (Mardling & Aarseth 1999). During regularisation, the code approximates binary star evolution with Roche-lobe mass transfer (Eggleton 1983; Hurley et al. 2002). Stellar evolution produces compact objects, and it also handles the products of stellar mergers (Kochanek 1992). We refer the reader to Aarseth (2003) for an in-depth description of relevant algorithms.

We assume that all stars are formed in embedded star clusters containing a substantial reservoir of unprocessed gas at the time of massive star formation. The gas is expelled due to feedback from massive stars, and the sudden change of the gravitational potential of the expelled gas affects the dynamics of the star cluster, which expands and loses a substantial fraction of its stars (Lada et al. 1984; Goodwin 1997; Geyer & Burkert 2001; Kroupa et al. 2001). To limit the number of free parameters, we adopted the gas expulsion model by Kroupa et al. (2001), where the gravitational potential ϕ_gas(r, t) of the gas is represented by a Plummer sphere that starts to exponentially decrease after a delay time t_d, that is,

$$\begin{array}{c}\hfill {\varphi }_{\mathrm{gas}}(r,t)=\{\begin{array}{c}-\frac{G{M}_{\mathrm{gas}}}{\sqrt{a_{\mathrm{gas}}^2+r^2}}t\le t_{\mathrm{d}},\hfill \\ -\frac{G{M}_{\mathrm{gas}}}{\sqrt{a_{\mathrm{gas}}^2+r^2}}exp\{-(t-t_{\mathrm{d}})/{\tau }_{\mathrm{M}}\}t>t_{\mathrm{d}}.\hfill \end{array}\end{array}$$(1)

In the gas model, G, M_gas, a_gas, and τ_M denote the gravitational constant, initial mass of the gas component of the embedded star cluster, the Plummer length parameter of the gas component, and the gas expulsion timescale, respectively. For the parameters of the gas expulsion model, we adopted the values that seem to be the most promising according to current understanding of embedded star clusters (Kroupa et al. 2001; Lada & Lada 2003; Kuhn et al. 2014; Banerjee & Kroupa 2017). We adopted a_gas = a_ecl, where a_ecl is the Plummer length parameter of stars, M_gas = 2 M_ecl, where M_ecl is the initial stellar mass of the embedded cluster (this implies the star formation efficiency SFE ≡ M_ecl/(M_ecl + M_gas) = 1/3), τ_M = a_gas/c_II ≈ 0.03 Myr, where c_II = 10 km s⁻¹ is the sound speed in ionised hydrogen, and t_d = 0.6 Myr is approximately the lifetime of the ultra-compact HII region in its deeply embedded phase.

The simulated clusters orbited a Milky Way-like galaxy on circular trajectories. The escaping stars were integrated in the gravitational field of the galaxy, so they form tidal tails, and stellar evolution was followed also in the tidal tails. The gravitational potential of the galaxy was represented by the model of Allen & Santillan (1991), which consists of three components: the bulge modelled as a Plummer potential of mass 1.4 × 10¹⁰ M_⊙ and scale length 0.38 kpc; the stellar disc modelled as a Miyamoto-Nagai (Miyamoto & Nagai 1975) potential of mass 8.5 × 10¹⁰ M_⊙ and radial and vertical scale lengths 5.3 kpc and 0.25 kpc, respectively; and a dark matter halo of a two-component potential that is nearly logarithmic at large galactocentric distances, and which encompasses a mass of 8 × 10¹¹ M_⊙ within a sphere of 100 kpc radius. The particular form of the galactic potential has likely negligible influence on the results, as we discuss in Sect. 6.3.

2.2. Initial conditions

The clusters were generated by the software package MCLUSTER (Küpper et al. 2011) as Plummer models of stellar mass M_ecl and Plummer parameter a_ecl. Employing the relation between the embedded cluster radius and its mass in stars (Marks & Kroupa 2012),

$$\begin{array}{c}\hfill \frac{a_{\mathrm{ecl}}}{\mathrm{pc}}=0.077{\left(\frac{M_{\mathrm{ecl}}}{M_{\odot }}\right)}^{0.13},\end{array}$$(2)

we reduced the two-parameter space (M_ecl, a_ecl) to a one-parameter family of models characterised by cluster mass M_ecl. The clusters were initially in virial equilibrium with the primordial gas. The stars were populated to the cluster according to the recipe of Aarseth et al. (1974), where the position is generated first from the cluster mass profile, whereupon the absolute value of the velocity is drawn from the velocity distribution at the given cluster radius from the Plummer distribution function, and the direction of the velocity vector is oriented randomly. The clusters were initially not mass segregated.

The initial masses of stars were sampled from the model of the stellar initial mass function (IMF) given by Kroupa (2001), where the maximum allowed stellar mass, m_max, within the cluster is a function of M_ecl according to the m_max − M_ecl relation (Weidner et al. 2010). Thus, only massive clusters harboured the most massive stars. The minimum stellar mass was 0.08 M_⊙. All stars formed in binaries with orbital parameters for stars of < 3 M_⊙ according to the initial distribution function for periods, eccentricities and mass ratios of Kroupa (1995b), and more massive stars (> 3 M_⊙) according to the binary properties of O stars (Sana et al. 2012; Moe & Di Stefano 2017). The method describing how to generate a realistic initial binary population without affecting the IMF is due to Oh et al. (2015) and Oh & Kroupa (2016). Since the evolution of the binaries is the focus of paper II, the details of the initial binary population and the evolution processes are provided there.

2.3. Assumptions about the embedded cluster mass function

There is observational evidence that the embedded cluster mass function (ECMF) in the Galaxy (e.g. Lada & Lada 2003; de la Fuente Marcos & de la Fuente Marcos 2004) as well as in more distant galaxies (e.g. Whitmore et al. 1999; Bik et al. 2003; Gieles et al. 2006) can be approximated by a power law of slope −2, that is,

$$\begin{array}{c}\hfill {\xi }_{\mathrm{ecl}}(M_{\mathrm{ecl}})\equiv \frac{\mathrm{d}{N}_{\mathrm{ecl}}(M_{\mathrm{ecl}})}{\mathrm{d}{M}_{\mathrm{ecl}}}\propto M_{\mathrm{ecl}}^{-2},\end{array}$$(3)

where dN_ecl(M_ecl) is the number of embedded star clusters forming in a mass bin of size dM_ecl.

Our models span a cluster mass range of M_ecl ∈ (50 M_⊙, 6400 M_⊙), where the lower mass limit follows from the m_max − M_ecl relation, which requires m_max ≳ m_min, Ceph, and the upper mass limit is comparable to the mass of the most massive young star clusters observed in the galaxy M 31 (≈8.5 × 10³ M_⊙; Johnson et al. 2017). The dependence ${\xi }_{\mathrm{ecl}}(M_{\mathrm{ecl}})\propto M_{\mathrm{ecl}}^{-2}$ has the convenient property that the total mass of clusters in each logarithmic mass bin of constant size contains an identical total mass of clusters regardless of M_ecl.

2.4. The library of simulations

The cluster masses, M_ecl, span the range of 50 M_⊙ to 6400 M_⊙, and they are sampled by multiples of 2; that is M_ecl = 50 M_⊙, 100 M_⊙, …, 6400 M_⊙. The models are listed in Table 1. We ran each model Nmod times with different random seeds for stellar masses, positions, and velocities, where Nmod = 2ⁱ and i is the index after the model name; for instance, models M4 have Nmod = 16. We chose these values of Nmod to have the same total mass of star clusters in each logarithmic cluster mass bin dM_ecl when used with the ECMF of Eq. (3). This facilitates the calculation of any quantity ⟨q⟩ for the whole population of clusters (and thus its observed value integrated over a galaxy) as $⟨q⟩=(\sum _{i=1}^8{q}_{\mathrm{i}}{w}_i)/(\sum _{i=1}^8{w}_i)$, where q_i is the quantity of interest and w_i the relevant weighting factor (e.g. number of all Cepheids originating in all clusters of index i¹).

We calculated the population of clusters described in the previous paragraph for three galactocentric radii: R_g = 4 kpc, R_g = 8 kpc, and R_g = 12 kpc for solar metallicity; that is, Z = 0.014. We also calculated the same population of clusters for a galactocentric radius of R_g = 8 kpc and sub-solar metallicities of Z = 0.006 and Z = 0.002. Thus, for each considered choice of R_g and Z, we calculated 510 simulations in total.

3. Identification of Cepheids

Isolated stars with zero age main sequence (ZAMS) mass in the appropriate mass range (m_min, Ceph, m_max, Ceph) enter the Cepheid instability strip in the course of their evolution. Present models take into account the metallicity dependence of both the stellar evolutionary tracks (Hurley et al. 2000) and the borders of the instability strip for Cepheid stars (Anderson et al. 2016). They provide the mass range for Cepheids as follows: (m_min, Ceph, m_max, Ceph) = (3.2 M_⊙, 11.0 M_⊙) for Z = 0.002, (m_min, Ceph, m_max, Ceph) = (4.0 M_⊙, 12.0 M_⊙) for Z = 0.006, and (m_min, Ceph, m_max, Ceph) = (4.7 M_⊙, 13.0 M_⊙) for Z = 0.014. For the purpose of this work, we term all stars with initial mass in the range of (m_min, Ceph, m_max, Ceph) prospective Cepheids (ProCephs), and stars with initial masses below m_min, Ceph as lower mass stars. The star is then referred to as ProCep or a lower mass star for the rest of the simulation regardless of its actual mass, which can change with time. If the ProCep becomes a Cepheid, it is denoted a canonical Cepheid (CanCeph). There are also Cepheids that originate from a coalescence in a binary with both components having a ZAMS mass below m_min, Ceph. We call these upscattered Cepheids (UpCeps). For the adopted evolutionary tracks, stars with mass m_min, Ceph terminate core Helium burning at time t_He, end of 300 Myr, 190 Myr and 146 Myr for Z = 0.002, Z = 0.006, and Z = 0.014, respectively.

We searched for Cepheids automatically in the HR diagrams produced by NBODY6. For each star in the simulation, we recorded its temperature and luminosity whenever its stellar evolution was updated, which, during core helium burning, typically occurs every 0.01 to 0.3 Myr depending on the Cepheid mass and exact position in the HR diagram. This way, we obtained densely sampled HR diagrams even during the rapid evolution on the instability strip as shown in Fig. 1. Then, the algorithm diagnosed each star and classified it as a Cepheid if the star fulfilled all of the following criteria. (i) The star’s initial mass is in the mass range (m_min, Ceph/2, m_max, Ceph). The lower value of the mass range is to include a possible merger of two 0.5 m_min, Ceph stars in a hard binary, which can evolve to a Cepheid (UpCeps). (ii) The star undergoes core Helium burning. (iii) The surface temperature at which core Helium burning starts is rightwards from the blue boundary of the instability strip to exclude the stars which do not have Cepheid loops (stars with m > m_max, Ceph typically have a monotonically decreasing temperature during core Helium burning). (iv) The maximum temperature during core helium burning is leftwards from the red boundary of the instability strip, which means that the star lies on the instability strip during its core Helium burning at least on one output file. We used metallicity-dependent boundaries of the instability strip from Anderson et al. (2016). Verification of this searching procedure for more complicated HR diagrams with merging stars is presented in Appendix A.

Fig. 1. HR diagrams outputted from NBODY6 for single stars of ZAMS mass of 4 M⊙ (left panel), 7 M⊙ (middle panel), and 10 M⊙ (right panel). The evolutionary paths are plotted after the stars left the main sequence. The colour indicates the metallicity, which is Z = 0.002 (blue), Z = 0.006 (green), and Z = 0.014 (red). Large crosses denote the time of core helium burning. For the purpose of this plot, we selected single stars, which evolve well in isolation from the other stars in the cluster. Thick solid lines show the Cepheid instability strip for the metallicities considered (Anderson et al. 2016). Thin solid lines represent ZAMS for stars of corresponding metallicity.

To obtain the properties of the Cepheid (e.g. its distance from the cluster centre or binarity), we first determined the times when the Cepheid crosses the instability strip and then evaluated the desired properties at a random time (the event when the Cepheid is observed) from the time interval when the star is located within the instability strip. If the Cepheid passes the instability strip twice (we ignored the first passage through the instability strip because of its very short duration), we determined the time randomly from one of the two intervals on the instability strip weighted by their duration. The dynamics of Cepheids usually changes very little when it is within the strip because the duration in the strip is significantly shorter than the cluster relaxation time.

We classify the Cepheid as a cluster Cepheid if it is located closer than 10 pc to the density centre (Casertano & Hut 1985) of the cluster. The majority of Cepheids bound to the cluster are located near the cluster centre (closer than ≈2 pc) as a result of the cluster half-mass radii (several parsec) and mass segregation of ProCeps. The distance of 10 pc is comparable to the tidal radius of the clusters, so if a star exceeds this distance, it is likely to become unbound and recede. If the cluster is already dissolved, all its Cepheids are classified as field Cepheids.

To measure the fraction of Cepheids located in clusters (relatively to the Cepheids in the field and in clusters combined) in different contexts, we introduce the following terminology. The fraction of Cepheids in clusters of initial mass M_ecl at given time t is denoted p_CC(M_ecl, t). The fraction of Cepheids in clusters of a particular cluster mass M_ecl throughout the time interval (0, 300 Myr) is denoted p_CC(M_ecl),

$$\begin{array}{c}\hfill p_{\mathrm{CC}}(M_{\mathrm{ecl}})=\frac{{\sum }_i{p}_{\mathrm{CC}}(M_{\mathrm{ecl}},t_{\mathrm{i}})\mathrm{\Delta }{t}_{\mathrm{i}}}{{\sum }_i\mathrm{\Delta }{t}_{\mathrm{i}}},\end{array}$$(4)

where Δt_i is the time interval near time t_i, and the summation goes over all time intervals. The galaxy-wide fraction of Cepheids in clusters; that is, Cepheids formed in all star clusters within the galaxy over their life-time are denoted f_CC,

$$\begin{array}{c}\hfill f_{\mathrm{CC}}=\frac{{\sum }_i{\xi }_{\mathrm{ecl}}(M_{\mathrm{ecl},i})p_{\mathrm{CC}}(M_{\mathrm{ecl},i})\mathrm{\Delta }{M}_{\mathrm{ecl},i}}{{\sum }_i{\xi }_{\mathrm{ecl}}(M_{\mathrm{ecl},i})\mathrm{\Delta }{M}_{\mathrm{ecl},i}},\end{array}$$(5)

where ξ_ecl(M_ecl) is the ECMF of Eq. (3), and the summation goes over all cluster mass intervals. We stress that throughout this work, M_ecl denotes the initial (not actual) cluster mass. As the cluster evolves, its actual mass can be significantly lower than its initial value. In addition, we study the fraction of Cepheids in clusters p_CC(m_Ceph), which are of a given stellar mass m_Ceph regardless of the mass of their birth star clusters.

To gain basic insight, it is useful to mention some typical numbers of ProCeps formed in a coeval stellar population (CSP). Consider a CSP that is populated by the Kroupa IMF from the mass interval of 0.08 M_⊙ to 50 M_⊙. For such a CSP, the total mass of ProCeps comprises 17%, 14%, and 13% of the total initial CSP mass for Z = 0.002, Z = 0.006, and Z = 0.014, respectively. By number, ProCeps constitute 1.7%, 1.2%, and 0.95% of the total initial CSP for Z = 0.002, Z = 0.006 and Z = 0.014, respectively.

4. Cepheids in dynamically evolving star clusters

In order to quantify the position of Cepheids relatively to their birth clusters, we study the influence of the dynamical evolution of star clusters on ProCeps. First, we review the most relevant features of the dynamics of young star clusters experiencing expulsion of their residual gas (Sect. 4.1). Then, we investigate the mechanisms that unbind ProCeps from clusters and therefore determine the value of f_CC as a function of cluster age and mass for clusters orbiting their host galaxy at a radius of R_g = 8 kpc and of metallicity Z = 0.014 (Sect. 4.2). We discuss the dependence of f_CC on R_g and Z in Sects. 4.3 and 4.4. The probability of finding a star cluster hosting a particular number of Cepheids at the same time is studied in Sect. 4.5, and the Cepheids in the field, but still close to their birth star clusters are dealt with in Sect. 4.6.

4.1. Overview of dynamical evolution of tidally limited star clusters experiencing early gas expulsion

In the present models, stars do not conserve their adiabatic invariants in the gravitational potential during gas expulsion because the gas expulsion timescale τ_M is substantially shorter (τ_M ≈ 0.02 Myr) than the stellar half-mass crossing time, t_cr (see Table 1). As a consequence, the cluster substantially expands while most stars (≈70% of the initial cluster population) escape from the cluster (Tutukov 1978; Lada et al. 1984; Geyer & Burkert 2001; Kroupa et al. 2001; Baumgardt & Kroupa 2007).

On a timescale of 1 Myr to several tens of Myr (depending on the cluster mass), the cluster revirialises, settling at a significantly larger half-mass radius than its initial value (Banerjee & Kroupa 2013, 2017). Further evolution is driven purely by stellar dynamics, where stars escape from the cluster mainly due to slow evaporation with a few energetic ejections, which mostly occur between more massive stars (Baumgardt & Makino 2003; Fujii & Portegies Zwart 2011; Oh et al. 2015; Schoettler et al. 2020). As the cluster mass decreases due to escaping stars and stellar mass loss, the tidal radius decreases, thereby further facilitating the escape of stars until the cluster eventually dissolves completely, leaving behind an unbound stellar stream.

Figure 2 shows the fraction of star clusters that remain gravitationally bound² at a given age t in the tidal field of a galaxy (taking the Galaxy as the specifically computed example). At a particular galactocentric radius, the dissolution time of a cluster is mainly determined by the cluster mass M_ecl. Particularly, at R_g = 8 kpc, all clusters with ≳1600 M_⊙ survive up to the age of 300 Myr (the black lines in Fig. 2), while half of all 400 M_⊙ and 100 M_⊙ clusters dissolve within 250 Myr (medium magenta line) and 90 Myr (medium cyan line), respectively.

Fig. 2. Fraction of star clusters $f_{\mathrm{cl},\mathrm{bound}}^{\mathrm{sim}}$ that are gravitationally bound by time t in the external tidal field of a galaxy as functions of the galactocentric radius Rg, initial cluster mass Mecl and metallicity Z. The galactocentric radius Rg of the orbit is indicated by the line thickness (thick – 4 kpc; medium – 8 kpc; thin – 12 kpc), the cluster mass is indicated by the line colour (cyan – 100 M⊙; magenta – 400 M⊙; black – 1600 M⊙), and the metallicity for the models with Rg = 8 kpc is indicated by the line style (solid – Z = 0.014; dashed – Z = 0.006; dotted – Z = 0.002).

Another important parameter strongly impacting the cluster dissolution timescale is the galactocentric radius of its orbit (Vesperini & Heggie 1997; Baumgardt & Makino 2003). Clusters orbiting closer to the galactic centre dissolve earlier due to the stronger tidal field: at R_g = 4 kpc (thick solid lines in Fig. 2), half of 400 M_⊙ and 100 M_⊙ clusters dissolve within 120 Myr and 50 Myr, respectively. At R_g = 12 kpc (thin solid lines), the majority (≈75%) of 400 M_⊙ clusters survive until at least 300 Myr, while half of 100 M_⊙ clusters dissolve within 150 Myr. All 1600 M_⊙ clusters are preserved for at least 300 Myr even at the smallest radius of 4 kpc. The influence of metallicity on cluster dissolution is negligible, as shown by the dashed and dotted lines for the model with R_g = 8 kpc.

4.2. Dynamical evolution of prospective Cepheids in clusters at R_g = 8 kpc and Z = 0.014

Here, we investigate the significance of different physical mechanisms that release ProCeps from clusters. Stars generally escape their birth clusters by three mechanisms: the rapid change of the gravitational potential during gas expulsion; gradual evaporation due to encounters within the cluster; and fast ejections in close encounters between several bodies.

To identify the mechanism responsible for the escape of the star, we state that any star escaping at speed v_esc > 3σ(0), where σ(0) is the initial 3D velocity dispersion of stars in the embedded cluster (Table 1), is an ejection. It is unlikely for a star to achieve such a high speed without interacting in an unstable system of three or more bodies. The cause of an escape of a star with v_esc < 3σ(0) is identified as the gas expulsion if the star escapes at time t < t_ee, where t_ee = r_t/(0.1σ(0)) and r_t is the tidal radius of the cluster; or, as dynamical evaporation if the star escapes at time t > t_ee. The choice of t_ee is motivated by the typical value of the velocity of stars released due to gas expulsion, which is of the order of σ(0). These stars reach the tidal boundary of the cluster at time ≈r_t/σ(0). The factor of 0.1 is a safety factor to take into account the low velocity tail of the velocity distribution. The same definition of t_ee was adopted by Dinnbier & Kroupa (2020).

The early evolution of embedded star clusters is dominated by gas expulsion, which in our models starts at t_d = 0.6 Myr after formation of the cluster. However, the clusters dynamically evolve prior to t_d. We compare mass segregation of ProCeps and of lower mass stars by measuring the half-mass radius, r_h, for both stellar populations, which is shown in Fig. 3. The simulation of the M_ecl = 100 M_⊙ cluster has r_h for ProCeps significantly smaller than r_h for lower mass stars at ≈t_d, which means that ProCeps are already concentrated near the cluster centre at the time of gas expulsion. The preferential presence of ProCeps near the cluster centre makes them more resilient to gas expulsion, which unbinds stars mainly from the outskirts of the cluster. These are mainly lower mass stars that are located at the outskirts. After gas expulsion, this cluster contains an excess of ProCeps relatively to lower mass stars. More massive clusters (M_ecl ≳ 400 M_⊙) do not have enough time to mass segregate prior t_d, which means that gas expulsion affects ProCeps and lower mass stars to a similar extent, meaning that more massive clusters lose ProCeps at a similar rate as lower mass stars. The mass segregation timescale for ProCeps seen in the simulations is close to the analytic value t_ms obtained from Spitzer (1969), which is listed in Table 1; for M_ecl = 100 M_⊙, the simulated cluster mass segregates before t_d (Table 1 provides t_ms = 0.28 Myr). For M_ecl ≳ 400 M_⊙, the simulated cluster mass segregates around or after t_d (Table 1 provides t_ms ≳ 0.50 Myr)³. Mass segregation continues after revirialisation, but it occurs on a significantly longer timescale (t_ms ≈ 30 Myr) than its initial value t_ms, because gas expulsion decreased the density of the cluster by a factor of 10³.

Fig. 3. Evolution of half-mass radius rh of ProCeps (solid lines) and lower mass stars (dotted lines) for clusters of mass Mecl = 100 M⊙, 400 M⊙, 1600 M⊙ and 6400 M⊙. The clusters orbit the galaxy at Rg = 8 kpc and have metallicity Z = 0.014. The vertical yellow line shows the time of gas expulsion td. The least massive clusters (Mecl = 100 M⊙; cyan) mass segregate their ProCeps before td, which results in ProCeps being more concentrated than lower mass stars. The following gas expulsion then preferentially unbinds stars from outskirts of the cluster, which are lower mass stars, resulting in a cluster that contains disproportionately more ProCeps. In contrast, the more massive clusters do not mass segregate before td, so the stellar loss due to gas expulsion affects ProCeps in a similar way to lower mass stars.

To what extent does the preferential escape of lower mass stars as a consequence of mass segregation increase the relative population of ProCeps and make the lower mass clusters ‘Cepheid rich’? To answer this question, we compare the fraction of ProCeps in clusters (i.e. N_PC, cl/N_PC, tot) to the fraction of lower mass stars in clusters (i.e. N_lm, cl/N_lm, tot) and define the relative abundance of ProCeps as

$$\begin{array}{c}\hfill \eta \equiv \frac{N_{\mathrm{PC},\mathrm{cl}}/N_{\mathrm{PC},\mathrm{tot}}}{N_{\mathrm{lm},\mathrm{cl}}/N_{\mathrm{lm},\mathrm{tot}}}·\end{array}$$(6)

We then plot it as a function of time (Fig. 4). If η is larger than 1 (lower than 1), the cluster has an excess (deficit) in ProCeps relative to the lower mass stars and thus to the overall stellar content because the lower mass stars dominate the overall population by number. At time t_ee, η is close to 1 for clusters with M_ecl ≳ 400 M_⊙, but it is significantly larger (η ≈ 1.8) for the 100 M_⊙ clusters which mass segregated before gas expulsion. After 50 Myr, η decreases in these clusters because of two effects. First, close encounters between ProCeps in the cluster centre are more likely than lower mass stars to eject ProCeps (see Fig. 5 below). Second, mass segregation retains more massive ProCeps in clusters, while the less massive ProCeps are more likely to escape to the field. Since more massive ProCeps evolve to Cepheids earlier than the less massive ones, ProCeps in clusters are preferentially depleted as a result of stellar evolution only. Thus, lower mass clusters contain more ProCeps than would be expected from their total mass.

Fig. 4. Abundance of ProCeps η relatively to the lower mass stars in clusters as a function of time t. The time tee, indicating the transition between stars escaping due to gas expulsion (t < tee) and evaporation (t > tee) is indicated by vertical lines. While the non-mass-segregated clusters (Mecl ≳ 400 M⊙) have only small changes in η, the mass-segregated cluster (Mecl = 100 M⊙) has an increase in η around time tee. To avoid rapid fluctuations due to small number statistics, the plots are terminated when there are fewer than 10 ProCeps in clusters in all realisations of clusters of given mass. The figure represents clusters at Rg = 8 kpc and Z = 0.014.

Clusters in the ≈400 M_⊙ to ≈1600 M_⊙ mass range gradually increase their η with time as the rarefied post-gas-expulsion clusters mass segregate (Fig. 4). The most massive clusters in our study (M_ecl = 6400 M_⊙) have η close to 1 almost independently of their age because of their very long mass segregation timescales. The differences in dynamical evolution of ProCeps between a lower mass (100 M_⊙; movie 1) and more massive (3200 M_⊙; movie 2) star cluster are illustrated in Fig. B.1, and in the accompanying online material⁴.

Figure 5 quantifies the mechanisms responsible for unbinding stars from clusters. For each cluster, we show two bars. The left bar represents ProCeps only, while the right bar represents lower mass stars only. The fraction of lower mass stars that become unbound due to gas expulsion (blue bars) decreases monotonically with increasing cluster mass, which is a well-known fact (e.g. Baumgardt & Kroupa 2007). However, ProCeps do not follow this trend, because they are preferentially retained in lower mass clusters due to dynamical mass segregation. Instead, the fraction of ProCeps that are unbound due to gas expulsion first increases with increasing M_ecl, reaching a maximum of 0.52 for M_ecl = 400 M_⊙ and it further decreases as M_ecl increases. The maximum is caused by mass segregation of ProCeps in clusters with M_ecl ≲ 400 M_⊙, which occurs before gas expulsion. The mechanism that unbinds most of the ProCeps is gas expulsion, with evaporation and ejections playing a secondary role; however, none of these processes are negligible. Dynamical ejections unbind fewer ProCeps than evaporation. In total, 2 to 10% of ProCeps are ejected, where lower mass clusters produce a higher fraction of ejected stars. In contrast, ejections are negligible for lower mass stars.

Fig. 5. Dynamical mechanisms responsible for stellar escapes as a function of the cluster mass. For each cluster mass, we plot two bar plots; the left bar plot represents ProCeps, the right bar plot represents lower mass stars. The colours indicate the relative number of stars in each group which escaped due to gas expulsion (blue), dynamical evaporation (green), dynamical ejection (red), and the rest which either terminated their core Helium burning within the cluster (for ProCeps) or are still present in the cluster by the age of 145 Myr (yellow). We note that in the least massive clusters (Mecl ≲ 200 M⊙), which mass segregate before td, lower mass stars and ProCeps react differently to gas expulsion: as Mecl decreases, more low-mass stars are unbound, while fewer ProCeps are unbound. The figure is for clusters with Rg = 8 kpc and Z = 0.014.

Figure 6 shows the fraction of Cepheids located in their birth clusters p_CC(M_ecl, t) as a function of the cluster age and mass. The values for clusters at R_g = 8 kpc are shown in the middle panel. The histograms are produced by the Bayesian block method to take into account the non-uniformly sampled time when Cepheids occur. The main strengths of the method, which uses variable bin sizes, is its independence of the assumptions about the smoothness and shape of the function and the time resolution of the data (Scargle et al. 2013; Astropy Collaboration 2013). The plots feature all Cepheids formed within the simulations including UpCeps, which are present after time t_He, end. Clusters with M_ecl ≳ 800 M_⊙ have almost time-independent p_CC(M_ecl, t) ≈ 0.5. Less massive clusters have lower p_CC(M_ecl, t), which decreases with time as these clusters dissolve more easily. The higher abundance η of Cepheids in clusters of lower masses, which is due to mass segregation, is partially compensated due to faster dissolution of these clusters; p_CC(M_ecl, t) would be lower in the absence of mass segregation. The fraction p_CC(M_ecl) of cluster Cepheids averaged over time is shown in the left panel of Fig. 7 (red circles).

Fig. 6. Fraction pCC(Mecl, t) of cluster Cepheids as a function of time for clusters orbiting the galaxy at different galactocentric radii: Rg = 4 kpc (upper panel), Rg = 8 kpc (middle panel), and Rg = 12 kpc (lower panel). The cluster mass is indicated by colour. Metallicity is Z = 0.014 for all models. The plot shows all Cepheids including UpCeps, which occur after tHe, end. Generally, pCC(Mecl, t) decreases with time for lower mass clusters (Mecl ≲ 800 M⊙), and it stays approximately constant for clusters more massive than that. The histogram is binned by the Bayesian blocks algorithm. The lines terminate at 300 Myr or at the occurrence of the last Cepheid.

The earlier dissolution of lower mass clusters implies that the Cepheids that occur in these clusters are on average of a larger mass ${\overline{m}}_{\mathrm{Ceph}}$ than those in more massive clusters, as shown in the right panel of Fig. 7, because lower mass Cepheids occur at the time when the majority of lower mass clusters have already dissolved, and the lower mass Cepheids have been released to the field.

Fig. 7. Left panel: fraction of cluster Cepheids pCC(Mecl) (including UpCeps) as a function of cluster mass Mecl integrated over the first 300 Myr of the cluster evolution. Clusters orbiting at different galactocentric radii Rg are denoted by different symbols; clusters with different metallicity Z are denoted by different colour of the symbol. The external field of the galaxy, which impacts the lower mass clusters more, is responsible for the decrease of pCC(Mecl) with Rg. pCC(Mecl) increases with Mecl because more massive clusters withstand gas expulsion and the tidal field of the galaxy better. For lower mass clusters, pCC(Mecl) decreases with metallicity because Cepheids occur later at lower Z, at which time more low-mass clusters have been dissolved. Right panel: mean mass of Cepheids in the cluster as a function of Mecl. At a given metallicity, ${\overline{m}}_{\mathrm{Ceph}}$ peaks at a cluster mass of 200 M⊙ because lower mass clusters do not live long enough to host lower mass Cepheids. The decrease of ${\overline{m}}_{\mathrm{Ceph}}$ with decreasing Z is caused by the decreasing lower mass limit mmin, Ceph for Cepheids with Z.

Conversely, lower mass Cepheids are preferentially hosted by more massive clusters, which survive for longer. It appears that the dependence of ${\overline{m}}_{\mathrm{Ceph}}$ on M_ecl is robust because it is found in clusters at different orbital radii (red squares and triangles in Fig. 7) and metallicities (green and blue circles). Finally, the least massive clusters (M_ecl ≲ 100 M_⊙) have smaller ${\overline{m}}_{\mathrm{Ceph}}$ because their ability to host high-mass ProCeps is limited by the IMF (cf. m_max in Table 1).

The occurrence of Cepheids in clusters therefore depends on several aspects. In particular, we note the dichotomy of p_CC(M_ecl, t) with age for high and low cluster masses, where M_ecl ≈ 800 M_⊙ separates the two regions.

The properties of ProCeps and Cepheids for the lower mass clusters show significant differences between individual realisations for clusters of the same mass. As an example, clusters of mass 200 M_⊙ contain, on average, only 3.5 ProCeps (cf. N_PC and Nmod in Table 1), which implies small number statistics for Cepheids in low mass clusters. Some low mass clusters do not follow the evolutionary path described in this section, where interaction of each individual ProCep with another ProCep or lower mass star decides whether these stars are retained in the cluster or not. For example, all ProCeps can be ejected from the cluster in a single dynamic encounter occurring before gas expulsion (i.e. ≲0.6 Myr). The coarseness is absent in more massive clusters (6400 M_⊙ clusters have ≈125 ProCeps). This should be taken into account when comparing present theoretical findings with observations. Our results describe averages over a large cluster sample, so they are representative, for example, for a cluster population of a galaxy rather than for the individual cluster.

Up to this point, our aim was to theoretically understand the fraction of Cepheids in clusters with little emphasis on quantities that can be directly compared with observations. One quantity of more observational interest is the fraction of Cepheids of given mass, m_Ceph, located in clusters (Fig. 8). In this plot, the Cepheids at any mass bin originate from the whole star cluster population of the ECMF of Eq. (3). The fraction p_CC(m_Ceph) increases with m_Ceph because more massive Cepheids occur earlier, when more lower mass star clusters are still gravitationally bound (it is the same effect which is responsible for the higher mean Cepheid mass ${\overline{m}}_{\mathrm{Ceph}}$ in lower mass clusters as seen in the right panel of Fig. 7).

Fig. 8. Fraction of cluster Cepheids pCC(mCeph) as a function of the Cepheid mass mCeph for Cepheids originating from clusters at different galactocentric radius Rg (indicated by symbols) and of different metallicity Z (indicated by colour). The plot represents Cepheids originating from the whole population of star clusters of the ECMF of Eq. (3). Cepheids of higher mass have a slightly higher probability of being found in clusters.

4.3. The role of galactocentric radius

The clusters at R_g = 4 kpc (upper panel of Fig. 6) and R_g = 12 kpc (lower panel) follow the same trends in p_CC(M_ecl, t) as the cluster at R_g = 8 kpc (middle panel), which was described in the previous section. The value of p_CC(M_ecl, t) decreases with time for lower mass clusters because they dissolve faster than more massive clusters (cf. Fig. 2), which have p_CC(M_ecl, t) almost independent of time. We note that Cepheids occur non-uniformly with time, with ≈80% of all Cepheids occurring before t_He, end ≈ 150 Myr (for Z = 0.014). This influences their statistics, which is relatively good at ≲150 Myr, but relatively poor afterwards, when only a small number of Cepheids can cause visible difference. While most of the clusters follow the expected trend of having larger p_CC(M_ecl, t) at a larger galactocentric radius, there are several exceptions because of low number statistics. An example is that of the M_ecl = 100 M_⊙ clusters (cyan lines), which have lower p_CC(M_ecl, t ≳ 150 Myr) for the R_g = 12 kpc models than for the R_g = 8 kpc models.

The time-averaged fraction of Cepheids in clusters for R_g = 4 kpc (red squares), R_g = 8 kpc (red circles) and R_g = 12 kpc (red triangles) is shown in the left panel of Fig. 7. For lower mass clusters (M_ecl ≲ 400 M_⊙), p_CC(M_ecl) generally increases with increasing R_g as the galactic tidal field becomes less important. Averaged over all the clusters with M_ecl ≲ 400 M_⊙, f_CC = 0.17 for R_g = 4 kpc, f_CC = 0.26 for R_g = 8 kpc, and f_CC = 0.27 for R_g = 12 kpc. On the other hand, more massive clusters (M_ecl ≳ 800 M_⊙) have f_CC ≈ 0.5 almost independently of R_g because they are not significantly affected by the tidal field of the galaxy.

4.4. The role of metallicity

Figure 9 shows the time evolution of p_CC(M_ecl, t) for clusters of different metallicities. At a lower metallicity, Cepheids occur at a lower mass m_min, Ceph and live longer (their t_He, end is larger). Consequently, p_CC(M_ecl, t) reaches lower values near t_He, end for clusters with lower metallicity because the clusters have more time to evaporate or dissolve; this trend is more pronounced in lower mass clusters (p_CC(M_ecl, t_He, end)≲0.1 for Z = 0.002) because they release stars to the field faster. When integrated over time (left panel of Fig. 7), the value of f_CC averaged over clusters with M_ecl ≲ 400 M_⊙ decreases from f_CC = 0.26 for Z = 0.014 via f_CC = 0.24 for Z = 0.006 to f_CC = 0.18 for Z = 0.002. Clusters with M_ecl ≳ 800 M_⊙ have f_CC ≈ 0.5 with no dependence on Z. This shows that the highest sensitivity of p_CC(M_ecl) on metallicity is for lower mass clusters (M_ecl ≲ 800 M_⊙).

Fig. 9. Fraction of cluster Cepheids pCC(Mecl, t) as a function of time for clusters of different metallicities: Z = 0.014 (upper panel), Z = 0.006 (middle panel), Z = 0.002 (lower panel). The cluster mass is indicated by colour. All clusters have Rg = 8 kpc.

This result is the consequence of the longer lifespan of Cepheids at a lower metallicity, at which time more lower mass clusters are dissolved, instead of differences in internal dynamics for clusters of a lower metallicity. This can be deduced from Fig. 2, where the fraction of gravitationally bound clusters at a given time, $f_{\mathrm{cl},\mathrm{bound}}^{\mathrm{sim}}$, is almost independent of metallicity (cf. the solid, dashed and dotted lines). This result is also supported by Hurley et al. (2004), who find a rather minor influence of metallicity on star cluster evolution even for a metallicity substantially lower than considered here (Z = 10⁻⁴). Additionally, the lower mass of Cepheids at a lower metallicity (right panel of Fig. 7) means that their progenitors need a longer timescale for mass segregation and are more likely to escape the cluster during gas expulsion, further reducing the population of cluster Cepheids.

4.5. The number of Cepheids hosted by a cluster at the same time

To interpret the results, we remember that M_ecl is the birth mass of the cluster. The present mass of the cluster hosting Cepheids is typically by a factor of two to three lower than this, depending on its age and initial M_ecl.

Figure 10 shows the time-averaged number of Cepheids, CC_mult, occurring at the same time within the same cluster as a function of the cluster age for clusters of four different masses. Since the number of Cepheids in the same cluster rapidly varies with time around small integer numbers, forming a teeth-like curve, we resort to taking a time average of this quantity. The histogram of the number of Cepheids is discussed below (see also Fig. 11 below). The value of CC_mult increases with time up to t_He, end (Fig. 10; the grey dashed line) because lower mass Cepheids, which occur later are more abundant (according to the adopted model for the IMF), and they spend a longer time within the instability strip; this occurs despite the continuing cluster dissolution with progressing time. More massive clusters, which contain more ProCeps, have higher CC_mult. CC_mult also increases with decreasing metallicity (particularly before t_He, end) because of the lower minimum mass of Cepheids (and thus their higher abundance) in this environment (right panel of Fig. 10). After the expected lifetime of ProCeps t_He, end, CC_mult suddenly drops, but it is non-zero; these Cepheids originate from UpCeps, which are shown by dashed lines (a more detailed study of UpCeps is available in Paper II).

Fig. 10. Time dependence of CCmult for CanCeps (solid lines) and all Cepheids (dashed lines) in clusters of different mass as shown by the colour. The grey vertical dashed lines indicate the time tHe, end for the given metallicity. The histograms are binned by the Bayesian block method. The lines are shown for the interval marked by the occurrence of the first and last Cepheid in clusters in present models (the discretisation to individual Cepheids is the reason why all models do not reach exactly 300 Myr).

How many Cepheids do we expect to see in a particular cluster at the same time? To answer this question, we counted the number of Cepheids in each of our simulated clusters every 30 kyr during the time interval of occurrence of CanCeps; that is, between 15 Myr and t_He, end. A histogram of this count is shown in Fig. 11. We first note that the histogram pertains to observations of a single star cluster, not the whole population of star clusters. For example, the probability of finding at least one Cepheid in a cluster with a birth stellar mass of 6400 M_⊙ is p = 0.317, while for an 400 M_⊙ cluster the probability is only p = 0.017. However, the galaxy forms ≈16× more 400 M_⊙ clusters than 6400 M_⊙ clusters, so the probability of finding at least one Cepheid among 400 M_⊙ clusters is not that different as among 6400 M_⊙ clusters (e.g. p = 0.017 × 16 = 0.272 vs. p = 0.317).

Fig. 11. Normalised histograms of the number of Cepheids located in the same cluster, CCmult, at a given time as a function of cluster mass and metallicity. The numbers above the columns indicate the number fraction of clusters which contained zero to four Cepheids at the same time. The Cepheids taken into account are the CanCeps only. None of the clusters contained more than six Cepheids at the same time.

The probability of finding more Cepheids at the same time decreases with the number of Cepheids we intend to find, and also with decreasing cluster mass (Fig. 11). Clusters harbouring two and more Cepheids at the same time are usually the more massive ones, with mass larger than 400 M_⊙. Clusters with initial mass ≲400 M_⊙ almost never host more than one Cepheid. The highest value of CC_mult found in our simulations is five, and this occurred for the most massive (M_ecl = 6400 M_⊙), and the least metallic (Z = 0.002) cluster. The event of CC_mult = 5 occurred only once in more than 2500 simulations in total, so observations of clusters of masses below ≈6400 M_⊙ with CC_mult > 5 are likely to be rare. This low probability of two or more Cepheids occurring simultaneously in a given cluster agrees very well with the observed statistics based on the MW, M 31, and the majority of Cepheid hosting clusters in the LMC and SMC (Anderson & Riess 2018).

However, some LMC clusters hosting many Cepheids are remarkable exceptions to this general rule, notably NGC 1866 and NGC 2031 that each host ≈20 or more Cepheids (Welch & Stetson 1993; Testa et al. 2007). We assessed whether the occurrence of such clusters is consistent with our simulations in the following way. We assume that a Cepheid in a cluster occurs at a time that is independent of the occurrence of other Cepheids, and that the probability p_ioc of finding a given ProCep as a Cepheid at any time is substantially lower than 1. These assumptions imply Poisson distribution for Cepheids with a mean value of

$$\begin{array}{c}\hfill \lambda =p_{\mathrm{ioc}}{N}_{\mathrm{Ceph}},\end{array}$$(7)

where N_Ceph is the number of Cepheids that occurred during the cluster’s life-time. From the values in the middle panel of Fig. 11 (metallicity similar to NGC 1866), we obtain λ = 0.97 and λ = 0.22 for models with M_ecl = 6400 M_⊙ and M_ecl = 1600 M_⊙, respectively. These values are consistent with the definition of λ, which is linearly proportional to N_Ceph. When scaling λ from the M_ecl = 1600 M_⊙ clusters to the M_ecl = 6400 M_⊙ clusters, we would obtain λ = 0.88 for M_ecl = 6400 M_⊙ clusters. Averaging these two points yields λ ≈ 0.93 for M_ecl = 6400 M_⊙ clusters. Adopting a mass of 1.35 × 10⁵ M_⊙ for NGC 1866 (Fischer et al. 1992), Eq. (7) yields λ = 19.6. Thus, the Poisson distribution yields an expected 20 Cepheids, which is not too far from the known 24 Cepheids in NGC 1866 (Welch & Stetson 1993) whose simultaneous occurrence has a probability of ≈13%. Thus, the wealth of Cepheid members of these two LMC clusters is explained by their high mass, low metallicity (more numerous lower mass stars can become Cepheids), and a cluster age consistent with the occurrence of comparatively old short-period (low-mass) Cepheids. For reference, the mean pulsation periods of Cepheids in NGC 1866 and NGC 2031 reported by Testa et al. (2007) are 2.9 d and 3.2 d, respectively, corresponding to ages of ≈100 Myr assuming period-age relations without rotation (Bono et al. 2005) and approximately 200 Myr for period-age relations assuming average main sequence rotation rates (Anderson et al. 2016).

4.6. Cepheids in the outskirts of star clusters

Figure 12 shows the cumulative distribution of Cepheids according to the distance r from the density centre of their birth cluster. Each star is plotted at the time when it becomes a Cepheid, which means that the stars are shown at different ages. This enables a direct comparison with observations. Approximately 20% of Cepheids whose birth cluster is still gravitationally bound are located at distances from 10 pc to 50 pc. Interestingly, Anderson et al. (2013) find that 5 of 23 (22%) cluster Cepheids in their sample are located at this distance from clusters. Recent searches for Cepheids in MW clusters have primarily reported such new possible ‘halo’ members (Medina et al. 2021; Zhou & Chen 2021). We find approximately 5% of Cepheids located far away (≳1 kpc) from their birth clusters as a result of dynamical ejection. The distribution of Cepheids according to r is independent of the galactocentric radius R_g or metallicity Z.

Fig. 12. Cumulative distribution of Cepheids according to the distance r to the density centre of their birth cluster for the whole population of star clusters from the ECMF of Eq. (3). We only plot the Cepheids whose birth cluster is still gravitationally bound at the time of the occurrence of the Cepheid.

5. The fraction of Cepheids in star clusters

In this section, we provide the likely upper and lower estimate on f_CC. The upper estimate on f_CC is taken from present simulations, where the main mechanisms responsible for cluster dissolution are the early gas expulsion and the tidal field of the galaxy. However, star clusters are subjected to additional dissolution mechanisms, such as encounters with giant molecular clouds (GMCs; Spitzer 1958; Terlevich 1987; Theuns 1992; Jerabkova et al. 2021), which is not included in present simulations. The lower estimate on f_CC is based on the difference between the population of observed open star clusters in the Galaxy and the population of bound star clusters in our models. We note that we do not, a priori, assume the physical mechanism for additional dissolution, we only assume the survival rate of star clusters.

We consider four different models of the galactic cluster population. In the first three models, we assume that at any galactocentric radius R_g, Cepheids are formed in star clusters of the ECMF and initial mass range as given in Sect. 2.3. Since the cluster number density often varies non-trivially as a function of the galactocentric radius with a significant fraction of star formation sometimes located in ring-like structures (e.g. Lewis et al. 2015), we assume rather extreme radial distributions of star clusters over the galaxy to constrain the possible variations in f_CC: (i) all star clusters are located at the same galactocentric radius of R_g = 4 kpc; (ii) all star clusters are located at the same galactocentric radius of R_g = 12 kpc; and (iii) star clusters follow the usual exponential stellar surface density profile of spiral galaxies in the form of Σ_ecl ∝ exp( − R_g/R_d) with the disc scale length of R_d = 2.5 kpc (Binney & Tremaine 2008).

The fourth model (iv) is based on the idea of Pflamm-Altenburg & Kroupa (2008), where the maximum cluster mass M_ecl, max decreases with increasing R_g in disc galaxies as M_ecl, max(R_g) = M_{ecl, max, 0}exp( − γR_g/R_d) because it is set by the local gas surface density. This dependence of the cluster masses on galactocentric radius was documented for the late-type galaxy M 33 by Pflamm-Altenburg et al. (2013), which also demonstrates the distribution of young clusters to be non-stochastic. Following Pflamm-Altenburg & Kroupa (2008), we adopted the slope γ = 1.5, the disc scale length R_d = 4 kpc, and the maximum cluster mass the galaxy can form M_{ecl, max, 0} = 8 × 10⁴ M_⊙ (Weidner et al. 2004). The final value of f_CC is calculated by weighting the ECMF of Eq. (3) (in the interval of 50 M_⊙ to M_ecl, max(R_g)) by the number of Cepheids forming in the particular annulus, which is proportional to 2πR_gexp( − γR_g/R_d).

5.1. The upper estimate on f_CC

For the cluster distributions (i), (ii), (iii), and (iv) and for solar metallicity, we obtain f_CC = 0.30, f_CC = 0.36, f_CC = 0.33, and f_CC = 0.34, respectively (cf. the sixth to ninth columns of Table 2). For the sub-solar metallicities, Z = 0.006 and Z = 0.002, we only have models for R_g = 8 kpc, for which we obtain f_CC = 0.35 and f_CC = 0.32, respectively. Since metallicity has a minor influence on dynamical evolution of our models (cf. Sect. 4.4), we estimate f_CC in sub-solar metallicity clusters by scaling their values from R_g = 8 kpc to R_g = 4 kpc and R_g = 12 kpc by the same factors as is the scaling for the model with Z = 0.014. We obtain f_CC in the range from 0.27 to 0.32 for Z = 0.006, and f_CC in the range from 0.25 to 0.29 for Z = 0.002 (Table 2).

5.2. The lower estimate on f_CC

In the solar vicinity, only ≈10% of embedded star clusters survive to the age of 10 Myr, and ≈4% of embedded star clusters survive to the age of 100 Myr (Lada & Lada 2003), which is mainly attributed to encounters with GMCs. The observed fraction of surviving clusters is smaller than what we obtain in our simulations (cf. Fig. 2), which implies that more Cepheids are released to the field in the real galaxy. However, the inclusion of more realistic (i.e. non-spherically symmetric and evolving with time) GMCs to this type of simulation would be a non-trivial task, which to our knowledge has not been accomplished yet, so we resort to a coarse approximation to obtain an order of magnitude estimate.

Based on Lada & Lada (2003), we approximate the fraction of star clusters which remain gravitationally bound at age t by

$$\begin{array}{c}\hfill f_{\mathrm{cl},\mathrm{bound}}^{\mathrm{obs}}(t)=0.1{(t/10\mathrm{Myr})}^{-2/5}.\end{array}$$(8)

We assume that Eq. (8) is independent of the cluster mass, and of R_g. Then, for each Cepheid in a cluster (which occurs at the cluster age t), we measure the fraction of bound clusters $f_{\mathrm{cl},\mathrm{bound}}^{\mathrm{sim}}(t)$ in our simulations. In the real galaxy, there is $f_{\mathrm{cl},\mathrm{bound}}^{\mathrm{obs}}(t)/f_{\mathrm{cl},\mathrm{bound}}^{\mathrm{sim}}(t)$ times Cepheids in clusters at t, so we take $min(f_{\mathrm{cl},\mathrm{bound}}^{\mathrm{obs}}(t)/f_{\mathrm{cl},\mathrm{bound}}^{\mathrm{sim}}(t),1)$ of each Cepheid as located within a star cluster.

The fraction of cluster Cepheids obtained by this correction for the additional dissolution processes is shown in the second to fifth columns of Table 2. It is smaller by a factor of 10 than without the correction (sixth to ninth columns), which indicates that molecular clouds or a similar disruptive mechanism play likely a crucial role in setting the value of f_CC. These values of f_CC are comparable to what is observed in the galaxy M 31 (f_CC = 0.025; Anderson & Riess 2018).

6. Discussion

6.1. Gas expulsion and mass segregation

The most important parameters of gas expulsion which influence cluster dynamics is the star formation efficiency SFE and the gas expulsion time-scale τ_M (Baumgardt & Kroupa 2007). Although we adopt what appear to be their most likely values (Kroupa et al. 2001; Banerjee & Kroupa 2013; Megeath et al. 2016), we discuss the possible influence on f_CC in case that these parameters are markedly different. For SFE higher than 1/2, gas expulsion unbinds only a small fraction of stars from the cluster. A small fraction of stars is also released when τ_M ≫ t_cr. Nevertheless, lower mass clusters dissolve by 300 Myr due to cluster dynamics only (de la Fuente Marcos & de la Fuente Marcos 2004), while the more massive clusters are likely to retain the majority of their ProCeps. In these cases, we expect f_CC to be somewhat larger than in our models.

On the other hand, if SFE ≲ 1/3 and τ_M ≲ t_cr, the majority of clusters is quickly dispersed forming unbound OB associations with a possible handful of small bound clusters, which are remainders from the most massive and compact clusters (Tutukov 1978). In this case, f_CC would be lower than in our models, possibly similar to our lower estimate on the f_CC.

The present models started as non-mass segregated, but the lower mass clusters (M_ecl ≲ 200 M_⊙) developed mass segregation dynamically before the onset of gas expulsion at t_d = 0.6 Myr. However, the value of t_d is not well constrained. If gas expulsion started earlier (later), mass segregation will develop in clusters less (more) massive than 200 M_⊙.

Mass segregation (either primordial or induced by dynamics) impacts p_CC(M_ecl) in several ways. As shown in Sect. 4.2, mass segregation retains more ProCeps in the cluster in comparison to lower mass stars, which tends to increase f_CC in total. On the other hand, this effect is, to some degree, counteracted by stellar evolution when mass-segregated massive stars explode as supernovae near the cluster centre, causing more pronounced cluster expansion, unbinds more stars and results in faster cluster dissolution (Brinkmann et al. 2017), which would decrease f_CC. The final outcome of the interplay of these two processes, which likely depends on both the cluster mass and Cepheid mass, is difficult to predict without further N-body calculations.

6.2. Initial binary fraction

The assumption of 100% binary fraction under the exemption of initial triple systems may have a significant impact on our results. Paper II in this series studies the dynamics of multiple systems in our simulation in detail. It suffices to state that our simulations suggest that a significant fraction (possibly the majority) of mid-B stars form in triples and quadruple systems instead of binaries. The presence of primordial triples can cause faster mass segregation (because these subsystems are more massive than single stars or binaries) and preserve more ProCeps in clusters increasing f_CC; on the other hand, the dynamical interactions within subsystems and between subsystems and cluster members are likely to cause more ejections, decreasing f_CC. Additional simulations spanning a wide parameter range of initial conditions would be required to study these effects in further detail.

6.3. The external potential of the galaxy

The gravitational potential of the galaxy was modelled using the approximation of Allen & Santillan (1991), even though more recent approximations are available (e.g. Irrgang et al. 2013; Pouliasis et al. 2017). We argue that the details of the potential are not critical for our results because the main influence of the potential is to define the tidal radius, beyond which the stars are removed from the cluster, thus setting the rate of evaporation. The clusters are on circular orbits, so the tidal radius does not change as the cluster revolves the galaxy (apart from cluster mass loss due to evaporation). A cluster of mass M_cl on a circular orbit of radius R_g has a tidal radius given by

$$\begin{array}{c}\hfill r_{\mathrm{t}}=R_{\mathrm{g}}{\left(\frac{M_{\mathrm{cl}}}{4M_{\mathrm{glx}}}\frac{{\gamma }^2}{({\gamma }^2-1)}\right)}^{1/3}\end{array}$$(9)

(Binney & Tremaine 2008, their Eqs. (3.80) and (8.106)), where M_glx is the mass of the galaxy within the cluster orbit and γ = 2Ω/κ is the ratio of the orbital and epicyclic frequency. Since M_glx follows from the magnitude of the rotational curve and γ follows from its slope, which is flat for most spiral galaxies (Bosma 1981; Rubin et al. 1982), implying γ ≈ 1.4, the value of the tidal radius is insensitive of given galaxy model for clusters located in the stellar disc. We note that in other theories of gravity (e.g. MOND), cluster evaporation can occur on a different timescale (faster in the case of MOND) and the tidal structure around the cluster can be asymmetric (Kroupa et al., in prep.).

This study is aimed at general spiral galaxies, which show a large variability in their gravitational potentials, so using a model particularly tuned for the Galaxy would not provide a general result.

6.4. The role of molecular clouds

The lower estimate on f_CC in our models (Table 2) is similar to the observed fraction of cluster Cepheids in M 31, f_CC, M31 = 0.025, (Senchyna et al. 2015; Anderson & Riess 2018). However, significant variations in f_CC can exist among galaxies, and f_CC ≈ 0.07 − 0.11 was estimated in the LMC, SMC, and the Milky Way (Anderson & Riess 2018, cf. their Sect. 3.1).

Here, we provide a tentative update on the Milky Way value of f_CC using the high-quality astrometry from the Gaia mission’s Early Data Release 3 (Gaia Collaboration 2021), updated literature on Galactic open clusters (Loktin & Popova 2017; Cantat-Gaudin et al. 2018, 2020; Ferreira et al. 2021; Hunt & Reffert 2021; Dias et al. 2021), recent investigations of Cepheid membership in Galactic clusters (Breuval et al. 2020; Zhou & Chen 2021; Medina et al. 2021), and the most extensive literature on Galactic Cepheids from the General Catalog of Variable Stars (Samus’ et al. 2017), the Variable Star Index (Watson et al. 2006), and the re-classification of MW Cepheids by Ripepi et al. (2019). A more detailed investigation of this fraction will be presented separately (Cruz Reyes et al., in prep.). We limited our estimate to a volume of 2 kpc to strike a useful balance between statistics and completeness of Cepheids and clusters. Following several quality cuts and detailed inspection of light curves and other available data, we find 217 classical Cepheids within 2 kpc. Of these, 15 are bona fide members of star clusters. We thus find f_CC, MW = 15/217 = 0.069, somewhat lower than the previous value of 11/130 = 0.085 in Anderson & Riess (2018). On one hand, the decrease in f_CC, MW is driven by the greater number of Cepheids that can be placed within 2 kpc thanks to Gaia astrometry. On the other hand, the increase in known clusters resulted in several new cases as well as a reassessment of cluster membership of some other cases. However, it is likely that the list of known clusters within 2 kpc remains significantly less complete than the list of known Cepheids, as evidenced by the recent discovery of more than 600 clusters, mostly at distances beyond 1 kpc (Castro-Ginard et al. 2022). We therefore consider our updated estimate of f_CC, MW = 15/217 = 0.069 a lower limit to the true fraction of cluster Cepheids within 2 kpc from the Sun.

Such large differences between different galaxies cannot be fully explained by lower metallicity (Z = 0.006 and Z = 0.002 was used for the LMC and SMC, respectively) nor by a different strength of the galactic tidal field. As shown in Sects. 4.3 and 4.4, decreasing metallicity from Z = 0.014 to Z = 0.002 decreases f_CC only by a factor of 1.2, and increasing galactocentric radius from R_g = 4 kpc to R_g = 12 kpc increases f_CC also by a factor of 1.2⁵. A possible explanation of the difference of f_CC for different galaxies might stem from different conditions in the interstellar medium in these galaxies, which forms molecular clouds of different masses and/or on different timescales, thus having a different impact on star clusters via the gravitational force. If true, this would imply a bias depending on the gaseous content of the galaxy.

7. Conclusions

We investigated the dynamical evolution of Cepheids in open star clusters with the main aim of estimating the fraction of Cepheids located in clusters and the field. We focused on the influence of the initial cluster mass, orbital radius within its galaxy, and metallicity. The stellar masses follow a realistic initial mass function, and all stars are assumed to form as binaries. At a young age, the clusters expel their non-star-forming gas; this unbinds a significant fraction of stars from the clusters. All the results are obtained directly from simulations.

Lower mass clusters dissolve more easily in the tidal field of the galaxy, releasing a higher fraction of their stars to the field than more massive clusters. For more massive stars, including ProCeps, this process is partially compensated by the rapid mass segregation in very low-mass clusters (M_ecl ≲ 200 M_⊙) occurring before gas expulsion. The fraction of Cepheids located in clusters increases with the cluster mass from ≈0.25 (for M_ecl = 50 M_⊙) to ≈0.60 (for M_ecl = 6400 M_⊙). Cepheids of higher mass are more likely (by ≈30%) to be found in clusters than Cepheids of lower mass because lower mass Cepheids appear at a later time when more clusters dissolve. The evolution of the lower mass (100 M_⊙) and the more massive cluster (3200 M_⊙) is shown as videos in the online material.

Clusters orbiting the galaxy at a larger galactocentric radius R_g take longer to dissolve in the galactic tidal field, retaining more stars (and thus more ProCeps). From R_g = 4 kpc to R_g = 12 kpc, the fraction of Cepheids in clusters f_CC increases by a factor of 1.2. At a lower metallicity, Cepheids occur later, at the time when more star clusters are dissolved. As Z changes from the Solar value of 0.014 to a SMC value of 0.002, f_CC decreases by a factor of 1.2.

At a given time, the majority of star clusters host at most one Cepheid. Generally, the probability of observing two or more Cepheids within the same cluster at the same time increases with the cluster mass and age, and slightly increases with decreasing metallicity (Figs. 10 and 11). For example, in an initially 6400 M_⊙ cluster at Z = 0.002, the probability to observe two, three, and four Cepheids simultaneously is 7.4%, 3.5%, and 2.5%, respectively. The highest number of Cepheids at the same time detected in our simulations in any cluster is five, indicating that observing more Cepheids in one cluster than this is rare. However, we find that the extremely Cepheid-rich LMC clusters NGC 1866 and NGC 2031 are consistent with our results once their higher mass is accounted for. Our simulations allow us to understand the existence of such a large number of Cepheids in these clusters as being due to the combination of very high cluster mass, lower LMC metallicity, and a cluster age near the age when the least massive (most numerous) Cepheid progenitors become Cepheids.

Assuming that star clusters are formed with the mass function of Eq. (3) and then evolve in a Milky Way-like galaxy, we obtain f_CC = 0.30 to 0.36 for the whole galaxy with Z = 0.014. This is an upper limit on f_CC because it neglects other cluster dissolving mechanisms than gas expulsion, cluster internal dynamics, and the galactic tidal field. Given that the observed difference of f_CC between the galaxy M 31 (f_CC ≈ 0.025) on the one hand, and the Galaxy, LMC, and SMC (f_CC ≈ 0.09) on the other hand, differs by a substantially larger factor than our simulations permit with a varying tidal field strength and metallicity, it is unlikely that the observed difference can be explained by variations of the galactic tidal field strength and metallicity only.

Using an analytical model to incorporate additional cluster disruption mechanisms to match the observed number of dissolved clusters (Lada & Lada 2003), we obtain substantially smaller values of f_CC = 0.015 to 0.026. This disruption is mainly attributed to interactions with giant molecular clouds, which is not included in the present N-body simulations. Thus, f_CC strongly depends on cluster disruption mechanisms, while its dependence on the galactic tidal field or metallicity is of secondary importance.

Movies

Movie 1 associated with Fig. B.1 (movie1_100_Msun) Access here

Movie 2 associated with Fig. B.1 (movie2_3200_Msun) Access here

Acknowledgments

We thank an anonymous referee for the useful comments and suggestions, which improved the quality of the paper. We would like to thank Sverre Aarseth for continuously developing the NBODY family of numerical integrators and Sambaran Banerjee for his advice about specifying integration parameters for the simulations. We appreciate the support of the ESO IT team, which was vital for performing the presented simulations. We thank Mauricio Cruz Reyes for assistance in obtaining an updated value for f_CC based on Gaia EDR3 data. This research made use of Astropy (http://www.astropy.org), a community developed core Python package for Astronomy (Astropy Collaboration 2013) and Matplotlib Python Package (Hunter 2007). F.D. and P.K. acknowledge the Scientific Visitors Programme of European Southern Observatory in Garching, which made this project possible. F.D. and P.K. acknowledge support through grant 20-21855S from the Czech Grant Agency. R.I.A. acknowledges funding provided by SNSF Eccellenza Professorial Fellowship PCEFP2_194638 and funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement No. 947660).

References

Appendix A: Detection of Cepheids originating from stellar mergers

Since we detect Cepheids automatically according to their evolution in the HR diagram (the description is in Sect. 3), it is important to make sure that the recipe is able to deal not only with stars evolving in isolation, but also with more complicated cases including binary star evolution and stellar mergers. Typical examples (for Z = 0.014) are shown in Figure A.1.

Fig. A.1. Examples of HR diagrams where the evolutionary path of the primary had been altered by a merger with the secondary. These plots show examples from the library of simulations calculated by NBODY6 for Z = 0.014. The evolution from the terminal age main-sequence to core He burning is represented by black points, during core He burning by orange points, and afterwards by green points. The coalescence causes a sudden increase of luminosity, which can be seen as the discontinuity in the evolutionary tracks. The instability strip and ZAMS are shown by the thick and thin black lines, respectively. Left panel: Coalescence with relatively massive secondary (7 M⊙), which increases the primary mass so much that the merger object does not become a Cepheid. Middle panel: Coalescence with a lower mass secondary, which does not prevent the star from becoming a Cepheid. Right panel: Lower mass primary (4.3 M⊙), which would not have evolved to a Cepheid in isolation, becomes a Cepheid after a coalescence.

The left panel details the evolution of a 11 M_⊙ and 7 M_⊙ binary with initial orbital period of 27 days, which coalesces when the primary increases in radius as it evolves out of the MS (the discontinuity before core He burning). The merger mass is 14.1 M_⊙. Core He burning (orange points) then commences leftwards of the blue boundary of the instability strip, and the star monotonically decreases in effective temperature during this stage, which is dissimilar to the loop experienced by Cepheids. For this reason, the star does not fulfil condition (iii) for the detection of Cepheids (Sect. 3).

The middle panel shows a 8.4 M_⊙ primary with a secondary of 1.1 M_⊙ and initial orbital period of 13 days. In this case, the coalescence only slightly increases the mass, and the ProCep becomes a Cepheid.

The right panel shows coalescence of a primary of ZAMS mass 4.3 M_⊙ with its 1.6 M_⊙ companion. For a metallicity of Z = 0.014, the initial mass of the primary is lower than the mass threshold for Cepheids, which is m_min, Ceph = 4.7 M_⊙. After the coalescence, the merged star has a mass of 5.3 M_⊙, and it becomes a Cepheid.

The cases in the middle and right panels are detected as Cepheids by the automatic routine. The case in the left panel is not detected as a Cepheid. This demonstrates that the automatic procedure is able to identify Cepheids also for more complicated evolutionary tracks including mergers.

Appendix B: Comparison of the evolution of ProCeps in a lower mass (100 M_⊙) and a more massive (3200 M_⊙) cluster

The lower mass cluster, which contains three ProCeps in total, segregates its more massive stars before the onset of gas expulsion at 0.6 Myr (upper row of Figure B.1). At this time (upper middle panel), two ProCeps are already near the cluster centre. The third ProCep, which is located at x ≈ −0.6 pc, returns to the cluster centre, where it interacts with the other two Cepheids, forming a densely bound system. This results in all three ProCeps being concentrated at the cluster centre, while lower mass stars are more dispersed or escape the cluster (upper right panel). The complicated interaction between the most massive bodies ejects one of the ProCep and a lower mass star at 7.2 Myr (shown only in the online material).

Fig. B.1. Evolutionary snapshots for an example of a lower mass (100 M⊙; upper row) and more massive star clusters (3200 M⊙; lower row). The plots are shown at the beginning of the simulations (left column), immediately before the beginning of gas expulsion (middle column), and after gas expulsion; the time is indicated above the panels. ProCeps are marked red, while stars with mass below mmin, Ceph or above mmax, Ceph are marked blue. Compact objects (seen only in the online material) are marked black. Stellar mass is proportional to the area of the dots. The magenta circle shows the cluster tidal radius. We note that the length scale in the right column is inflated by a factor of 10. The evolution of both star clusters up to 12 Myr can be seen as movies in the online material.

In contrast, the more massive star cluster (lower row of Figure B.1) does not mass segregate ProCeps before gas expulsion, which results in the population of ProCeps being dispersed throughout the cluster to a similar degree as the lower mass stars. Gas expulsion then unbinds a comparable fraction of ProCeps as low mass stars from the cluster.

All Tables

All Figures

Fig. 1. HR diagrams outputted from NBODY6 for single stars of ZAMS mass of 4 M⊙ (left panel), 7 M⊙ (middle panel), and 10 M⊙ (right panel). The evolutionary paths are plotted after the stars left the main sequence. The colour indicates the metallicity, which is Z = 0.002 (blue), Z = 0.006 (green), and Z = 0.014 (red). Large crosses denote the time of core helium burning. For the purpose of this plot, we selected single stars, which evolve well in isolation from the other stars in the cluster. Thick solid lines show the Cepheid instability strip for the metallicities considered (Anderson et al. 2016). Thin solid lines represent ZAMS for stars of corresponding metallicity.

In the text

Fig. 2. Fraction of star clusters $f_{\mathrm{cl},\mathrm{bound}}^{\mathrm{sim}}$ that are gravitationally bound by time t in the external tidal field of a galaxy as functions of the galactocentric radius Rg, initial cluster mass Mecl and metallicity Z. The galactocentric radius Rg of the orbit is indicated by the line thickness (thick – 4 kpc; medium – 8 kpc; thin – 12 kpc), the cluster mass is indicated by the line colour (cyan – 100 M⊙; magenta – 400 M⊙; black – 1600 M⊙), and the metallicity for the models with Rg = 8 kpc is indicated by the line style (solid – Z = 0.014; dashed – Z = 0.006; dotted – Z = 0.002).

In the text

Fig. 3. Evolution of half-mass radius rh of ProCeps (solid lines) and lower mass stars (dotted lines) for clusters of mass Mecl = 100 M⊙, 400 M⊙, 1600 M⊙ and 6400 M⊙. The clusters orbit the galaxy at Rg = 8 kpc and have metallicity Z = 0.014. The vertical yellow line shows the time of gas expulsion td. The least massive clusters (Mecl = 100 M⊙; cyan) mass segregate their ProCeps before td, which results in ProCeps being more concentrated than lower mass stars. The following gas expulsion then preferentially unbinds stars from outskirts of the cluster, which are lower mass stars, resulting in a cluster that contains disproportionately more ProCeps. In contrast, the more massive clusters do not mass segregate before td, so the stellar loss due to gas expulsion affects ProCeps in a similar way to lower mass stars.

In the text

Fig. 4. Abundance of ProCeps η relatively to the lower mass stars in clusters as a function of time t. The time tee, indicating the transition between stars escaping due to gas expulsion (t < tee) and evaporation (t > tee) is indicated by vertical lines. While the non-mass-segregated clusters (Mecl ≳ 400 M⊙) have only small changes in η, the mass-segregated cluster (Mecl = 100 M⊙) has an increase in η around time tee. To avoid rapid fluctuations due to small number statistics, the plots are terminated when there are fewer than 10 ProCeps in clusters in all realisations of clusters of given mass. The figure represents clusters at Rg = 8 kpc and Z = 0.014.

In the text

Fig. 5. Dynamical mechanisms responsible for stellar escapes as a function of the cluster mass. For each cluster mass, we plot two bar plots; the left bar plot represents ProCeps, the right bar plot represents lower mass stars. The colours indicate the relative number of stars in each group which escaped due to gas expulsion (blue), dynamical evaporation (green), dynamical ejection (red), and the rest which either terminated their core Helium burning within the cluster (for ProCeps) or are still present in the cluster by the age of 145 Myr (yellow). We note that in the least massive clusters (Mecl ≲ 200 M⊙), which mass segregate before td, lower mass stars and ProCeps react differently to gas expulsion: as Mecl decreases, more low-mass stars are unbound, while fewer ProCeps are unbound. The figure is for clusters with Rg = 8 kpc and Z = 0.014.

In the text

Fig. 6. Fraction pCC(Mecl, t) of cluster Cepheids as a function of time for clusters orbiting the galaxy at different galactocentric radii: Rg = 4 kpc (upper panel), Rg = 8 kpc (middle panel), and Rg = 12 kpc (lower panel). The cluster mass is indicated by colour. Metallicity is Z = 0.014 for all models. The plot shows all Cepheids including UpCeps, which occur after tHe, end. Generally, pCC(Mecl, t) decreases with time for lower mass clusters (Mecl ≲ 800 M⊙), and it stays approximately constant for clusters more massive than that. The histogram is binned by the Bayesian blocks algorithm. The lines terminate at 300 Myr or at the occurrence of the last Cepheid.

In the text

Fig. 7. Left panel: fraction of cluster Cepheids pCC(Mecl) (including UpCeps) as a function of cluster mass Mecl integrated over the first 300 Myr of the cluster evolution. Clusters orbiting at different galactocentric radii Rg are denoted by different symbols; clusters with different metallicity Z are denoted by different colour of the symbol. The external field of the galaxy, which impacts the lower mass clusters more, is responsible for the decrease of pCC(Mecl) with Rg. pCC(Mecl) increases with Mecl because more massive clusters withstand gas expulsion and the tidal field of the galaxy better. For lower mass clusters, pCC(Mecl) decreases with metallicity because Cepheids occur later at lower Z, at which time more low-mass clusters have been dissolved. Right panel: mean mass of Cepheids in the cluster as a function of Mecl. At a given metallicity, ${\overline{m}}_{\mathrm{Ceph}}$ peaks at a cluster mass of 200 M⊙ because lower mass clusters do not live long enough to host lower mass Cepheids. The decrease of ${\overline{m}}_{\mathrm{Ceph}}$ with decreasing Z is caused by the decreasing lower mass limit mmin, Ceph for Cepheids with Z.

In the text

Fig. 8. Fraction of cluster Cepheids pCC(mCeph) as a function of the Cepheid mass mCeph for Cepheids originating from clusters at different galactocentric radius Rg (indicated by symbols) and of different metallicity Z (indicated by colour). The plot represents Cepheids originating from the whole population of star clusters of the ECMF of Eq. (3). Cepheids of higher mass have a slightly higher probability of being found in clusters.

In the text

	Fig. 9. Fraction of cluster Cepheids pCC(Mecl, t) as a function of time for clusters of different metallicities: Z = 0.014 (upper panel), Z = 0.006 (middle panel), Z = 0.002 (lower panel). The cluster mass is indicated by colour. All clusters have Rg = 8 kpc.
In the text

Fig. 10. Time dependence of CCmult for CanCeps (solid lines) and all Cepheids (dashed lines) in clusters of different mass as shown by the colour. The grey vertical dashed lines indicate the time tHe, end for the given metallicity. The histograms are binned by the Bayesian block method. The lines are shown for the interval marked by the occurrence of the first and last Cepheid in clusters in present models (the discretisation to individual Cepheids is the reason why all models do not reach exactly 300 Myr).

In the text

Fig. 11. Normalised histograms of the number of Cepheids located in the same cluster, CCmult, at a given time as a function of cluster mass and metallicity. The numbers above the columns indicate the number fraction of clusters which contained zero to four Cepheids at the same time. The Cepheids taken into account are the CanCeps only. None of the clusters contained more than six Cepheids at the same time.

In the text

	Fig. 12. Cumulative distribution of Cepheids according to the distance r to the density centre of their birth cluster for the whole population of star clusters from the ECMF of Eq. (3). We only plot the Cepheids whose birth cluster is still gravitationally bound at the time of the occurrence of the Cepheid.
In the text

Fig. A.1. Examples of HR diagrams where the evolutionary path of the primary had been altered by a merger with the secondary. These plots show examples from the library of simulations calculated by NBODY6 for Z = 0.014. The evolution from the terminal age main-sequence to core He burning is represented by black points, during core He burning by orange points, and afterwards by green points. The coalescence causes a sudden increase of luminosity, which can be seen as the discontinuity in the evolutionary tracks. The instability strip and ZAMS are shown by the thick and thin black lines, respectively. Left panel: Coalescence with relatively massive secondary (7 M⊙), which increases the primary mass so much that the merger object does not become a Cepheid. Middle panel: Coalescence with a lower mass secondary, which does not prevent the star from becoming a Cepheid. Right panel: Lower mass primary (4.3 M⊙), which would not have evolved to a Cepheid in isolation, becomes a Cepheid after a coalescence.

In the text

Fig. B.1. Evolutionary snapshots for an example of a lower mass (100 M⊙; upper row) and more massive star clusters (3200 M⊙; lower row). The plots are shown at the beginning of the simulations (left column), immediately before the beginning of gas expulsion (middle column), and after gas expulsion; the time is indicated above the panels. ProCeps are marked red, while stars with mass below mmin, Ceph or above mmax, Ceph are marked blue. Compact objects (seen only in the online material) are marked black. Stellar mass is proportional to the area of the dots. The magenta circle shows the cluster tidal radius. We note that the length scale in the right column is inflated by a factor of 10. The evolution of both star clusters up to 12 Myr can be seen as movies in the online material.

In the text