11institutetext:  Astronomical Institute, Faculty of Mathematics and Physics,
Charles University in Prague, V Holešovičkách 2, 180 00 Praha 8, Czech
Republic
11email<EMAIL_ADDRESS>
11email<EMAIL_ADDRESS>22institutetext: Institute of Physics,
Laboratory of Astrophysics, École Polytechnique Fédérale de Lausanne (EPFL),
Observatoire de Sauverny, 1290 Versoix, Switzerland
22email<EMAIL_ADDRESS>33institutetext: Helmholtz-Institut für
Strahlen- und Kernphysik, University of Bonn, Nussallee 14-16, 53115 Bonn,
Germany

# On the dynamical evolution of Cepheids in star clusters

František Dinnbier 11 Richard I. Anderson 22 Pavel Kroupa 1133

(Received ; accepted ??)

We investigate the occurrence of classical (type-I) Cepheid variable stars
(henceforth: Cepheids) in dynamically evolving star clusters from birth to an
age of approximately 300 Myr. The clusters are modelled by the Aarseth code
nbody6, and they feature a realistic stellar initial mass function and initial
binary star population, single star and binary star evolution, expulsion of
the primordial gas, and the tidal field of the galaxy. Our simulations provide
the first detailed dynamical picture of how frequently Cepheids remain
gravitationally bound to their birth clusters versus how frequently they occur
in the field. They allow us to quantify the relevance of various cluster
escape mechanisms and how they depend on stellar mass. Overall, the
simulations agree with the empirical picture that a small fraction ($\approx
10\%$) of Cepheids reside in clusters, that cluster halo membership is
relatively common, and that the majority of Cepheid hosting clusters have only
a single Cepheid member. Additionally, the simulations predict that a) Cepheid
progenitors are much more likely to escape from low-mass than higher-mass
clusters; b) higher-mass (long-period) Cepheids are $\approx 30\%$ more likely
to be found in clusters than low-mass (short-period) Cepheids; c) the
clustered Cepheid fraction increases with galactocentric radius since cluster
dispersal is less efficient at greater radii; d) a lower metallicity reduces
the overall clustered Cepheid fraction because the lower minimum mass of
Cepheids leaves more time for cluster dispersal (this primarily affects short-
period Cepheids); e) high-mass clusters are much more likely to have more than
one Cepheid member at any given time, in particular at a lower metallicity. We
interpret the results as outcomes of various aspects of star cluster dynamics.
The comparison of predicted and observed clustered Cepheid fractions, $f_{\rm
CC}$, highlights the need for additional cluster disruption mechanisms, most
likely encounters with giant molecular clouds, to explain the observed
fractions.

###### Key Words.:

Galaxies: stellar content, open clusters and associations: general, Stars:
variables: Cepheids, Stars: kinematics and dynamics

##  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 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 et al. 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. 2018a; 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 location. 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,\rm{true}}\lesssim 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); d) $f_{CC,\rm{true}}$ is generally higher for the younger
long-period Cepheids (Anderson & Riess 2018). The most significant exceptions
to this rule are very massive ($\approx 5\times 10^{4}\,\,\mathrm{M}_{\odot}$)
clusters in the Large Magellanic Cloud, such as NGC 1866 or NGC 2031, that
each contains 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 of
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,\rm{eff}}$, compared to the true number of clustered
extragalactic Cepheids, and any uncertainties in $f_{CC,\rm{eff}}$ affect
stellar association bias corrections for the distance ladder linearly
(Anderson & Riess 2018). Understanding the difference between
$f_{CC,\rm{eff}}$ and $f_{CC,\rm{true}}$ and possible variations in
$f_{\rm{CC}}$ among different galaxies will be very 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 1995b; 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
and 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 Milky Way, LMC, SMC, and M31
(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 have 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 1995a; 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 forward modeling 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, Anderson,
and Kroupa, in prep.).

This _article_ is structured as follows. Sect. 2 describes the set up of the
simulations, initial conditions, and other choices made to create the cluster
library. Sect. 3 presents how Cepheids are identified and traced in the NBODY6
simulations. Sect. 4 investigates the evolution and occurrence of Cepheids in
dynamically evolving star clusters as a function of the cluster mass,
metallicity and galactocentric radius. Sect. 5 seeks to estimate the number of
Cepheids occurring in clusters, $f_{CC,\mathrm{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 Sect. 7 summarizes 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 4th 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 time-scales, 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 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 gets 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 adopt the gas expulsion
model by Kroupa et al. (2001), where the gravitational potential $\phi_{\rm
gas}(r,t)$ of the gas is represented by a Plummer sphere that starts to
exponentially decrease after a delay time $t_{\rm d}$, i.e.

$\phi_{\rm gas}(r,t)=\cases{-}\frac{GM_{\rm gas}}{\sqrt{a_{\rm
gas}^{2}+r^{2}}}\;\;t\leq t_{\rm d},\\\ -\frac{GM_{\rm gas}}{\sqrt{a_{\rm
gas}^{2}+r^{2}}}\exp{\\{-(t-t_{\rm d})/\tau_{\rm M}\\}}\;\;t>t_{\rm d}.{}$ (1)

In the gas model, $G$, $M_{\rm gas}$, $a_{\rm gas}$, and $\tau_{\rm 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 time-scale, respectively. For the parameters of the gas expulsion
model, we adopt the values which 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 adopt $a_{\rm
gas}=a_{\rm ecl}$, where $a_{\rm ecl}$ is the Plummer length parameter of
stars, $M_{\rm gas}=2M_{\rm ecl}$, where $M_{\rm ecl}$ is the initial stellar
mass of the embedded cluster (this implies the star formation efficiency
$\mathrm{SFE}\equiv M_{\rm ecl}/(M_{\rm ecl}+M_{\rm gas})=1/3$), $\tau_{\rm
M}=a_{\rm gas}/c_{\rm II}\approx 0.03\,\mathrm{Myr}$, where $c_{\rm
II}=10\,\mathrm{km}\,\,\mathrm{s}^{-1}$ is the sound speed in ionised
hydrogen, and $t_{\rm d}=0.6\,\mathrm{Myr}$ is approximately the lifetime of
the ultra-compact HII region and the prior deeply embedded phase.

The simulated clusters orbit a Milky Way-like galaxy on circular trajectories.
The escaping stars are integrated in the gravitational field of the galaxy, so
they form tidal tails, and stellar evolution is followed also in the tidal
tails. The gravitational potential of the galaxy is represented by the model
of Allen & Santillan (1991), which consists of three components: the bulge
modelled as a Plummer potential of mass $1.4\times
10^{10}\,\mathrm{M}_{\odot}$ and scale-length $0.38\,\mathrm{kpc}$; the
stellar disc modelled as a Miyamoto-Nagai (Miyamoto & Nagai 1975) potential of
mass $8.5\times 10^{10}\,\mathrm{M}_{\odot}$ and radial and vertical scale-
lengths $5.3\,\mathrm{kpc}$ and $0.25\,\mathrm{kpc}$, respectively; and a dark
matter halo of a two-component potential which is nearly logarithmic at large
galactocentric distances, and which encompasses a mass of $8\times
10^{11}\,\mathrm{M}_{\odot}$ within a sphere of $100\,\mathrm{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 are generated by the software package MCLUSTER (Küpper et al.
2011) as Plummer models of stellar mass $M_{\rm ecl}$ and Plummer parameter
$a_{\rm ecl}$. Employing the relation between the embedded cluster radius and
its mass in stars (Marks & Kroupa 2012),

$\frac{a_{\rm ecl}}{\rm{pc}}=0.077\left(\frac{M_{\rm
ecl}}{\,\mathrm{M}_{\odot}}\right)^{0.13},$ (2)

we reduce the two parameter space $(M_{\rm ecl},a_{\rm ecl})$ to one parameter
family of models characterised by cluster mass $M_{\rm ecl}$. The clusters are
initially in virial equilibrium with the primordial gas. The stars are
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 are
initially not mass segregated.

The initial masses of stars are sampled from the model of the stellar initial
mass function (IMF) due to Kroupa (2001), where the maximum allowed stellar
mass, $m_{\rm max}$, within the cluster is a function of $M_{\rm ecl}$
according to the $m_{\rm max}-M_{\rm ecl}$ relation (Weidner et al. 2010).
Thus, only massive clusters harbour the most massive stars. The minimum
stellar mass is $0.08\,\mathrm{M}_{\odot}$. All stars form in binaries with
orbital parameters for stars of $<3\,\mathrm{M}_{\odot}$ according to the
initial distribution function for periods, eccentricities and mass ratios of
Kroupa (1995a), and more massive stars ($>3\,\mathrm{M}_{\odot}$) according to
the binary properties of O stars (Sana et al. 2012; Moe & Di Stefano 2017).
The method how to generate a realistic initial binary population without
affecting the IMF is described by Oh et al. (2015) and Oh & Kroupa (2016).
Since the evolution of the binaries is at the focus of paper II, the details
of the initial binary population and the evolution processes are provided
there.

Model | $\rm{Nmod}$ | $M_{\rm ecl}$ | $r_{\rm h}$ | $m_{\rm max}$ | $N_{\rm PC}$ | $\sigma(0)$ | $t_{\rm cr}$ | $t_{\rm rlx}$ | $t_{\rm ms}$ |
---|---|---|---|---|---|---|---|---|---|---
|  | [$\,\mathrm{M}_{\odot}$] | [$\,\mathrm{pc}$] | [$\,\mathrm{M}_{\odot}$] |  | [$\,\mathrm{km}\,\,\mathrm{s}^{-1}$] | [$\,\mathrm{Myr}$] | [$\,\mathrm{Myr}$] | [$\,\mathrm{Myr}$] |
M1 | 2 | 6400 | 0.299 | 92 | 252 | 10.3 | 0.028 | 58 | 2.10 |
M2 | 4 | 3200 | 0.274 | 63 | 277 | 7.6 | 0.035 | 40 | 1.44 |
M3 | 8 | 1600 | 0.251 | 42 | 225 | 5.6 | 0.044 | 28 | 1.00 |
M4 | 16 | 800 | 0.231 | 26 | 233 | 4.1 | 0.054 | 19 | 0.70 |
M5 | 32 | 400 | 0.211 | 17 | 214 | 3.1 | 0.068 | 14 | 0.50 |
M6 | 64 | 200 | 0.194 | 10 | 226 | 2.3 | 0.084 | 10 | 0.37 |
M7 | 128 | 100 | 0.178 | 6.6 | 217 | 1.7 | 0.104 | 7.7 | 0.28 |
M8 | 256 | 50 | 0.163 | 4.1 | 151 | 1.2 | 0.129 | 6.1 | 0.22 |
Table 1: List of star cluster models. Model name and number $\rm{Nmod}$ of
realisations of the model (cluster of the same global properties but with
initial conditions from a different random seed) are listed in the first and
second columns. The following columns list the cluster stellar mass $M_{\rm
ecl}$, initial half-mass radius $r_{\rm h}$, mass of the most massive star in
the cluster $m_{\rm max}$ (mean value from all realisations), the number of
ProCeps $N_{\rm PC}$ (total for all realisations), the initial 3D velocity
dispersion $\sigma(0)$, the half-mass radius crossing time $t_{\rm cr}$, the
half-mass radius relaxation time $t_{\rm rlx}$ (Binney & Tremaine 2008, their
eqs. 7.106 and 7.108), and the mass segregation time-scale $t_{\rm ms}$
(Spitzer 1969, his eq. 9) for $10\,\mathrm{M}_{\odot}$ stars. To calculate the
velocity dispersion, $\sigma(0)$, and time-scales $t_{\rm cr}$ and $t_{\rm
rlx}$ for embedded clusters, we multiply their values for gas-free clusters by
$(1/\mathrm{SFE})^{1/2}$, $(\mathrm{SFE})^{1/2}$ and $(1/\mathrm{SFE})^{3/2}$,
respectively.

###  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$, i.e.

$\xi_{\rm ecl}(M_{\rm ecl})\equiv\frac{\mathrm{d}N_{\rm ecl}(M_{\rm
ecl})}{\mathrm{d}M_{\rm ecl}}\propto M_{\rm ecl}^{-2},$ (3)

where $\mathrm{d}N_{\rm ecl}(M_{\rm ecl})$ is the number of embedded star
clusters forming in a mass bin of size $\mathrm{d}M_{\rm ecl}$.

Our models span cluster mass range of $M_{\rm
ecl}\in(50\,\mathrm{M}_{\odot},6400\,\mathrm{M}_{\odot})$, where the lower
mass limit follows from the $m_{\rm max}-M_{\rm ecl}$ relation, which requires
$m_{\rm max}\gtrsim m_{\rm 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 ($\approx 8.5\times 10^{3}\,\mathrm{M}_{\odot}$; Johnson et al. 2017). The
dependence $\xi_{\rm ecl}(M_{\rm ecl})\propto M_{\rm 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_{\rm ecl}$.

###  2.4 The library of simulations

The cluster masses, $M_{\rm ecl}$, span the range of $50\,\mathrm{M}_{\odot}$
to $6400\,\mathrm{M}_{\odot}$, and they are sampled by multiples of $2$, i.e.
$M_{\rm
ecl}=50\,\mathrm{M}_{\odot},100\,\mathrm{M}_{\odot},..,6400\,\mathrm{M}_{\odot}$.
The models are listed in Table 1. We run each model Nmod times with different
random seed for stellar masses, positions and velocities, where Nmod$=2^{i}$,
and $i$ is the index after the model name, e.g. models M4 have $\rm{Nmod}=16$.
We choose these values of Nmod to have the same total mass of star clusters in
each logarithmic cluster mass bin $\mathrm{d}M_{\rm ecl}$ when used with the
ECMF of eq. (3). This facilitates calculating any quantity $\langle q\rangle$
for the whole population of clusters (and thus its observed value integrated
over a galaxy) as $\langle
q\rangle=(\sum_{i=1}^{8}q_{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$ 111Although
the total mass of clusters in each logarithmic mass bin is the same, the
number of stars which evolve to Cepheids is different in each bin because of
the $m_{\rm max}-M_{\rm ecl}$ relation. For example, the least massive
clusters do not populate the stellar IMF up to the maximum mass $m_{\rm
max,Ceph}$ of stars which evolve to Cepheids because they have $m_{\rm
max}<m_{\rm max,Ceph}$. ).

We calculate the population of clusters described in the previous paragraph
for three galactocentric radii: $R_{\rm g}=4\,\mathrm{kpc}$, $R_{\rm
g}=8\,\mathrm{kpc}$, and $R_{\rm g}=12\,\mathrm{kpc}$ for solar metallicity,
i.e. $Z=0.014$. We also calculate the same population of clusters for
galactocentric radius $R_{\rm g}=8\,\mathrm{kpc}$ and subsolar metallicity of
$Z=0.006$ and $Z=0.002$. Thus, for each considered choice of $R_{\rm g}$ and
$Z$, we calculate 510 simulations in total.

##  3 Identification of Cepheids

Isolated stars with zero age main sequence (ZAMS) mass in the appropriate mass
range $(m_{\rm min,Ceph},m_{\rm 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_{\rm
min,Ceph},m_{\rm
max,Ceph})=(3.2\,\mathrm{M}_{\odot},11.0\,\mathrm{M}_{\odot})$ for $Z=0.002$,
$(m_{\rm min,Ceph},m_{\rm
max,Ceph})=(4.0\,\mathrm{M}_{\odot},12.0\,\mathrm{M}_{\odot})$ for $Z=0.006$,
and $(m_{\rm min,Ceph},m_{\rm
max,Ceph})=(4.7\,\mathrm{M}_{\odot},13.0\,\mathrm{M}_{\odot})$ for $Z=0.014$.
For the purpose of this work, we term all stars with initial mass in the range
of $(m_{\rm min,Ceph},m_{\rm max,Ceph})$ prospective Cepheids (ProCephs), and
stars with initial mass below $m_{\rm 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 which originate from a coalescence in a binary with both
components having a ZAMS mass below $m_{\rm min,Ceph}$; we call these Cepheids
as upscattered Cepheids (UpCeps). For the adopted evolutionary tracks, stars
with mass $m_{\rm min,Ceph}$ terminate core Helium burning at time $t_{\rm
He,end}$ of $300\,\mathrm{Myr}$, $190\,\mathrm{Myr}$ and $146\,\mathrm{Myr}$
for $Z=0.002$, $Z=0.006$, and $Z=0.014$, respectively.

Figure 1: HR diagrams outputted from nbody6 for single stars of ZAMS mass of
$4\,\mathrm{M}_{\odot}$ (left panel), $7\,\mathrm{M}_{\odot}$ (middle panel),
and $10\,\mathrm{M}_{\odot}$ (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.

We search Cepheids automatically in the HR diagrams produced by nbody6. For
each star in the simulation, we record its temperature and luminosity whenever
its stellar evolution is updated, which, during core helium burning, typically
occurs every $0.01$ to $0.3\,\mathrm{Myr}$ depending on the Cepheid mass and
exact position in the HR diagram. In this way, we obtain densely sampled HR
diagrams even during the rapid evolution on the instability strip as shown in
Fig. 1. Then, the algorithm diagnoses each star and classifies it as a Cepheid
if the star fulfils all of the following criteria: (i) The star’s initial mass
is in the mass range $(m_{\rm min,Ceph}/2,m_{\rm max,Ceph})$, the lower value
of the mass range is to include a possible merger of two $0.5m_{\rm 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_{\rm 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, i.e. the
star lies on the instability strip during its core Helium burning at least on
one output file. We use 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.

To obtain the properties of the Cepheid (e.g. its distance from the cluster
centre or binarity), we first determine the times when the Cepheid crosses the
instability strip and then evaluate 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 ignore the first passage through the instability strip because
of its very short duration), we determine the time randomly from one of the
two intervals on the instability strip weighted by their duration. The
dynamics of Cepheids usually changes little when being 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\,\mathrm{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 $\approx 2\,\mathrm{pc}$) as a result of the cluster half-
mass radii (several parsec) and mass segregation of ProCeps. The distance of
$10\,\mathrm{pc}$ is comparable to the tidal radius of the clusters, so if a
star exceeds this distance, it is likely to become unbound and recedes. 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_{\rm ecl}$ at given time $t$ is denoted $p_{\rm CC}(M_{\rm
ecl},t)$. The fraction of Cepheids in clusters of a particular cluster mass
$M_{\rm ecl}$ throughout the time interval $(0,300\,\mathrm{Myr})$ is denoted
$p_{\rm CC}(M_{\rm ecl})$,

$p_{\rm CC}(M_{\rm ecl})=\frac{\sum_{i}\;p_{\rm CC}(M_{\rm ecl},t_{i})\Delta
t_{i}}{\sum_{i}\;\Delta t_{i}},$ (4)

where $\Delta 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, i.e. formed in all star clusters within the galaxy over their life-
time is denoted $f_{\rm CC}$,

$f_{\rm CC}=\frac{\sum_{i}\;\xi_{\rm ecl}(M_{{\rm ecl},i})p_{\rm CC}(M_{{\rm
ecl},i})\Delta M_{{\rm ecl},i}}{\sum_{i}\;\xi_{\rm ecl}(M_{{\rm ecl},i})\Delta
M_{{\rm ecl},i}},$ (5)

where $\xi_{\rm ecl}(M_{\rm ecl})$ is the ECMF of eq. (3), and the summation
goes over all cluster mass intervals. We stress that throughout this work,
$M_{\rm 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_{\rm CC}(m_{\rm
Ceph})$ which are of a given stellar mass $m_{\rm 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, which is
populated by the Kroupa IMF from the mass interval of
$0.08\,\mathrm{M}_{\odot}$ to $50\,\mathrm{M}_{\odot}$. 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 relative 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 which unbind ProCeps from clusters and
therefore determine the value of $f_{\rm CC}$ as a function of cluster age and
mass for clusters orbiting their host galaxy at a radius of $R_{\rm
g}=8\,\mathrm{kpc}$ and of metallicity $Z=0.014$ (Sect. 4.2). We discuss the
dependence of $f_{\rm CC}$ on $R_{\rm g}$ and $Z$ in Sections 4.3 and 4.4. The
probability to find 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

Figure 2: The fraction of star clusters $f_{\rm cl,bound}^{\rm sim}$ which are
gravitationally bound by time $t$ in the external tidal field of a galaxy as a
function of the galactocentric radius $R_{\rm g}$, initial cluster mass
$M_{\rm ecl}$ and metallicity $Z$. The galactocentric radius $R_{\rm g}$ of
the orbit is indicated by the line thickness (thick - $4\,\mathrm{kpc}$;
medium - $8\,\mathrm{kpc}$; thin - $12\,\mathrm{kpc}$), the cluster mass is
indicated by the line colour (cyan - $100\,\mathrm{M}_{\odot}$; magenta -
$400\,\mathrm{M}_{\odot}$; black - $1600\,\mathrm{M}_{\odot}$), and the
metallicity for the models with $R_{\rm g}=8\,\mathrm{kpc}$ is indicated by
the line style (solid - $Z=0.014$; dashed - $Z=0.006$; dotted - $Z=0.002$).

In the present models, stars do not conserve their adiabatic invariants in the
gravitational potential during gas expulsion because the gas expulsion time
scale $\tau_{\rm M}$ is substantially shorter ($\tau_{\rm M}\approx
0.02\,\mathrm{Myr}$) than the stellar half-mass crossing time, $t_{\rm cr}$
(see Table 1). As a consequence, the cluster substantially expands while most
stars ($\approx 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 time-scale of $1\,\mathrm{Myr}$ to several tens of Myr (depending on the
cluster mass), the cluster revirialises, settling at a significantly larger
half-mass radius than was 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 occur mostly 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
bound222In this work, we consider a star cluster to be any group of at least
ten c.m. bodies (a c.m. body is either a single star or a regularised binary
or hierarchical system), which is located within the cluster tidal radius.  at
a given age $t$ in the tidal field of a galaxy (taking here 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_{\rm
ecl}$. Particularly, at $R_{\rm g}=8\,\mathrm{kpc}$, all clusters with
$\gtrsim 1600\,\mathrm{M}_{\odot}$ survive up to the age of
$300\,\mathrm{Myr}$ (the black lines in Fig. 2), while half of all
$400\,\mathrm{M}_{\odot}$ and $100\,\mathrm{M}_{\odot}$ clusters dissolve
within $250\,\mathrm{Myr}$ (medium magenta line) and $90\,\mathrm{Myr}$
(medium cyan line), respectively.

Another important parameter strongly impacting the cluster dissolution time-
scale 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_{\rm
g}=4\,\mathrm{kpc}$ (thick solid lines in Fig. 2), half of
$400\,\mathrm{M}_{\odot}$ and $100\,\mathrm{M}_{\odot}$ clusters dissolve
within $120\,\mathrm{Myr}$ and $50\,\mathrm{Myr}$, respectively. At $R_{\rm
g}=12\,\mathrm{kpc}$ (thin solid lines), the majority ($\approx$ 75 %) of
$400\,\mathrm{M}_{\odot}$ clusters survive until at least $300\,\mathrm{Myr}$,
while half of $100\,\mathrm{M}_{\odot}$ clusters dissolve within
$150\,\mathrm{Myr}$. All $1600\,\mathrm{M}_{\odot}$ clusters are preserved at
least for $300\,\mathrm{Myr}$ even at the smallest radius of
$4\,\mathrm{kpc}$. The influence of metallicity on cluster dissolution is
negligible as shown by the dashed and dotted lines for the model with $R_{\rm
g}=8\,\mathrm{kpc}$.

Figure 3: Evolution of the half-mass radius $r_{\rm h}$ of ProCeps (solid
lines) and lower mass stars (dotted lines) for clusters of mass $M_{\rm
ecl}=100\,\mathrm{M}_{\odot}$, $400\,\mathrm{M}_{\odot}$,
$1600\,\mathrm{M}_{\odot}$ and $6400\,\mathrm{M}_{\odot}$. The clusters orbit
the galaxy at $R_{\rm g}=8\,\mathrm{kpc}$ and have metallicity $Z=0.014$. The
vertical yellow line shows the time of gas expulsion $t_{\rm d}$. The least
massive clusters ($M_{\rm ecl}=100\,\mathrm{M}_{\odot}$; cyan) mass segregate
their ProCeps before $t_{\rm d}$, 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 which contains disproportionately more
ProCeps. In contrast, the more massive clusters do not mass segregate before
$t_{\rm d}$, so the stellar loss due to gas expulsion affects ProCeps in a
similar way as lower mass stars.

###  4.2 Dynamical evolution of prospective Cepheids in clusters at $R_{\rm
g}=8\,\mathrm{kpc}$ and $Z=0.014$

Here, we investigate the significance of different physical mechanisms which
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 define
that any star escaping at speed $v_{\rm esc}>3\sigma(0)$, where $\sigma(0)$ is
the initial 3D velocity dispersion of stars in the embedded cluster (Table 1),
is an ejection. It is unlikely that a star achieves 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_{\rm esc}<3\sigma(0)$ is identified as the gas
expulsion if the star escapes at time $t<t_{\rm ee}$, where $t_{\rm ee}=r_{\rm
t}/(0.1\sigma(0))$ and $r_{\rm t}$ is the tidal radius of the cluster; or as
dynamical evaporation if the star escapes at time $t>t_{\rm ee}$. The choice
of $t_{\rm ee}$ is motivated by the typical value of the velocity of stars
released due to gas expulsion, which is of the order of $\sigma(0)$. These
stars reach the tidal boundary of the cluster at time $\approx r_{\rm
t}/\sigma(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_{\rm
ee}$ was adopted by Dinnbier & Kroupa (2020).

Figure 4:  The abundance of ProCeps $\eta$ relative to the lower mass stars in
clusters as a function of time $t$. The time $t_{\rm ee}$, indicating the
transition between stars escaping due to gas expulsion ($t<t_{\rm ee}$) and
evaporation ($t>t_{\rm ee}$) is indicated by vertical lines. While the non-
mass segregated clusters ($M_{\rm ecl}\gtrsim 400\,\mathrm{M}_{\odot}$) have
only small changes in $\eta$, the mass segregated cluster ($M_{\rm
ecl}=100\,\mathrm{M}_{\odot}$) has an increase in $\eta$ around time $t_{\rm
ee}$. To avoid rapid fluctuations due to small number statistics, the plots
are terminated when there is less than $10$ ProCeps in clusters in all
realisations of clusters of given mass. The figure represents clusters at
$R_{\rm g}=8\,\mathrm{kpc}$ and $Z=0.014$.  Figure 5:  The 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\,\mathrm{Myr}$
(yellow). Note that in the least massive clusters ($M_{\rm ecl}\lesssim
200\,\mathrm{M}_{\odot}$), which mass segregate before $t_{\rm d}$, lower mass
stars and ProCeps react differently on gas expulsion: as $M_{\rm ecl}$
decreases, more low mass stars get unbound, while fewer ProCeps get unbound.
The figure is for clusters with $R_{\rm g}=8\,\mathrm{kpc}$ and $Z=0.014$.

The early evolution of embedded star clusters is dominated by gas expulsion,
which in our models starts at $t_{\rm d}=0.6\,\mathrm{Myr}$ after formation of
the cluster. However, the clusters dynamically evolve prior to $t_{\rm d}$. We
compare mass segregation of ProCeps and of lower mass stars by measuring the
half-mass radius, $r_{\rm h}$, for both stellar populations, which is shown in
Fig. 3. The simulation of the $M_{\rm ecl}=100\,\mathrm{M}_{\odot}$ cluster
has $r_{\rm h}$ for ProCeps significantly smaller than $r_{\rm h}$ for lower
mass stars at $\approx t_{\rm d}$, which means that ProCeps are concentrated
near the cluster centre already 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_{\rm ecl}\gtrsim 400\,\mathrm{M}_{\odot}$) do
not have enough time to mass segregate prior $t_{\rm 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 time-scale for ProCeps seen in the simulations is close
to the analytic value $t_{\rm ms}$ obtained from Spitzer (1969), which is
listed in Table 1; for $M_{\rm ecl}=100\,\mathrm{M}_{\odot}$, the simulated
cluster mass segregates before $t_{\rm d}$ (Table 1 provides $t_{\rm
ms}=0.28\,\mathrm{Myr}$). For $M_{\rm ecl}\gtrsim 400\,\mathrm{M}_{\odot}$,
the simulated cluster mass segregates around or after $t_{\rm d}$ (Table 1
provides $t_{\rm ms}\gtrsim 0.50\,\mathrm{Myr}$) 333The formula used for
calculating $t_{\rm ms}$ according to Spitzer (1969) was derived under
idealised assumptions, so it might vary by a factor of several. The close
agreement between our simulations is probably coincidental for our particular
cluster profile, IMF and the choice of stars with the mass
$>10\,\mathrm{M}_{\odot}$ to be mass segregated. . Mass segregation continues
after revirialisation, but it occurs on a significantly longer time-scale
($t_{\rm ms}\approx 30\,\mathrm{Myr}$) than its initial value $t_{\rm ms}$,
because gas expulsion decreased the density of the cluster by a factor of
$10^{3}$.

To which 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_{\rm PC,cl}/N_{\rm
PC,tot}$) to the fraction of lower mass stars in clusters (i.e. $N_{\rm
lm,cl}/N_{\rm lm,tot}$), and define the relative abundance of ProCeps as

$\eta\equiv\frac{N_{\rm PC,cl}/N_{\rm PC,tot}}{N_{\rm lm,cl}/N_{\rm lm,tot}},$
(6)

and plot it as a function of time (Fig. 4). If $\eta$ 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_{\rm ee}$,
$\eta$ is close to $1$ for clusters with $M_{\rm ecl}\gtrsim
400\,\mathrm{M}_{\odot}$, but it is significantly larger ($\eta\approx 1.8$)
for the $100\,\mathrm{M}_{\odot}$ clusters which mass segregated before gas
expulsion. After $50\,\mathrm{Myr}$, $\eta$ decreases in these clusters
because of two effects. First, close encounters between ProCeps in the cluster
centre eject ProCeps more likely than lower mass stars (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 what would be expected
from their total mass.

Clusters in the mass range $\approx 400\,\mathrm{M}_{\odot}$ to $\approx
1600\,\mathrm{M}_{\odot}$ gradually increase their $\eta$ with time as the
rarefied post gas expulsion clusters mass segregate (Fig. 4). The most massive
clusters in our study ($M_{\rm ecl}=6400\,\mathrm{M}_{\odot}$) have $\eta$
close to $1$ almost independent of their age because of their very long mass
segregation time-scale. The differences in dynamical evolution of ProCeps
between a lower mass ($100\,\mathrm{M}_{\odot}$; movie 1) and more massive
($3200\,\mathrm{M}_{\odot}$; movie 2) star cluster are illustrated in Fig. 14
of Appendix B, and in the accompanying online material
444https://zenodo.org/record/5786512 .

Figure 6:  The fraction $p_{\rm CC}(M_{\rm ecl},t)$ of cluster Cepheids as a
function of time for clusters orbiting the galaxy at different galactocentric
radii: $R_{\rm g}=4\,\mathrm{kpc}$ (upper panel), $R_{\rm g}=8\,\mathrm{kpc}$
(middle panel), and $R_{\rm g}=12\,\mathrm{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 $t_{\rm He,end}$.
Generally, $p_{\rm CC}(M_{\rm ecl},t)$ decreases with time for lower mass
clusters ($M_{\rm ecl}\lesssim 800\,\mathrm{M}_{\odot}$), 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\,\mathrm{Myr}$ or at the occurrence of the last Cepheid.  Figure 7: Left
panel: The fraction of cluster Cepheids $p_{\rm CC}(M_{\rm ecl})$ (including
UpCeps) as a function of cluster mass $M_{\rm ecl}$ integrated over the first
$300\,\mathrm{Myr}$ of the cluster evolution. Clusters orbiting at different
galactocentric radii $R_{\rm g}$ 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 more the lower mass clusters,
is responsible for the decrease of $p_{\rm CC}(M_{\rm ecl})$ with $R_{\rm g}$.
$p_{\rm CC}(M_{\rm ecl})$ increases with $M_{\rm ecl}$ because more massive
clusters withstand better gas expulsion and the tidal field of the galaxy.
$p_{\rm CC}(M_{\rm ecl})$ for lower mass clusters decreases with metallicity
because Cepheids occur later at lower $Z$, at which time more low mass
clusters has been dissolved. Right panel: The mean mass of Cepheids in the
cluster as a function of $M_{\rm ecl}$. At a given metallicity,
$\overline{m}_{\rm Ceph}$ peaks at a cluster mass of $200\,\mathrm{M}_{\odot}$
because lower mass clusters do not live long enough to host lower mass
Cepheids. The decrease of $\overline{m}_{\rm Ceph}$ with decreasing $Z$ is
caused by the decreasing lower mass limit $m_{\rm min,Ceph}$ for Cepheids with
$Z$.

Fig. 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 which get 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 which get unbound due to gas
expulsion first increases with increasing $M_{\rm ecl}$, reaching a maximum of
$0.52$ for $M_{\rm ecl}=400\,\mathrm{M}_{\odot}$ and further decreases as
$M_{\rm ecl}$ increases. The maximum is caused by mass segregation of ProCeps
in clusters with $M_{\rm ecl}\lesssim 400\,\mathrm{M}_{\odot}$, which occurs
before gas expulsion. The mechanism which unbinds most of the ProCeps is gas
expulsion, with evaporation and ejections playing a secondary role; however
none of these processes is negligible. Dynamical ejections unbind fewer
ProCeps than evaporation. In total, $2$ to $10$ % of ProCeps get ejected,
where lower mass clusters produce a higher fraction of ejected stars. In
contrast, ejections are negligible for lower mass stars.

Figure 8:  The fraction of cluster Cepheids $p_{\rm CC}(m_{\rm Ceph})$ as a
function of the Cepheid mass $m_{\rm Ceph}$ for Cepheids originating from
clusters at different Galactocentric radius $R_{\rm g}$ (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 to be
found in clusters.

Fig. 6 shows the fraction of Cepheids located in their birth clusters $p_{\rm
CC}(M_{\rm ecl},t)$ as a function of the cluster age and mass. The values for
clusters at $R_{\rm g}=8\,\mathrm{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 et al. 2013). The plots
feature all Cepheids formed within the simulations including UpCeps, which are
present after time $t_{\rm He,end}$. Clusters with $M_{\rm ecl}\gtrsim
800\,\mathrm{M}_{\odot}$ have almost time independent $p_{\rm CC}(M_{\rm
ecl},t)\approx 0.5$. Less massive clusters have lower $p_{\rm CC}(M_{\rm
ecl},t)$, which decreases with time as these clusters dissolve more easily.
The higher abundance $\eta$ of Cepheids in clusters of lower masses, which is
due to mass segregation, is partially compensated due to faster dissolution of
these clusters; $p_{\rm CC}(M_{\rm ecl},t)$ would be lower in the absence of
mass segregation. The fraction $p_{\rm CC}(M_{\rm ecl})$ of cluster Cepheids
averaged over time is shown in the left panel of Fig. 7 (red circles).

Figure 9:  The fraction of cluster Cepheids $p_{\rm CC}(M_{\rm ecl},t)$ as a
function of time for clusters of different metallicity: $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 $R_{\rm g}=8\,\mathrm{kpc}$.

The earlier dissolution of lower mass clusters implies that the Cepheids which
occur in these clusters are on average of a larger mass $\overline{m}_{\rm
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. Conversely, lower mass Cepheids are preferentially
hosted by more massive clusters, which survive for longer. It appears that the
dependence of $\overline{m}_{\rm Ceph}$ on $M_{\rm 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_{\rm ecl}\lesssim 100\,\mathrm{M}_{\odot}$) have smaller
$\overline{m}_{\rm Ceph}$ because their ability to host high-mass ProCeps is
limited by the IMF (cf. $m_{\rm max}$ in Table 1).

The occurrence of Cepheids in clusters therefore depends on several aspects.
In particular, we note the dichotomy of $p_{\rm CC}(M_{\rm ecl},t)$ with age
for high and low cluster masses, where $M_{\rm ecl}\approx
800\,\mathrm{M}_{\odot}$ 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\,\mathrm{M}_{\odot}$ contain
on average only $3.5$ ProCeps (c.f. $N_{\rm PC}$ and $\rm{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. $\lesssim
0.6\,\mathrm{Myr}$). The coarseness is absent in more massive clusters
($6400\,\mathrm{M}_{\odot}$ clusters have $\approx 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.

So far, we aimed at theoretically understanding of the fraction of Cepheids in
clusters with little emphasis on quantities which can be directly compared
with observations. One quantity of more observational interest is the fraction
of Cepheids of given mass, $m_{\rm Ceph}$, which is 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_{\rm CC}(m_{\rm
Ceph})$ increases with $m_{\rm 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}_{\rm Ceph}$ in lower mass clusters as seen in the right panel of
Fig. 7).

###  4.3 The role of galactocentric radius

The clusters at $R_{\rm g}=4\,\mathrm{kpc}$ (upper panel of Fig. 6) and
$R_{\rm g}=12\,\mathrm{kpc}$ (lower panel) follow the same trends in $p_{\rm
CC}(M_{\rm ecl},t)$ as the cluster at $R_{\rm g}=8\,\mathrm{kpc}$ (middle
panel), which was described in the previous section: The value of $p_{\rm
CC}(M_{\rm ecl},t)$ decreases with time for lower mass clusters because they
dissolve faster than more massive clusters (c.f. Fig. 2), which have $p_{\rm
CC}(M_{\rm ecl},t)$ almost independent of time. We note that Cepheids occur
non-uniformly with time, with $\approx 80$% of all Cepheids occurring before
$t_{\rm He,end}\approx 150\,\mathrm{Myr}$ (for $Z=0.014$). This influences
their statistics, which is relatively good at $\lesssim 150\,\mathrm{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_{\rm CC}(M_{\rm ecl},t)$ at a larger galactocentric radius,
there are several exceptions because of low number statistics. An example are
the $M_{\rm ecl}=100\,\mathrm{M}_{\odot}$ clusters (cyan lines), which have
lower $p_{\rm CC}(M_{\rm ecl},t\gtrsim 150\,\mathrm{Myr})$ for the $R_{\rm
g}=12\,\mathrm{kpc}$ models than for the $R_{\rm g}=8\,\mathrm{kpc}$ models.

The time-averaged fraction of Cepheids in clusters for $R_{\rm
g}=4\,\mathrm{kpc}$ (red squares), $R_{\rm g}=8\,\mathrm{kpc}$ (red circles)
and $R_{\rm g}=12\,\mathrm{kpc}$ (red triangles) is shown in the left panel of
Fig. 7. For lower mass clusters ($M_{\rm ecl}\lesssim
400\,\mathrm{M}_{\odot}$), $p_{\rm CC}(M_{\rm ecl})$ generally increases with
increasing $R_{\rm g}$ as the galactic tidal field becomes less important.
Averaged over all the clusters with $M_{\rm ecl}\lesssim
400\,\mathrm{M}_{\odot}$, $f_{\rm CC}=0.17$ for $R_{\rm g}=4\,\mathrm{kpc}$,
$f_{\rm CC}=0.26$ for $R_{\rm g}=8\,\mathrm{kpc}$, and $f_{\rm CC}=0.27$ for
$R_{\rm g}=12\,\mathrm{kpc}$. On the other hand, more massive clusters
($M_{\rm ecl}\gtrsim 800\,\mathrm{M}_{\odot}$), have $f_{\rm CC}\approx 0.5$
almost independently of $R_{\rm g}$ because they are not much affected by the
tidal field of the galaxy.

###  4.4 The role of metallicity

Fig. 9 shows the time evolution of $p_{\rm CC}(M_{\rm ecl},t)$ for clusters of
different metallicity. At a lower metallicity, Cepheids occur at a lower mass
$m_{\rm min,Ceph}$ and live longer (their $t_{\rm He,end}$ is larger).
Consequently, $p_{\rm CC}(M_{\rm ecl},t)$ reaches lower values near $t_{\rm
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_{\rm CC}(M_{\rm ecl},t_{\rm He,end})\lesssim 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_{\rm CC}$ averaged over clusters with
$M_{\rm ecl}\lesssim 400\,\mathrm{M}_{\odot}$ decreases from $f_{\rm CC}=0.26$
for $Z=0.014$ via $f_{\rm CC}=0.24$ for $Z=0.006$ to $f_{\rm CC}=0.18$ for
$Z=0.002$. Clusters with $M_{\rm ecl}\gtrsim 800\,\mathrm{M}_{\odot}$ have
$f_{\rm CC}\approx 0.5$ with no dependence on $Z$. This shows that the highest
sensitivity of $p_{\rm CC}(M_{\rm ecl})$ on metallicity is for lower mass
clusters ($M_{\rm ecl}\lesssim 800\,\mathrm{M}_{\odot}$).

This result is the consequence of the longer life-span of Cepheids at a lower
metallicity, at which time more lower mass clusters is 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_{\rm cl,bound}^{\rm sim}$, is almost independent of
metallicity (c.f. 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^{-4}$). Additionally, the lower mass of
Cepheids at a lower metallicity (right panel of Fig. 7) means that their
progenitors need a longer time-scale for mass segregation, and are more likely
to escape the cluster during gas expulsion, further reducing the population of
cluster Cepheids.

Figure 10:  The time dependence of $CC_{\rm mult}$ 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 $t_{\rm He,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\,\mathrm{Myr}$).  Figure 11:  Normalised histograms of the number of
Cepheids located in the same cluster, $CC_{\rm mult}$, 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.  Figure 12:
The 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 plot only the Cepheids whose birth
cluster is still gravitationally bound at the time of the occurrence of the
Cepheid.

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

To interpret the results, we remember that $M_{\rm 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_{\rm ecl}$.

Fig. 10 shows the time-averaged number of Cepheids, $CC_{\rm 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_{\rm mult}$ increases with time up to $t_{\rm
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_{\rm mult}$. $CC_{\rm
mult}$ also increases with decreasing metallicity (particularly before $t_{\rm
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
life-time of ProCeps $t_{\rm He,end}$, $CC_{\rm 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).

How many Cepheids do we expect to see in a particular cluster at the same
time? To answer this question, we count the number of Cepheids in each of our
simulated clusters every $30\,\mathrm{kyr}$ during the time interval of
occurrence of CanCeps, i.e. between $15\,\mathrm{Myr}$ and $t_{\rm He,end}$.
An histogram of this count is shown in Fig. 11. We note first 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\,\mathrm{M}_{\odot}$ is $p=0.317$, while for an
$400\,\mathrm{M}_{\odot}$ cluster the probability is only $p=0.017$. However,
the galaxy forms $\approx 16\times$ more $400\,\mathrm{M}_{\odot}$ clusters
than $6400\,\mathrm{M}_{\odot}$ clusters, so the probability of finding at
least one Cepheid among $400\,\mathrm{M}_{\odot}$ clusters is not that
different as among $6400\,\mathrm{M}_{\odot}$ clusters (e.g. $p=0.017\times
16=0.272$ vs. $p=0.317$).

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\,\mathrm{M}_{\odot}$. Clusters with initial mass $\lesssim
400\,\mathrm{M}_{\odot}$ almost never host more than one Cepheid. The highest
value of $CC_{\rm mult}$ found in our simulations is five, and this occurred
for the most massive ($M_{\rm ecl}=6400\,\mathrm{M}_{\odot}$), and the least
metallic ($Z=0.002$) cluster. The event of $CC_{\rm mult}=5$ occurred only
once in more than $2500$ simulations in total, so observations of clusters of
mass below $\approx 6400\,\mathrm{M}_{\odot}$ with $CC_{\rm mult}>5$ are
likely to be rare. This low probability for two or more Cepheids to occur
simultaneously in a given cluster agrees very well with the observed
statistics based on the MW, M31, 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 $\approx 20$
or more Cepheids (Welch & Stetson 1993; Testa et al. 2007). We assess 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 which is
independent of the occurrence of other Cepheids, and that the probability
$p_{\rm ioc}$ to find a given ProCep as a Cepheid at any time is substantially
smaller than 1. These assumptions imply Poisson distribution for Cepheids with
mean value

$\lambda=p_{\rm ioc}N_{\rm Ceph},$ (7)

where $N_{\rm Ceph}$ is the number of Cepheids which occurred during the
cluster’s life-time. From the values in the middle panel of Fig. 11
(metallicity similar to NGC 1866), we obtain $\lambda=0.97$ and $\lambda=0.22$
for models with $M_{\rm ecl}=6400\,\mathrm{M}_{\odot}$ and $M_{\rm
ecl}=1600\,\mathrm{M}_{\odot}$, respectively. These values are consistent with
the definition of $\lambda$, which is linearly proportional to $N_{\rm Ceph}$.
When scaling $\lambda$ from the $M_{\rm ecl}=1600\,\mathrm{M}_{\odot}$
clusters to the $M_{\rm ecl}=6400\,\mathrm{M}_{\odot}$ clusters, we would
obtain $\lambda=0.88$ for $M_{\rm ecl}=6400\,\mathrm{M}_{\odot}$ clusters.
Averaging these two points yields $\lambda\approx 0.93$ for $M_{\rm
ecl}=6400\,\mathrm{M}_{\odot}$ clusters. Adopting a mass of $1.35\times
10^{5}\,\mathrm{M}_{\odot}$ for NGC 1866 (Fischer et al. 1992), eq. (7) yields
$\lambda=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 $\approx 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 2103 reported by Testa et al. (2007)
are $2.9\,$d and $3.2$ d, respectively, corresponding to ages of $\approx 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

Fig. 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 the distance from $10\,\mathrm{pc}$ to $50\,\mathrm{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 ($\gtrsim 1\,\mathrm{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_{\rm g}$ or
metallicity $Z$.

##  5 The fraction of Cepheids in star clusters

In this Section, we provide the likely upper and lower estimate on $f_{\rm
CC}$. The upper estimate on $f_{\rm 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, e.g. 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_{\rm 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. 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_{\rm 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_{\rm CC}$:
(i) all star clusters are located at the same galactocentric radius of $R_{\rm
g}=4\,\mathrm{kpc}$; (ii) all star clusters are located at the same
galactocentric radius of $R_{\rm g}=12\,\mathrm{kpc}$; (iii) star clusters
follow the usual exponential stellar surface density profile of spiral
galaxies in the form of $\Sigma_{\rm ecl}\propto\exp{(-R_{\rm g}/R_{\rm d})}$
with the disc scale length $R_{\rm d}=2.5\,\mathrm{kpc}$ (Binney & Tremaine
2008).

The fourth model (iv) is based on the idea of Pflamm-Altenburg & Kroupa
(2008), where the maximum cluster mass $M_{\rm ecl,max}$ decreases with
increasing $R_{\rm g}$ in disc galaxies as $M_{\rm ecl,max}(R_{\rm g})=M_{\rm
ecl,max,0}\exp{(-\gamma R_{\rm g}/R_{\rm d})}$ because it is set by the local
gas surface density. This dependence of the cluster masses on galactocentric
radius has been documented for the late-type galaxy M33 by Pflamm-Altenburg et
al. (2013), which also demonstrates the distribution of young clusters to be
non-stochastic. Following Pflamm-Altenburg & Kroupa (2008), we adopt the slope
$\gamma=1.5$, the disc scale length $R_{\rm d}=4\,\mathrm{kpc}$, and the
maximum cluster mass which the galaxy can form $M_{\rm
ecl,max,0}=8.10^{4}\,\mathrm{M}_{\odot}$ (Weidner et al. 2004). The final
value of $f_{\rm CC}$ is calculated by weighting the ECMF of eq. (3) (in the
interval of $50\,\mathrm{M}_{\odot}$ to $M_{\rm ecl,max}(R_{\rm g})$) by the
number of Cepheids forming in the particular annulus, which is proportional to
$2\pi R_{\rm g}\exp{(-\gamma R_{\rm g}/R_{\rm d})}$.

$Z$ | lower $f_{\rm CC}$ | upper $f_{\rm CC}$
---|---|---
| (i) | (ii) | (iii) | (iv) | (i) | (ii) | (iii) | (iv)
0.014 | 0.026 | 0.015 | 0.024 | 0.021 | 0.30 | 0.36 | 0.33 | 0.34
0.006 | 0.024 | 0.014 | 0.023 | 0.019 | 0.27 | 0.32 | 0.29 | 0.30
0.002 | 0.024 | 0.014 | 0.022 | 0.017 | 0.25 | 0.29 | 0.27 | 0.27
Table 2: The total fraction $f_{\rm CC}$ of Cepheids in star clusters for the
whole galaxy derived under different assumptions about the cluster
distribution within the galaxy (indicated by Roman numerals; they are defined
in Sect. 5), cluster dissolution processes and metallicity $Z$. Columns 2-5
assume pronounced cluster dissolution with the fraction $f_{\rm cl,bound}^{\rm
obs}(t)$ of clusters remaining gravitationally bound by the age $t$ given by
eq. (8). This is likely the lower estimate on $f_{\rm CC}$. Columns 6-9 assume
cluster dissolution by the galactic tidal field only, which provides an upper
estimate on $f_{\rm CC}$.

###  5.1 The upper estimate on $f_{\rm CC}$

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

###  5.2 The lower estimate on $f_{\rm CC}$

In the solar vicinity, only $\approx 10$% of embedded star clusters survive to
the age of $10\,\mathrm{Myr}$, and $\approx 4$% of embedded star clusters
survive to the age of $100\,\mathrm{Myr}$ (Lada & Lada 2003), which is
attributed mainly to encounters with GMCs. The observed fraction of surviving
clusters is smaller than what we obtain in our simulations (c.f. 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

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

We assume that eq. (8) is independent of the cluster mass, and of $R_{\rm g}$.
Then, for each Cepheid in a cluster (which occurs at the cluster age $t$), we
measure the fraction of bound clusters $f_{\rm cl,bound}^{\rm sim}(t)$ in our
simulations. In the real galaxy, there is $f_{\rm cl,bound}^{\rm
obs}(t)/f_{\rm cl,bound}^{\rm sim}(t)$ times Cepheids in clusters at $t$, so
we take $\min(f_{\rm cl,bound}^{\rm obs}(t)/f_{\rm cl,bound}^{\rm 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 column of
Table 2. It is by a factor of $10$ smaller than without the correction (sixth
to ninth column), which indicates that molecular clouds or a similar
disruptive mechanism play likely a crucial role in setting the value of
$f_{\rm CC}$. These values of $f_{\rm CC}$ are comparable to what is observed
in the galaxy M 31 ($f_{\rm 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 $\mathrm{SFE}$ and the gas expulsion
time-scale $\tau_{\rm 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_{\rm 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 released also when $\tau_{\rm M}\gg t_{\rm cr}$.
Nevertheless, lower mass clusters dissolve by $300\,\mathrm{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_{\rm CC}$ to be somewhat larger than in our
models.

On the other hand, if $\mathrm{SFE}\lesssim 1/3$ and $\tau_{\rm M}\lesssim
t_{\rm 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_{\rm CC}$ would be lower than in our models, possibly similar to our
lower estimate on the $f_{\rm CC}$.

The present models started as non-mass segregated, but the lower mass clusters
($M_{\rm ecl}\lesssim 200\,\mathrm{M}_{\odot}$) developed mass segregation
dynamically before the onset of gas expulsion at $t_{\rm
d}=0.6\,\mathrm{Myr}$. However, the value of $t_{\rm d}$ is not well
constrained. If gas expulsion started earlier (later), mass segregation will
develop in clusters less (more) massive than $200\,\mathrm{M}_{\odot}$.

Mass segregation (either primordial or induced by dynamics) impacts $p_{\rm
CC}(M_{\rm 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_{\rm 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_{\rm 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 significant impact on our results. Paper II in this
series studies in detail the dynamics of multiple systems in our simulation.
It suffices here 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_{\rm
CC}$; on the other hand, the dynamical interactions within subsystems and
between subsystems and cluster members are likely to cause more ejections,
decreasing $f_{\rm 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_{\rm cl}$ on a circular orbit of radius $R_{\rm g}$ has a
tidal radius given by

$r_{\rm t}=R_{\rm g}\left(\frac{M_{\rm cl}}{4M_{\rm
glx}}\frac{\gamma^{2}}{(\gamma^{2}-1)}\right)^{1/3},$ (9)

(Binney & Tremaine 2008, their eqs. 3.80 and 8.106), where $M_{\rm glx}$ is
the mass of the galaxy within the cluster orbit and $\gamma=2\Omega/\kappa$ is
the ratio of the orbital and epicyclic frequency. Since $M_{\rm glx}$ follows
from the magnitude of the rotational curve and $\gamma$ follows from its
slope, which is flat for most of spiral galaxies (Bosma 1981; Rubin et al.
1982) implying $\gamma\approx 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 another theories of gravity (e.g. MOND), cluster evaporation can
occur on a different time-scale (faster in the case of MOND) and the tidal
structure around the cluster can be asymmetric (Kroupa et al. in preparation).

This study aims 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_{\rm CC}$ in our models (Table 2) is similar to the
observed fraction of cluster Cepheids in M 31, $f_{\rm CC,M31}=0.025$,
(Senchyna et al. 2015; Anderson & Riess 2018). However, significant variations
in $f_{\rm CC}$ can exist among galaxies, and $f_{\rm CC}\approx 0.07-0.11$
was estimated in the LMC, SMC, and the Milky Way (Anderson & Riess 2018, cf.
their Sect. 3.1).

We here provide a tentative update on the Milky Way value of $f_{\rm CC}$
using the high-quality astrometry from the Gaia mission’s Early Data Release 3
(Gaia Collaboration et al. 2021), updated literature on Galactic open clusters
(Loktin & Popova 2017; Cantat-Gaudin et al. 2018b, 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 limit 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_{\rm CC,MW}=15/217=0.069$, somewhat lower than
the previous value of $11/130=0.085$ in Anderson & Riess (2018). On the one
hand, the decrease in $f_{\rm 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. 2021). We therefore consider our updated estimate of $f_{\rm
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 Sect. 4.3 and Sect. 4.4, decreasing metallicity from $Z=0.014$ to
$Z=0.002$ decreases $f_{\rm CC}$ only by a factor of $1.2$, and increasing
galactocentric radius from $R_{\rm g}=4\,\mathrm{kpc}$ to $R_{\rm
g}=12\,\mathrm{kpc}$ increases $f_{\rm CC}$ also by a factor of $1.2$ 555The
present clusters are calculated for the external gravitational field of the
Galaxy. If the clusters were positioned in the LMC and SMC, their tidal radii
would be only marginally larger by a factor of $1.3$ and $1.1$, respectively,
(the estimate is based on rotation curves from Stanimirović et al. 2004 and
van der Marel & Kallivayalil 2014), than if they were located in the Galaxy at
$R_{\rm g}=8\,\mathrm{kpc}$. Given the modest variation of $f_{\rm CC}$ on
$R_{\rm g}$, we do not expect that the tidal field of the Magellanic Clouds
causes the large difference in the observed value of $f_{\rm CC}$.. A possible
explanation of the difference of $f_{\rm 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 time-
scales, 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 to estimate 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_{\rm
ecl}\lesssim 200\,\mathrm{M}_{\odot}$) occurring before gas expulsion. The
fraction of Cepheids located in clusters increases with the cluster mass from
$\approx 0.25$ (for $M_{\rm ecl}=50\,\mathrm{M}_{\odot}$) to $\approx 0.60$
(for $M_{\rm ecl}=6400\,\mathrm{M}_{\odot}$). Cepheids of higher mass are more
likely (by $\approx 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\,\mathrm{M}_{\odot}$) and the
more massive cluster ($3200\,\mathrm{M}_{\odot}$) is shown as videos in the
online material.

Clusters orbiting the galaxy at a larger galactocentric radius $R_{\rm g}$
take longer to dissolve in the galactic tidal field, retaining more stars (and
thus more ProCeps). From $R_{\rm g}=4\,\mathrm{kpc}$ to $R_{\rm
g}=12\,\mathrm{kpc}$, the fraction of Cepheids in clusters $f_{\rm 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_{\rm 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 to observe two and 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\,\mathrm{M}_{\odot}$ 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 2103 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_{\rm CC}=0.30$ to $0.36$
for the whole galaxy with $Z=0.014$. This is an upper limit on $f_{\rm 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_{\rm CC}$ between the galaxy M 31 ($f_{\rm
CC}\approx 0.025$) on the one handside, and the Galaxy, LMC and SMC ($f_{\rm
CC}\approx 0.09$) on the other handside, differs by a substantially larger
factor than what 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_{\rm CC}=0.015$ to
$0.026$. This disruption is attributed mainly to interactions with giant
molecular clouds, which is not included in the present N-body simulations.
Thus, $f_{\rm CC}$ strongly depends on cluster disruption mechanisms, while
its dependence on the galactic tidal field or metallicity is of secondary
importance.

###### Acknowledgements.

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_{\rm CC}$ based on Gaia EDR3
data. This research made use of Astropy 666http://www.astropy.org, a community
developed core Python package for Astronomy (Astropy Collaboration et al.
2013) and Matplotlib Python Package (Hunter 2007). FD and PK acknowledge the
Scientific Visitors Programme of European Southern Observatory in Garching,
which made this project possible. FD and PK acknowledge support through grant
20-21855S from the Czech Grant Agency. RIA 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

  * Aarseth (1971) Aarseth, S. J. 1971, Ap&SS, 13, 324
  * Aarseth (1999) Aarseth, S. J. 1999, PASP, 111, 1333
  * Aarseth (2003) Aarseth, S. J. 2003, Gravitational N-Body Simulations (Cambridge: Cambridge University Press)
  * Aarseth et al. (1974) Aarseth, S. J., Henon, M., & Wielen, R. 1974, A&A, 37, 183
  * Aarseth & Zare (1974) Aarseth, S. J. & Zare, K. 1974, Celestial Mechanics, 10, 185
  * Ahmad & Cohen (1973) Ahmad, A. & Cohen, L. 1973, Journal of Computational Physics, 12, 389
  * Allen & Santillan (1991) Allen, C. & Santillan, A. 1991, Rev. Mexicana Astron. Astrofis., 22, 255
  * Anderson et al. (2013) Anderson, R. I., Eyer, L., & Mowlavi, N. 2013, MNRAS, 434, 2238
  * Anderson & Riess (2018) Anderson, R. I. & Riess, A. G. 2018, ApJ, 861, 36
  * Anderson et al. (2016) Anderson, R. I., Saio, H., Ekström, S., Georgy, C., & Meynet, G. 2016, A&A, 591, A8
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Banerjee & Kroupa (2013) Banerjee, S. & Kroupa, P. 2013, ApJ, 764, 29
  * Banerjee & Kroupa (2017) Banerjee, S. & Kroupa, P. 2017, A&A, 597, A28
  * Baumgardt et al. (2002) Baumgardt, H., Hut, P., & Heggie, D. C. 2002, MNRAS, 336, 1069
  * Baumgardt & Kroupa (2007) Baumgardt, H. & Kroupa, P. 2007, MNRAS, 380, 1589
  * Baumgardt & Makino (2003) Baumgardt, H. & Makino, J. 2003, MNRAS, 340, 227
  * Bettis & Szebehely (1971) Bettis, D. G. & Szebehely, V. 1971, Ap&SS, 14, 133
  * Bik et al. (2003) Bik, A., Lamers, H. J. G. L. M., Bastian, N., Panagia, N., & Romaniello, M. 2003, A&A, 397, 473
  * Binney & Tremaine (2008) Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press)
  * Bonnell & Davies (1998) Bonnell, I. A. & Davies, M. B. 1998, MNRAS, 295, 691
  * Bono et al. (2005) Bono, G., Marconi, M., Cassisi, S., et al. 2005, ApJ, 621, 966
  * Bosma (1981) Bosma, A. 1981, AJ, 86, 1825
  * Breuval et al. (2020) Breuval, L., Kervella, P., Anderson, R. I., et al. 2020, A&A, 643, A115
  * Brinkmann et al. (2017) Brinkmann, N., Banerjee, S., Motwani, B., & Kroupa, P. 2017, A&A, 600, A49
  * Cantat-Gaudin et al. (2020) Cantat-Gaudin, T., Anders, F., Castro-Ginard, A., et al. 2020, A&A, 640, A1
  * Cantat-Gaudin et al. (2018a) Cantat-Gaudin, T., Jordi, C., Vallenari, A., et al. 2018a, A&A, 618, A93
  * Cantat-Gaudin et al. (2018b) Cantat-Gaudin, T., Jordi, C., Vallenari, A., et al. 2018b, A&A, 618, A93
  * Casertano & Hut (1985) Casertano, S. & Hut, P. 1985, ApJ, 298, 80
  * Castro-Ginard et al. (2020) Castro-Ginard, A., Jordi, C., Luri, X., et al. 2020, A&A, 635, A45
  * Castro-Ginard et al. (2021) Castro-Ginard, A., Jordi, C., Luri, X., et al. 2021, arXiv e-prints, arXiv:2111.01819
  * de la Fuente Marcos & de la Fuente Marcos (2004) de la Fuente Marcos, R. & de la Fuente Marcos, C. 2004, New A, 9, 475
  * Dékány et al. (2015) Dékány, I., Minniti, D., Hajdu, G., et al. 2015, ApJ, 799, L11
  * Dias et al. (2021) Dias, W. S., Monteiro, H., Moitinho, A., et al. 2021, MNRAS, 504, 356
  * Dinnbier & Kroupa (2020) Dinnbier, F. & Kroupa, P. 2020, A&A, 640, A84
  * Eggleton (1983) Eggleton, P. P. 1983, ApJ, 268, 368
  * Evans et al. (2020) Evans, N. R., Günther, H. M., Bond, H. E., et al. 2020, ApJ, 905, 81
  * Ferrarese et al. (2000) Ferrarese, L., Silbermann, N. A., Mould, J. R., et al. 2000, PASP, 112, 177
  * Ferreira et al. (2021) Ferreira, F. A., Corradi, W. J. B., Maia, F. F. S., Angelo, M. S., & Santos, J. F. C., J. 2021, MNRAS, 502, L90
  * Fischer et al. (1992) Fischer, P., Welch, D. L., Cote, P., Mateo, M., & Madore, B. F. 1992, AJ, 103, 857
  * Fujii & Portegies Zwart (2011) Fujii, M. S. & Portegies Zwart, S. 2011, Science, 334, 1380
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1
  * Geyer & Burkert (2001) Geyer, M. P. & Burkert, A. 2001, MNRAS, 323, 988
  * Gieles et al. (2006) Gieles, M., Larsen, S. S., Bastian, N., & Stein, I. T. 2006, A&A, 450, 129
  * Goodwin (1997) Goodwin, S. P. 1997, MNRAS, 284, 785
  * Goodwin & Kroupa (2005) Goodwin, S. P. & Kroupa, P. 2005, A&A, 439, 565
  * Heggie (1975) Heggie, D. C. 1975, MNRAS, 173, 729
  * Holl et al. (2018) Holl, B., Audard, M., Nienartowicz, K., et al. 2018, A&A, 618, A30
  * Hunt & Reffert (2021) Hunt, E. L. & Reffert, S. 2021, A&A, 646, A104
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Hurley et al. (2000) Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543
  * Hurley et al. (2004) Hurley, J. R., Tout, C. A., Aarseth, S. J., & Pols, O. R. 2004, MNRAS, 355, 1207
  * Hurley et al. (2002) Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897
  * Irrgang et al. (2013) Irrgang, A., Wilcox, B., Tucker, E., & Schiefelbein, L. 2013, A&A, 549, A137
  * Irwin (1955) Irwin, J. B. 1955, Monthly Notes of the Astronomical Society of South Africa, 14, 38
  * Jerabkova et al. (2021) Jerabkova, T., Boffin, H. M. J., Beccari, G., et al. 2021, A&A, 647, A137
  * Johnson et al. (2017) Johnson, L. C., Seth, A. C., Dalcanton, J. J., et al. 2017, ApJ, 839, 78
  * Kholopov (1956) Kholopov, P. N. 1956, Peremennye Zvezdy, 11, 325
  * Kochanek (1992) Kochanek, C. S. 1992, ApJ, 385, 604
  * Kroupa (1995a) Kroupa, P. 1995a, MNRAS, 277, 1491
  * Kroupa (1995b) Kroupa, P. 1995b, MNRAS, 277, 1507
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
  * Kroupa et al. (2001) Kroupa, P., Aarseth, S., & Hurley, J. 2001, MNRAS, 321, 699
  * Kuhn et al. (2014) Kuhn, M. A., Feigelson, E. D., Getman, K. V., et al. 2014, ApJ, 787, 107
  * Küpper et al. (2010) Küpper, A. H. W., Kroupa, P., Baumgardt, H., & Heggie, D. C. 2010, MNRAS, 401, 105
  * Küpper et al. (2008) Küpper, A. H. W., MacLeod, A., & Heggie, D. C. 2008, MNRAS, 387, 1248
  * Küpper et al. (2011) Küpper, A. H. W., Maschberger, T., Kroupa, P., & Baumgardt, H. 2011, MNRAS, 417, 2300
  * Kustaanheimo & Stiefel (1965) Kustaanheimo, P. & Stiefel, E. 1965, Reine Angew. Math., 218, 204
  * Lada & Lada (2003) Lada, C. J. & Lada, E. A. 2003, ARA&A, 41, 57
  * Lada et al. (1984) Lada, C. J., Margulis, M., & Dearborn, D. 1984, ApJ, 285, 141
  * Lada et al. (1991) Lada, E. A., Depoy, D. L., Evans, Neal J., I., & Gatley, I. 1991, ApJ, 371, 171
  * Leavitt & Pickering (1912) Leavitt, H. S. & Pickering, E. C. 1912, Harvard College Observatory Circular, 173, 1
  * Lewis et al. (2015) Lewis, A. R., Dolphin, A. E., Dalcanton, J. J., et al. 2015, ApJ, 805, 183
  * Loktin & Popova (2017) Loktin, A. V. & Popova, M. E. 2017, Astrophysical Bulletin, 72, 257
  * Majaess et al. (2013) Majaess, D., Carraro, G., Moni Bidin, C., et al. 2013, A&A, 560, A22
  * Makino (1991) Makino, J. 1991, ApJ, 369, 200
  * Makino & Aarseth (1992) Makino, J. & Aarseth, S. J. 1992, PASJ, 44, 141
  * Mardling & Aarseth (1999) Mardling, R. & Aarseth, S. 1999, in NATO Advanced Study Institute (ASI) Series C, Vol. 522, The Dynamics of Small Bodies in the Solar System, A Major Key to Solar System Studies, ed. B. A. Steves & A. E. Roy, 385
  * Marks & Kroupa (2012) Marks, M. & Kroupa, P. 2012, A&A, 543, A8
  * McMillan et al. (2007) McMillan, S. L. W., Vesperini, E., & Portegies Zwart, S. F. 2007, ApJ, 655, L45
  * Medina et al. (2021) Medina, G. E., Lemasle, B., & Grebel, E. K. 2021, MNRAS, 505, 1342
  * Megeath et al. (2016) Megeath, S. T., Gutermuth, R., Muzerolle, J., et al. 2016, AJ, 151, 5
  * Meingast & Alves (2019) Meingast, S. & Alves, J. 2019, A&A, 621, L3
  * Mikkola & Aarseth (1990) Mikkola, S. & Aarseth, S. J. 1990, Celestial Mechanics and Dynamical Astronomy, 47, 375
  * Miyamoto & Nagai (1975) Miyamoto, M. & Nagai, R. 1975, PASJ, 27, 533
  * Moe & Di Stefano (2017) Moe, M. & Di Stefano, R. 2017, ApJS, 230, 15
  * Mouri & Taniguchi (2002) Mouri, H. & Taniguchi, Y. 2002, ApJ, 580, 844
  * Musella et al. (2016) Musella, I., Marconi, M., Stetson, P. B., et al. 2016, MNRAS, 457, 3084
  * Oh & Kroupa (2016) Oh, S. & Kroupa, P. 2016, A&A, 590, A107
  * Oh et al. (2015) Oh, S., Kroupa, P., & Pflamm-Altenburg, J. 2015, ApJ, 805, 92
  * Pawlak et al. (2019) Pawlak, M., Pejcha, O., Jakubčík, P., et al. 2019, MNRAS, 487, 5932
  * Perets & Šubr (2012) Perets, H. B. & Šubr, L. 2012, ApJ, 751, 133
  * Pflamm-Altenburg et al. (2013) Pflamm-Altenburg, J., González-Lópezlira, R. A., & Kroupa, P. 2013, MNRAS, 435, 2604
  * Pflamm-Altenburg & Kroupa (2006) Pflamm-Altenburg, J. & Kroupa, P. 2006, MNRAS, 373, 295
  * Pflamm-Altenburg & Kroupa (2008) Pflamm-Altenburg, J. & Kroupa, P. 2008, Nature, 455, 641
  * Porras et al. (2003) Porras, A., Christopher, M., Allen, L., et al. 2003, AJ, 126, 1916
  * Pouliasis et al. (2017) Pouliasis, E., Di Matteo, P., & Haywood, M. 2017, A&A, 598, A66
  * Riess et al. (2020) Riess, A. G., Yuan, W., Casertano, S., Macri, L. M., & Scolnic, D. 2020, ApJ, 896, L43
  * Ripepi et al. (2019) Ripepi, V., Molinaro, R., Musella, I., et al. 2019, A&A, 625, A14
  * Röser et al. (2019) Röser, S., Schilbach, E., & Goldman, B. 2019, A&A, 621, L2
  * Rubin et al. (1982) Rubin, V. C., Ford, W. K., J., Thonnard, N., & Burstein, D. 1982, ApJ, 261, 439
  * Samus’ et al. (2017) Samus’, N. N., Kazarovets, E. V., Durlevich, O. V., Kireeva, N. N., & Pastukhova, E. N. 2017, Astronomy Reports, 61, 80
  * Sana et al. (2012) Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444
  * Sandage (1958) Sandage, A. 1958, ApJ, 128, 150
  * Scargle et al. (2013) Scargle, J. D., Norris, J. P., Jackson, B., & Chiang, J. 2013, ApJ, 764, 167
  * Schoettler et al. (2020) Schoettler, C., de Bruijne, J., Vaher, E., & Parker, R. J. 2020, MNRAS, 495, 3104
  * Senchyna et al. (2015) Senchyna, P., Johnson, L. C., Dalcanton, J. J., et al. 2015, ApJ, 813, 31
  * Sollima (2020) Sollima, A. 2020, MNRAS, 495, 2222
  * Soszyński et al. (2020) Soszyński, I., Udalski, A., Szymański, M. K., et al. 2020, Acta Astron., 70, 101
  * Spera et al. (2016) Spera, M., Mapelli, M., & Jeffries, R. D. 2016, MNRAS, 460, 317
  * Spitzer (1958) Spitzer, Lyman, J. 1958, ApJ, 127, 17
  * Spitzer (1969) Spitzer, Jr., L. 1969, ApJ, 158, L139
  * Stanimirović et al. (2004) Stanimirović, S., Staveley-Smith, L., & Jones, P. A. 2004, ApJ, 604, 176
  * Sterzik & Durisen (1998) Sterzik, M. F. & Durisen, R. H. 1998, A&A, 339, 95
  * Stiefel & Scheifele (1975) Stiefel, E. L. & Scheifele, G. 1975, Linear and regular celestial mechanics. Perturbed two-body motion. Numerical methods. Canonical theory.
  * Tanikawa et al. (2012) Tanikawa, A., Hut, P., & Makino, J. 2012, New A, 17, 272
  * Terlevich (1987) Terlevich, E. 1987, MNRAS, 224, 193
  * Testa et al. (2007) Testa, V., Marconi, M., Musella, I., et al. 2007, A&A, 462, 599
  * Theuns (1992) Theuns, T. 1992, A&A, 259, 493
  * Tout et al. (1996) Tout, C. A., Pols, O. R., Eggleton, P. P., & Han, Z. 1996, MNRAS, 281, 257
  * Turner (1976) Turner, D. G. 1976, AJ, 81, 1125
  * Turner & Burke (2002) Turner, D. G. & Burke, J. F. 2002, AJ, 124, 2931
  * Tutukov (1978) Tutukov, A. V. 1978, A&A, 70, 57
  * van der Marel & Kallivayalil (2014) van der Marel, R. P. & Kallivayalil, N. 2014, ApJ, 781, 121
  * Vesperini & Heggie (1997) Vesperini, E. & Heggie, D. C. 1997, MNRAS, 289, 898
  * Wang et al. (2019) Wang, L., Kroupa, P., & Jerabkova, T. 2019, MNRAS, 484, 1843
  * Watson et al. (2006) Watson, C. L., Henden, A. A., & Price, A. 2006, Society for Astronomical Sciences Annual Symposium, 25, 47
  * Weidner et al. (2010) Weidner, C., Kroupa, P., & Bonnell, I. A. D. 2010, MNRAS, 401, 275
  * Weidner et al. (2004) Weidner, C., Kroupa, P., & Larsen, S. S. 2004, MNRAS, 350, 1503
  * Welch & Stetson (1993) Welch, D. L. & Stetson, P. B. 1993, AJ, 105, 1813
  * Whitmore et al. (1999) Whitmore, B. C., Zhang, Q., Leitherer, C., et al. 1999, AJ, 118, 1551
  * Zhou & Chen (2021) Zhou, X. & Chen, X. 2021, MNRAS, 504, 4768

##  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 Fig. 13.

The left panel details the evolution of a $11\,\mathrm{M}_{\odot}$ and
$7\,\mathrm{M}_{\odot}$ 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\,\mathrm{M}_{\odot}$. Core He burning (orange points) then commences
leftwards from 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 detection of Cepheids (Sect. 3).

The middle panel shows a $8.4\,\mathrm{M}_{\odot}$ primary with a secondary of
$1.1\,\mathrm{M}_{\odot}$ and initial orbital period of $13$ days. In this
case, the coalescence increases the mass only slightly, and the ProCep becomes
a Cepheid.

The right panel shows coalescence of a primary of ZAMS mass
$4.3\,\mathrm{M}_{\odot}$ with its $1.6\,\mathrm{M}_{\odot}$ companion. For
metallicity $Z=0.014$, the initial mass of the primary is lower than the mass
threshold for Cepheids, which is $m_{\rm min,Ceph}=4.7\,\mathrm{M}_{\odot}$.
After the coalescence, the merged star has a mass of
$5.3\,\mathrm{M}_{\odot}$, and it becomes a Cepheid.

The cases in the middle and right panel 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.

Figure 13:  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\,\mathrm{M}_{\odot}$), 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\,\mathrm{M}_{\odot}$), which would not
have evolved to a Cepheid in isolation, becomes a Cepheid after a coalescence.

##  Appendix B Comparison of the evolution of ProCeps in a lower mass
($100\,\mathrm{M}_{\odot}$) and a more massive ($3200\,\mathrm{M}_{\odot}$)
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\,\mathrm{Myr}$
(upper row of Fig. 14). At this time (upper middle panel), two ProCeps are
already near the cluster centre. The third ProCep, which is located at
$x\approx-0.6\,\mathrm{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 escaping 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\,\mathrm{Myr}$ (shown only in the online
material).

In contrast, the more massive star cluster (lower row of Fig. 14) does not
mass segregate ProCeps before gas expulsion, which results in the population
of ProCeps to be 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.

Figure 14:  Evolutionary snapshots for an example of a lower mass
($100\,\mathrm{M}_{\odot}$; upper row) and more massive star clusters
($3200\,\mathrm{M}_{\odot}$; 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 $m_{\rm
min,Ceph}$ or above $m_{\rm max,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. Note
that the length-scale on the right column is inflated by a factor of $10$. The
evolution of both star clusters up to $12\,\mathrm{Myr}$ can be seen as movies
in the online material.