© The Authors 2024

1 Introduction

Classical Cepheids (henceforth, Cepheids) are evolved intermediate-mass stars best known for their high-amplitude pulsations and application as standard candles, thanks to Leavitt’s law (Leavitt & Pickering 1912, period-luminosity relation). The B-star progenitors of Cepheids feature a frequency of companions greater than unity and recent astrometric studies based on the ESA mission Gaia (Gaia Collaboration 2016, 2018, 2021) imply a Cepheid binary fraction on the order of 80% (Kervella et al. 2019b,a). By including other types of binaries (UV observations, X-ray observations), Evans et al. (2022) have provided a comprehensive lower estimate of the binary fraction to be 57 ± 12%. Qualitatively, it is clear that short-period binaries are influenced (mass transfer or common envelope evolution) by the secular evolution of the primary star towards the red giant phase, which precedes the core He-burning blue loop phase where the vast majority of Cepheids are observed. However, observational selection effects related to the generally large distances of Cepheids (hundreds to thousands of pc), long orbital periods (≳1 yr), and high contrast ratios at optical and near infrared wavelengths are challenging with respect to devising a detailed characterization of the Cepheid companion population and its orbital characteristics.

Substantial efforts to observationally determine the binary fraction of Cepheids have been made, resulting in a large body of accessible literature on the subject (cf. the binary Cepheids database Szabados 2003). For the sake of concision, it is sufficient to mention only a few of the most relevant recent studies based on high-quality data: line-of-sight evidence for orbital motion comes primarily from radial velocity (RV) observations (cf. Anderson 2018 for a recent review), and, in a couple of cases, from pulsation timing (Csörnyei et al. 2022). Evans et al. (2015) determined a magnitude-limited binary fraction for stars with orbital periods ≲20 yr of 29 ± 8% based on 40 Milky Way Cepheids brighter than a magnitude of 9. The systematic search for spectroscopic binary Cepheids (Shetye et al. 2024) based on data from the VELOCE project (Anderson et al. 2024) has reassessed the evidence for binary candidates and identified many new binaries. This has sharpened the spectroscopic binary fraction to 29.6 ± 3.4%, with 15.2 ± 2.4% having orbital periods shorter than 10 yr. Shetye et al. (2024) further found the shortest-period binary Cepheid in the Milky Way, R Crucis, with a period of 237 d, provided 30 homogeneously determined orbital solutions, and identified 17 triple systems with distant, visual companions of unknown period.

Astrometry offers a completely complementary and very detailed view of Cepheid multiplicity. For example, Kervella et al. (2019b) used varying proper motion vectors across the epochs from HIPPARCOS and Gaia DR2 to infer the presence of unseen companions and determine a binary fraction of ≳80% for Cepheids, more consistent with the high multiplicity fraction of B-stars. In a companion paper, Kervella et al. (2019a) further identified likely wide (resolved) companions via proper motion pairs.

The search for resolved companions requires techniques of resolving very fine angular scales due to the generally large distances of Cepheids (hundreds of pc to kpc). From the ground, the most successful techniques employed include lucky imaging (Gallenne et al. 2014) and long-baseline interferometry (e.g. Kervella et al. 2019a), which have allowed detection of companions with H-band magnitude difference of up to 5.5. From space, HST’s WFC3 camera has been used to identify close-in (<5″ separation) companions of a magnitude-limited sample of 70 Cepheids and revealed that all the 7 Cepheids with resolved companions closer than 2″ were also spectroscopic binaries, often with significant eccentricity (Evans et al. 2020). A companion survey of more distant companions (>5″ separation) (Evans et al. 2016) matched with XMM-Newton observations to identify young companion stars found that X-ray active (i.e. young) companions of the 14 Cepheids with available X-ray observations were all closer than 4000 AU. More distant resolved companions identified in the vicinity of δ Cephei (which is also a low-amplitude spectroscopic binary; cf. Anderson et al. 2015) and S Normae (Irwin 1955) were explained as coeval stars originating in the same birth cluster or association.

The properties of companions to Cepheids and their multiplicity frequency reflect the history of their evolution since their birth as mid-B stars within clusters embedded in molecular clouds (Lada & Lada 2003; Zinnecker & Yorke 2007), followed by cloud dispersal due to stellar feedback and the dynamical evolution of the exposed cluster in the galactic tidal field. This includes stellar encounters with other stars inside the natal cluster and stellar evolution of the binary components, whose importance increases when the stars evolve beyond the main sequence (MS). The influence of stellar evolution was estimated by Neilson et al. (2015) by means of binary star population synthesis, who found that the presence of a companion leads to a lower limit in Cepheid orbital periods of approximately one year and that the destruction of closer-in companions implies that the mid-B star binary fractions should be about 15–30% larger than the measured Cepheid binary fraction. More recently, Karczmarek et al. (2022) investigated the companion properties of a population of Cepheids in environments of different metallicity.

However, the influence of the dynamical environment of the Cepheid birth cluster has not been investigated yet. Based on theoretical estimates and general results of N-body modelling, it has been firmly established that stellar encounters of binaries with other single and binary stars in the cluster can change the orbital period and eccentricity of the binary by hardening, softening and dissolution, according to Heggie’s law (Heggie 1975; Hills 1975; Heggie & Hut 2003). These interactions also involve resonances, exchanges, or collisions (e.g. Saslaw et al. 1974; Hut & Bahcall 1983; Sigurdsson & Phinney 1993; Fregeau et al. 2004; Ginat & Perets 2021) and they often result in fast ejections of single or binary stars, which are observed after their escape from the cluster as walkaway (de Mink et al. 2014; Schoettler et al. 2019) or runaway stars (e.g. Leonard & Duncan 1988, 1990; Hoogerwerf et al. 2001; Eldridge et al. 2011; Fujii & Portegies Zwart 2011; Maíz Apellániz et al. 2018). The ability of a cluster to eject stars increases with the cluster density (Hut & Bahcall 1983; Binney & Tremaine 2008; Oh et al. 2015) and with the presence of primordial binaries (Perets & Šubr 2012; Oh et al. 2015; Oh & Kroupa 2016). Some of the binaries containing Cepheids are likely to capture another star as an outer companion, thereby forming a hierarchical triple. The outer orbits in triple stars are able to influence the orbital parameters of the inner binary via Kozai-Lidov oscillations (Kozai 1962; Lidov 1962), bringing the binary components closer together at their pericentre, whereby the orbit is circularised by tidal forces. This mechanism produces short-period, low-eccentricity inner binaries (Kiseleva et al. 1998; Fabrycky & Tremaine 2007).

In Paper I of this series (Dinnbier et al. 2022), we investigated the occurrence of Cepheids in star clusters and the field using N-body simulations. Here, we analyse the same suite of simulations with respect to the properties of Cepheid multiplicity. Stellar evolution in NBODY6, in particular with respect to binary evolution (e.g. tidal circularisation, common envelope phase), has been implemented as described by Hurley et al. (2000) and Hurley et al. (2002), which also form the basis of stellar evolution calculations for the analysis due to Neilson et al. (2015) and Karczmarek et al. (2022). A key difference between these two works and our study is that our N-body simulations consider the full cluster dynamical picture for a large range of birth environments, spanning from large bound clusters to field Cepheids mostly originating from dispersed low-mass clusters. For example, this also allows us to assess the occurrence frequency of triples and higher-order systems that form dynamically after the zero-age main sequence (ZAMS).

2 Initial conditions for binary stars

This work utilises the extensive library of simulations documented in Paper I, which also describes the initial conditions for star clusters and numerical methods used in the code NBODY6 (Aarseth 1999, 2003). Here, we provide a detailed description only of the initial conditions for the binary population and its subsequent integration, as these features are key in the context of the present work.

All stars are formed in binaries at the beginning of the simulations. This means that all triple and higher-order systems in our simulations are formed only by capture events. Single stars, which are later paired to binaries, are drawn from the initial mass function (IMF) of Kroupa (2001). To reflect the observed differences in the distribution functions for the orbital period, eccentricity, and mass ratio between early-type and late-type stars (Moe & Di Stefano 2017; Offner et al. 2023), we used two types of the initial distributions for binaries according to their primary mass, m₁; below and above the primary threshold mass, m_1,thr = 3 M_⊙. The choice of m_1,thr, which is slightly lower than the usual value of ≈5 M_⊙ (Kroupa et al. 2013), ensures that all stars that evolve into Cepheids (minimum mass of 3.2 M_⊙ for the considered metallicity range) have their orbital parameters and mass ratios drawn from the same distribution functions (i.e. those for early-type stars), and they avoid the sudden transition to the other set of distribution functions, which occurs at m_1,thr. The initial distributions for orbital periods for binaries with m₁ < m_1.thr and m₁ > m_1,thr are shown in Fig. 1 by the solid black and magenta lines, respectively.

Binaries with the primary of mass m₁ < m_1,thr are generated with the birth distribution of orbital periods and eccentricities according to Kroupa (1995a). The interaction of the two protostars with each other and with the protostellar material that remains in the system (pre-MS eigenevolution) is approximated by the method of Kroupa (1995b) with the update attributed to Belloni et al. (2017). Pre-MS eigenevolution takes place on the time-scale of 10⁵ years, transforming the birth binary population to the initial binary population. Pre-main-sequence (pre-MS) eigenevolution is applied on each binary individually before the start of the simulations; the initial binary population is then taken as the initial conditions for the integration in NBODY6, which starts at time t₀. For simplicity, it is assumed that all stars are already at their ZAMS at a time, t₀; whereas, in nature, it takes tens of Myr for solar type stars to reach the ZAMS (for the dynamical consequences of this approximation, we refer to Railton et al. 2014).

The birth distribution of binaries is a theoretical concept (analogically as the IMF), which is not likely to be observed directly because of the large extinction and rapid evolution of these systems. However, the adopted form of the birth distribution, after being processed via pre-MS eigenevolution and then cluster dynamics (stimulated evolution), closely reproduces the observed properties (binarity, orbital period distribution) of stars in clusters and in the field (e.g. Kroupa 1995b; Kroupa et al. 2001; Marks & Kroupa 2012; Belloni et al. 2018), which is accessible to observations.

Binaries with a primary of mass m₁ > m_1,thr were generated with the birth distribution of the mass ratio, orbital periods, and eccentricity taken from Belloni et al. (2017), which differs from that of their lower mass counterparts. These distributions are based mainly on the observations due to Sana et al. (2012); Kobulnicky et al. (2014); Sana et al. (2014); Aldoretta et al. (2015). In particular, the orbital period distribution function, f_p, for massive stars is given by (Sana et al. 2012): $$f_{\mathrm{p}}({\log }_{10}P)=0.23{({\log }_{10}P)}^{-0.55},$$(1)

where orbital periods occupy the range log₁₀ P ∈ (0.15, 6.7), where P is in days. In contrast to the lower mass binaries, the more massive binaries have no pre-MS eigenevolution, namely, their birth distribution is the same as the initial distribution. This approximation is motivated by the much faster pre-MS evolution of massive stars (e.g. Stahler & Palla 2004), uncertainties in massive star formation (e.g. Zinnecker & Yorke 2007; Peters et al. 2010), and the fact that a study addressing the influence of cluster environment on properties of massive binaries (an analogy to the iterative solution of Kroupa 1995a for late type binaries) has not been undertaken so far. Instead, we assume that massive binaries form with the same properties as the observed ones.

When initialising the binary population, it is of paramount importance (as stressed in Kroupa & Jerabkova 2018) to not affect the IMF. This is achieved by first sampling stars from the stellar IMF and only then distributing the already sampled stars into the binaries. We note that initializing the binary population by first sampling only the primaries from the IMF and then generating the secondary mass according to an assumed mass-ratio distribution (as is often done) is likely to affect the shape of the IMF, which can lead to contradictions with the star counts.

The orbital parameters of binaries change not only due to the stellar evolution in binary stars, but also due to dynamical encounters with other stars and binaries within the cluster. According to the Heggie’s law (Heggie 1975; Binney & Tremaine 2008), sufficiently close binaries harden (reduce their orbital period P) due to dynamical encounters, bringing their companion closer so that the binary merges in some cases during post-MS evolution, although it would not have merged had it evolved in the absence of the other stars in the same cluster. A sufficiently large mass increase during such a merger might increase the mass of the star so much that it does not become a Cepheid. In contrast, wider binaries can soften with the opposite consequence of being dissociated. These examples point out that the fraction of stars which actually become Cepheids is given not only by the distribution function of their initial orbital periods, but also by the internal stellar dynamics of their host clusters.

Here, the binaries are investigated in star clusters orbiting the galaxy on circular trajectories at galactocentric radii R_g = 4 kpc, R_g = 8 kpc, and R_g = 12 kpc for a metallicity of Z = 0.014. Positioning the cluster orbits at different radii results in a different degree of confinement due to the gravitational potential of the galaxy (i.e. different tidal radius of the cluster), which is represented by the Milky-Way like gravitational potential (Allen & Santillan 1991). At each galactocentric radius, we have studied clusters across a wide mass spectrum: from M_ecl = 50 M_⊙ to 6400 M_⊙ (M_ecl is the total stellar mass of the initially gas embedded cluster) in steps increasing M_ecl by a factor of 2. Less massive clusters are realised more often than the more massive ones to obtain more robust statistics for Cepheids; clusters of mass 6400 × 2⁻ⁱ M_⊙ are realised 2ⁱ⁺¹ times (i being an integer from 0 to 7) with a different random seed for initial conditions. We also studied the possible influence of metallicity by calculating additional sets of models at R_g = 8 kpc with Z = 0.002 and Z = 0.006. These models, which all include internal cluster dynamics, are referred to as “standard models”. We refer to Paper I for more details.

To separate the influence of stellar evolution and cluster dynamics on binaries, we performed additional simulations with cluster radii large enough (half mass radius ≈100 pc) for stellar dynamical processing of the binary population to be negligible, but with a stellar evolution identical to that in the standard models. These models, which we call “control models”, are calculated for the same metallicity as the standard models (Z = 0.014, Z = 0.006, and Z = 0.002). Control models consist of the same cluster masses M_ecl as standard models (from 50 M_⊙ to 6400 M_⊙ in steps increasing M_ecl by a factor of 2) and with the same number of realisations (e.g. we calculated two times as many clusters of initial mass M_ecl than clusters of initial mass 2 M_ecl). The reason of using the same cluster mass in control models is to generate the initial binary population of identical properties as in the standard models; the least massive clusters (M_ecl ≲ 200 M_⊙) do not form the most massive Cepheids because we assume that the mass of the most massive star is a function of the cluster mass (Weidner et al. 2010; Yan et al. 2023; see also Table 1 of Paper I for maximum stellar mass, m_max, forming in a cluster of mass M_ecl). Control simulations follow the evolution not only for Cepheids that evolve from stars of the appropriate initial mass range (m_min,Ceph, m_max,Ceph), but also for Cepheids which originate from a star below this mass range, which gains mass by merging with its companion. For this purpose, the control runs are calculated up to 2 Gyr, which is long enough for the least massive stars which can evolve to Cepheids in a mass conserving merger of two m_min,Ceph/2 ≈ 1.6 M_⊙ stars.

The typical mass of stars that evolve to Cepheids (mainly mid-B stars) is lower than the mass of O stars, whose binary period distribution was found by Sana et al. (2012) (c.f. Eq. (1)). Generally, observations show that the median of the orbital period distribution of early-type stars increases with decreasing primary mass (Duchêne & Kraus 2013; Moe & Di Stefano 2017; Offner et al. 2023; Fig. 1), so the stars that evolve into Cepheids have a wider orbital period distribution than O-stars; the assumption of the same orbital period distribution for all stars with m > 3 M_⊙ according to Eq. (1) is a simplification. To access the influence of this approximation on the binary evolution of stars with 3 M_⊙ ≲ m ≲ 12 M_⊙, we perform an additional control run of binaries with a wider orbital period distribution¹, $$f_{\mathrm{p}}({\log }_{10}P)=0.23{({\log }_{10}P+\frac{3}{2}{\log }_{10}(A))}^{-0.55},$$(2)

with the parameter A = 3. The purpose of this distribution is to approximate the orbital period distribution for stars in the (2 M_⊙, 5 M_⊙) range (cyan line in Fig. 1). According to the modified Sana’s distribution with A = 3, 50% of binaries have log₁₀(P[days]) < 3.5, while according to the original Sana distribution Eq. (1), 50% of binaries have log₁₀(P[days]) < 2.1. Models with binary period distribution of Eq. (2) are investigated in Sec. 4.4.

During the simulation, NBODY6 regularises compact binaries, and integrates them according to the Kustaanheimo-Stiefel method (Kustaanheimo & Stiefel 1965). Hierarchical systems, which form dynamically in the course of the simulation, are treated either by the three-body regularization method of Aarseth & Zare (1974), if the inner binary is perturbed by the outer body, or are treated as single objects (Mardling & Aarseth 1999) if the inner binary is sufficiently isolated from the outer body, which considerably speeds up the calculation. NBODY6 takes into account single star (Tout et al. 1996; Hurley et al. 2000) as well as binary star evolution (Aarseth 1996; Hurley et al. 2002), which is approximated by synthetic stellar evolutionary tracks. Collisions and subsequent merges between stars are evaluated according to the criterion of Kochanek (1992). The code enables modelling of the common envelope evolution and Roche-lobe mass transfer, where the Roche radius is estimated by the formula of Eggleton (1983). The compact remnant after a supernova explosion receives a kick in a randomly selected direction and of the velocity drawn from a Maxwell distribution with 1D dispersion σ = 40 km s⁻¹ (Aarseth & Mardling 1996; Hansen & Phinney 1997). The kick might result in binary disruption.

Fig. 1 Cumulative distribution functions of orbital periods of binary stars according to the primary mass. The observed distributions for primaries in specific mass intervals are shown by bars (Moe & Di Stefano 2017): O-type stars (m > 16 M⊙; blue bars), mid B-type stars (5 M⊙ < m < 9 M⊙; green), late B- and A-type stars (2 M⊙ < m < 5 M⊙; yellow), and solar type stars (0.8 M⊙ < m < 1.2 M⊙; red). The theoretical birth distribution of Kroupa (1995a) for late type stars (m < 3 M⊙; dashed black line) is transformed by pre-MS eigenevolution to the initial distribution (solid black line). The initial distribution for O-stars according to Eq. (1) (Sana et al. 2012) and the wider distribution of Eq. (2) with A = 3, which is the widest estimate on the Cepheid orbital periods, is plotted by the magenta and cyan line, respectively.

3 Evolutionary channels of Cepheids in binaries

For the sake of clarity, we introduce terminology for Cepheids originating from stars of different initial masses or with a different history of interaction with their companions. Several terms of the present terminology differ from the terminology adopted in Paper I; to avoid confusion, we present a comparison of the terminologies in Table A.1. We refer to the stars formed within the mass range from which Cepheids evolve in isolation as prospective Cepheids (hereafter, ProCeps). The ZAMS mass of ProCeps slightly increases with metallicity from the range (m_min,Ceph, m_max,Ceph) ∈ (3.2, 11) M_⊙ for Z = 0.002 to (m_min,Ceph, m_max,Ceph) ∈ (4.7, 13) M_⊙ for Z = 0.014 (Anderson et al. 2016b). Given that a significant fraction of ProCeps are in close binaries Eq. (1), these stars interact with their companions as they evolve out of the MS. Close binaries (log₁₀(P[days]) ≲2.6; Neilson et al. 2015; Moe & Di Stefano 2017) merge or undergo Roche-lobe mass transfer or common envelope evolution, which in some cases changes the evolution of the final object to such an extent that it does not become a Cepheid or it becomes a Cepheid at a significantly different time than expected from its mass at the ZAMS.

The considered mass changes are sketched in Fig. 2, where ProCeps are located between the blue-dashed horizontal lines (mass is on the ordinate) at their ZAMS. Some ProCeps evolve to the Cepheid stage, which is shown by the hatched area. Some of the ProCeps that evolve to Cepheids do not interact with their companions, they effectively evolve in isolation (CanCeps; violet line), while others experience mass gain (InterCeps-I; grey line) or mass loss (InterCeps-II; grey line) to their companions. The mass gain or mass loss can be due to merger or mass transfer with their companion. The Cepheids that evolve from ProCeps (i.e. CanCeps and InterCeps of both types) are refereed to as Cepheids originating from the expected mass range: RanCeps. The more massive Cepheids have the tips of their blue loops outside the instability strip (IS).

Other ProCeps experience more significant interactions with their companions (e.g. merger, Roche-lobe mass transfer or common envelope evolution) that prevent them from reaching the Cepheid stage (ExCeps). If the ExCep gains mass during the interaction, it is defined to be of type I (brown line); if the ExCep loses mass during the interaction, it is defined to be of type II (orange line). Analogously, some stars with the ZAMS mass outside the interval (m_min,Ceph, m_max,Ceph) become Cepheids (AdCeps) either due to mass gain (type I; dark green) or mass loss (type II; light green) with respect to their companion. The method of automatic identification of Cepheids in NBODY6 simulations is described in sect. 3 of Paper I.

Before analysing the process of binary evolution leading to Cepheids, it is instructive to estimate the percentage of ProCeps that form as primary or secondary in a representative sample of a coeval stellar population. Assuming the IMF of Kroupa (2001), the mass ratio distribution between the primary and secondary of Belloni et al. (2017) and a star-burst producing embedded clusters distributed according to the mass function of Eq. (4) with the cluster mass range spanning from 50 M_⊙ to 6400 M_⊙, approximately 10% of ProCeps are secondaries where the primary is more massive than m_max,Ceph, 50% of ProCeps are in a binary with another ProCep, and 40% of ProCeps are primaries where the secondary is less massive than m_min.Ceph.

Fig. 2 Schematic illustration of the considered mass changes due to binary evolution. The abscissa represents time; the ordinate represents the initial stellar mass, m. ProCeps, which have an initial mass range of (mmin,Ceph, mmax,Ceph) evolve towards Cepheids (without mass exchange: violet line or with mass exchange: grey lines), the state of which is marked by the hatched area. The beginning of core helium burning is indicated by the yellow dotted line. See Sect. 3 for full explanation of the evolutionary channels and the terminology used in this paper.

4 Evolution of binary and multiple Cepheids in star clusters

In star clusters, binary stars interact with the other cluster members, which can change the initial orbital parameters of a binary before it starts evolving out of the MS. Depending on the nature of the interaction, the binary can harden (a decrease of orbital period P), soften (an increase of P), or be ionised (disrupted). For an idealised star cluster of half-mass radius r_h and consisting of N equal mass stars, it is possible to derive an approximate watershed for the binary semi-major axis a_hard (Heggie & Hut 2003): $$a_{\text{hard }}=r_{\mathrm{h}}/N\text{, }$$(3)

where binaries with the semi major axis smaller than a_hard harden and binaries with the semi major axis larger than a_hard soften. For comparison, in the control sets of simulations, binaries are subjected to their stellar (binary) evolution only, being isolated from the possible influence of other binaries or single stars.

When the quantities of interest are evaluated for a large sample of stars formed in a population of star clusters of different masses, we assume that the embedded cluster mass function (ECMF) takes the form $$\frac{\mathrm{d}{N}_{\mathrm{e}\mathrm{c}\mathrm{l}}\left(M_{\mathrm{e}\mathrm{c}\mathrm{l}}\right)}{\mathrm{d}{M}_{\mathrm{e}\mathrm{c}\mathrm{l}}}\propto M_{\mathrm{e}\mathrm{c}\mathrm{l}}^{-\beta },$$(4)

where M_ecl is the initial stellar mass of the embedded cluster and dN_ecl is the number of clusters in the cluster mass bin of size dM_ecl. We adopt the slope β = 2 (Lada & Lada 2003; Bik et al. 2003; de la Fuente Marcos & de la Fuente Marcos 2004; Gieles et al. 2006). The mass of the embedded clusters spans the range of (50 M_⊙, 6400 M_⊙). The typical number of Cepheids formed in the population of clusters in one set of models is ≈1700 (for Z = 0.014; the number increases for models with sub-solar metallicities). The ECMF mass slope and the range of cluster masses are the same as adopted in Paper I.

4.1 Cepheids in binaries

We consider a Cepheid to be a member of a binary system (hereafter binary Cepheid; BiCep) if it is a component of a binary of an orbital period shorter than 10¹⁰ days (corresponding to a semi major axis of 1 pc and orbital period of 28 Myr for an 8 M_⊙ binary) both at the terminal-age main sequence (TAMS) and at the end of core helium burning (t_He,end); the latter event occurs immediately after the Cepheid stage. This constitutes rather an upper limit on the binary orbital period because the orbital time is comparable to the time for a star to evolve to a Cepheid and also the Jacobi radius of the binary (≈3 pc) in the Galactic potential is not much larger than its semi-major axis. Although some of the widest binaries are only temporarily bound, they do not increase the binary fraction significantly because most of the binaries have log₁₀(P[days]) ≲ 8 (Fig. 3).

4.1.1 Fraction of binary Cepheids in clusters and in the field

We adopt the definition of the fraction of binary Cepheids p_BC as the ratio of all Cepheids with at least one companion to the total number of Cepheids². We determine the fraction of BiCeps separately for Cepheids located in star clusters (fraction p_BC,in)³ and in field (p_BC,out). For a population of stars originating from clusters of the ECMF of Eq. (4), p_BC,in ≈ 0.61 is larger than p_BC,out ≈ 0.36 (Table 1). The binary fraction both in the clusters and in the field decreases with increasing cluster mass as shown in the left panel of Fig. 4. The influence of pure binary evolution is obtained from the control models; the coalescence with companions (in 44% of cases) and the disruption due to the supernova kick to primaries where Cepheids are secondaries (in 6.5% of cases), reduce the fraction of BiCeps to p_BC = 0.50 (for Z = 0.014; Table 1). Cluster dynamics slightly decreases the fraction of BiCeps, because the standard models have p_BC ≈ 0.42 (for Z = 0.014) in clusters and field combined.

Fig. 3 Cumulative distribution of orbital periods of ProCeps (taken at t0; black solid line) and BiCeps (taken at tHe,end; dashed lines). The orbital period distribution of BiCeps is shown separately for Cepheids in the lower (magenta dashed line) and more massive clusters (cyan dashed line), in the whole cluster population (yellow dashed line), in the field originating from the whole cluster population (black dashed line), and in the control model (brown dashed line). The orbital period distribution of Cepheids with a compact companion is indicated by the black dotted line. The distributions are normalised to the total number of BiCeps in each group. There are almost no BiCeps with log10(P[days]) ≲ 2 because of coalescence with their companion before the onset of core helium burning. The vertical dotted line represents the upper limit on the initial orbital periods (log10(P[days]) ≲ 6.7). All BiCeps with an orbital period longer than that formed due to dynamical encounters in clusters.

4.1.2 Changes in orbital period distribution

In order to obtain the fraction of binaries which soften (p_soft) or harden (p_hard), we compare the orbital periods between time t₀ and the terminal age main-sequence t_TAMS of the binaries in which the ProCeps reside. We choose the TAMS as the upper time limit because it is the time when the star starts increasing its radius substantially, which increases the probability of the interaction with its companion, so the evolution before the TAMS is mostly influenced by cluster dynamics (apart from the ProCeps which are secondaries), while the evolution after TAMS is mostly influenced by stellar evolution and possible mass transfer. The time interval from t₀ to t_TAMS is also substantially longer than the interval from t_TAMS to the end of core helium burning t_He,end, so any change of orbital parameters due to cluster dynamics is more likely to occur between time t₀ and t_TAMS than between t_TAMS and t_He,end.

For softening, we require that the orbital period of the binary either increases by at least a factor of 2 between time t₀ and t_TAMS, or the binary gets dissociated. For hardening, we require that the orbital period of the binary either decreases by at least a factor of 2 between t₀ and t_TAMS, or the two components have merged. To separate the hardening events caused by cluster dynamics (i.e. by the assistance of other cluster member(s)) from those due to stellar evolution, we take into account only the wider binaries (with log₁₀(P[days]) ≳ 3 at t₀), which provides a lower estimate on the total number of hardening events. We choose the minimum orbital period for hard binaries (log₁₀(P[days]) ≳ 3) to be longer than the minimum period for a Cepheid with a companion on a circular orbit (log₁₀(P[days]) ≈ 2.6) to exclude the majority of the cases where the orbital period shrunk because of tidal circularisation during stellar binary evolution; only ≈15% of the BiCeps with log₁₀(P[days]) = 3 have eccentricity high enough (e ≳ 0.67; Kobulnicky et al. 2014) to circularise to orbital period of log₁₀(P[days]) ≲ 2.6.

We note that the adopted definition of soft and hard binaries differs from the classical definition of soft and hard binaries (Heggie 1975; Hills 1975), where the change of the orbital period is purely due to dynamical encounters, while the present simulations include also the orbital period change due to stellar evolution (which is caused by mass loss, mass transfer or the supernova kick). We adopt this definition because the standard models do not allow us to directly distinguish the influence of dynamical encounters and stellar evolution, but a comparison with control models provides an estimate of the influence of stellar evolution only.

Figure 5 contrasts the outcome of binary evolution of the given initial orbital period distribution in star clusters (upper row) with the control models (lower row) for the clusters at the galactocentric radius R_g = 8 kpc and metallicity Z = 0.014. More than 97% of ProCeps are identified as being a binary member at t₀ (the remaining 3% are immediately disrupted by the presence of other stars in the cluster), and all the binaries are shown in the histograms. The outcomes in dependence of the orbital period distribution are separated by log₁₀(P[days]) ≈ 3, where practically all binaries with log₁₀(P[days]) < 3 interact with their companion (see Sec. 4.3 below), while cluster dynamics plays a more significant role for the BiCeps with log₁₀(P[days]) > 3, which interact with their companions infrequently.

As expected, both hardening (magenta regions) and softening events (cyan regions) are more pronounced in models with cluster dynamics (p_hard = 0.05 and p_soft = 0.18; Table 1) with respect to the control models (p_hard = 0.02 and p_soft = 0.12). The orbital period P_sep for the watershed semi-major axis a_hard Eq. (3) which separates soft binaries from hard binaries, ranges from log₁₀(P_sep[days]) ≈ 3.2 to log10(P_sep[days]) ≈ 6.2 as M_ecl decreases from 6400 M_⊙ to 50 M_⊙ (grey bands in Fig. 5). The clusters behave in a more complex way than expected from this simple picture because there are softening events for P < P_sep, and hardening events for P > P_sep, but on average softening events increase with P, hardening events decrease with P, and hardenings are more common in less massive clusters, which have a larger value of P_sep. We note that some of the softenings and hardenings are due to stellar evolution only (lower row). In this case, the most of the softenings is produced by stellar mass loss of the primary, and the hardenings are produced when the former primary receives a supernova kick in a suitable direction.

Fig. 3 compares the distribution of orbital periods of BiCeps in more massive clusters (cyan line), low mass clusters (magenta), and in the field (black dashed line). BiCeps in more massive clusters (M_ecl ≳ 1000 M_⊙) have shorter orbital periods (50% of BiCeps have log₁₀(P[days]) < 4.1) than BiCeps in less massive clusters (50% of BiCeps have log₁₀(P[days]) < 5.9 in clusters with M_ecl ≲ 500 M_⊙), which indicates that the dynamical environment of more massive clusters preferentially destroys Cepheids with longer orbital periods. The orbital period distribution of BiCeps in the field lies between these two distributions, but closer to the distribution of lower mass clusters for longer orbital periods because lower mass clusters contribute more Cepheids to the field because of their higher evaporation rate and earlier disintegration. The median orbital period of BiCeps, Q₅₀ = med(log₁₀(P[days])), and its lower (Q₂₅) and upper (Q₇₅) quartile in different environments are listed in Table 2. According to the Table, BiCeps in clusters tend to have longer orbital periods (med(log₁₀(P[days])) = 5.3) than those in the field (med(log₁₀(P[days])) = 4.5).

Companions to Cepheids with orbital periods shorter than ≈10 years are of particular interest because they can be inferred from long-term variations of radial velocity of Cepheids. Shetye et al. (2024) find that 15 ± 2% of Cepheids have a companion with log₁₀(P[days]) < 3.6, which is higher than in our simulations (7% for the set of models with Z = 0.014 and R_g = 8 kpc)⁴.

The orbital period distribution for Cepheids in clusters and the field combined (yellow line in Fig. 3) has ≈25% of binaries with P larger than the initial upper limit of the period distribution, which is log₁₀(P[days]) = 6.7. These BiCeps are not primordial, but they form dynamically during the cluster evolution; they are absent in the control models (brown line). A closer inspection shows that the BiCeps with log₁₀(P[days]) ≳ 6.7 have different origins in lower mass and more massive clusters: lower mass clusters form the wide BiCeps mainly during complicated interactions between the cluster stars and binaries throughout the cluster life-time, often including several consecutive exchanges; while more massive clusters form wide BiCeps mainly during gas expulsion at the outskirts of the clusters, with many of the BiCeps escaping the cluster in the process. The latter mechanism was studied in a general setting by Kouwenhoven et al. (2007, 2010). The dichotomy in the formation mechanism probably accounts for the higher ratio of p_BC,out/p_BC,in in more massive clusters (left panel of Fig. 4). Most of the long period binaries are disrupted when the Cepheid becomes a white dwarf or neutron star, with the compact object receiving a supernova kick.

Fig. 4 Fraction of Cepheids in binaries and triples as a function of the initial stellar mass of the embedded cluster where the Cepheids formed. Left panel: fraction of Cepheids in binaries of orbital period shorter than 1010 days. The fraction for Cepheids occurring in clusters and in the field is plotted separately by solid and dashed lines, respectively. Clusters at different galactocentric radii are shown by symbols, and clusters of different metallicity are shown by colours. Some of the Cepheids in binaries are orbited by an outer body, namely, they are the inner binary of a triple. Right panel: fraction of Cepheids in triples where the orbital period of the outer body is shorter than 1010 days. The meaning of the symbols is the same as in the left panel.

Fig. 5 Outcome of binary evolution in ProCeps as a function of their initial orbital period. Upper row: the panels from left to right show the future of ProCeps formed in a population of star clusters of the assumed ECMF of Eq. (4) for clusters of mass Mecl < 800 M⊙, clusters with Mecl ≳ 800 M⊙, and all clusters combined. BiCeps which harden, soften or exchange their companions are shown by the magenta, cyan and yellow regions, respectively. Cepheids which have coalesced with their companion are shown by green regions, and ExCeps are shown by red regions. These evolutionary outcomes are shown as stacked histograms. The initial orbital period distribution for all Cepheids is shown by the thick black histogram (the hatched regions represent the rest of BiCeps, which have not underwent any of the aforementioned evolutionary changes). We note that the colour-coded histograms are not additive; for instance an exchange that hardens at the same time is shown in both areas as if it were two separate binaries. The orbital period of a binary can change between t0 and the time when the star becomes a Cepheid (see e.g. Fig. 6). The vertical black dashed lines (at log10(P[days]) = 2.6) show the upper orbital period for physical collision with their companion during the Cepheid phase. The vertical grey bands represent the threshold orbital period Psep between the soft and hard binaries for clusters of the particular mass range. The plots are for a population of star clusters with Rg = 8 kpc and Z = 0.014. Lower row: same as the upper row but for the control models with the same metallicity.

Fig. 6 Evolution of the orbital period and eccentricity between the ZAMS (grey triangles) and TAMS (coloured bars) for the binaries which contain at least one ProCep. Although the plot shows all the binaries at the ZAMS, it shows the state at the TAMS only of the binaries which originate with (log10(P), e) from the interval (3, 7) × (0, 0.25) (the green rectangle). The lines connect the states at the ZAMS with those at the TAMS. The purpose of the dot and line colour is only to aid clarity. The black lines indicate the maximum allowed eccentricity of a binary of given orbital period to avoid a physical collision with its companion for m1 + m2 = 9 M⊙ and stellar radii of R1 = 5 R⊙ (on the MS) and R1 = 60 R⊙ (before the second passage through the instability strip). Left panel: control models. Only the minority of binaries evolve; these are the secondaries where the primary mass loss results in an increase of the orbital period. Right panel: models with cluster environment, which result in a pronounced increase of eccentricity. Both cases are calculated for Z = 0.014, Rg = 8 kpc and for the cluster mass range from Mecl = 50 to Mecl = 6400 M⊙. Note: the binaries which either form or are destroyed between the ZAMS and TAMS are not shown.

4.1.3 Dynamical exchanges

Although some Cepheids are in binaries, their companions are not always the same stars with which the Cepheids formed. This happened due to a dynamical exchange involving at least three stars. Alternatively, some of these BiCeps originate from the primordial binaries which capture an outer body forming a temporal triple, and where the primordial (inner) binary merges later. The outer body becomes a new companion to the merged object, with the possibility of the Kozai-Lidov mechanism triggering some of the coalescences (Perets & Fabrycky 2009; Naoz & Fabrycky 2014); in both cases, the identity of the companion has changed. The fraction p_chg of BiCeps which have exchanged their companion between time t₀ and TAMS is ≈0.10 (Table 1). This quantity increases with decreasing cluster mass from ≈0.05 for the 6400 M_⊙ clusters to ≈0.20 for the 50 M_⊙ clusters. Exchanges occur in BiCeps over the whole orbital period range with a mild preference for longer orbital periods (the yellow area in the upper row of panels in Fig. 5).

4.1.4 Cepheids with compact companions

Next, we study the fraction of Cepheids p_cmp which have a compact companion. This estimate is subject to large uncertainties because it depends not only on the details of orbital period changes during possible Roche-lobe mass transfer, but also on the supernova kick velocity distribution, which is not well-known (e.g. Hansen & Phinney 1997; Hobbs et al. 2005; Fragos et al. 2009; Beniamini & Piran 2016; Giacobbo & Mapelli 2020). Assuming that the supernova kicks follow a Maxwell distribution with dispersion σ = 40 km s⁻¹, ≈3 to 5% of all Cepheids have a compact companion (Table 1). The compact companion is a white dwarf in ≈80% of the cases, a neutron star in ≈16% of the cases and a black hole in ≈4% of the cases. The orbital period distribution of Cepheids with compact companions is more narrow (interquartile range of Q₇₅ − Q₂₅ = 1.4; Table 2; the dotted line in Fig. 3) than for Cepheids with non-degenerate companions with practically no binaries of log₁₀(P[days]) ≳ 7.

As expected, BiCeps with compact companions have a smaller mass ratio than the BiCeps with gaseous companions; the median and interquartile range is log₁₀(m₂/m₁) of −0.67 and 0.07 for the former, and −0.32 and 0.54 for the latter. The mass ratio is significantly more narrow for compact companions because of the relatively small mass range allowed for neutron stars and white dwarfs originating from stars of mass ≳4 M_⊙.

Fig. 7 Change in orbital period and eccentricity between t0 and tTAMS (upper row), and between t0 and the end of the core helium burning tHe,end (lower row). The colour-scale shows the ratio between the number of ProCeps (at tTAMS) or Cepheids (at tHe,end) and the number of ProCeps at t0, which are located at the nearby position in the log10(P) − e plane. Values larger (lower) than 1 indicate that the number of ProCeps or Cepheids increased (decreased) since t0. The dashed and dotted white lines show the period threshold for coalescence (at log10(P[days]) = 2.6) and the upper limit on the birth period distribution (log10(P[days]) = 6.7), respectively. The black line shows the maximum allowed eccentricity for given orbital period at t0 (assuming m1 + m2 = 9 M⊙ and stellar radius of R1 = 5 R⊙). The absence of binaries with log10(P[days]) ≲ 2.6 at tHe.end is due to stellar collisions. The depletion of binaries with log10(P[days]) ≳ 3 and e ≲ 0.5 (the white rectangle with an indication of eccentricity evolution) seen already at the tTAMS in the standard model is due to dynamical encounters with other cluster members (the depletion is absent in the control model). The binaries with log10(P[days]) ≳ 6.7 in the standard model (and absent in the control model) are mostly outer orbits in triples.

4.1.5 Cepheids originating from secondaries

In the present model of binary stars at metallicity Z = 0.014, approximately 26% of ProCeps formed as secondaries (the remaining 74% of ProCeps formed as primaries). The ProCeps which formed as secondaries are more impacted (by collisions or mass transfer) with their more massive and faster evolving companions than the ProCeps which formed as primaries: Only 38% of the ProCeps which formed as secondaries become Cepheids, while 80% of the ProCeps which formed as primaries become Cepheids. Among all RanCeps, only 14% were initially secondary ProCeps (the rest were primaries). The RanCeps which originated from secondaries have a substantially smaller binary fraction (p_BC ≈ 0.15) than the RanCeps which originated from primaries (p_BC ≈ 0.5). The lower binary fraction for the secondary Cepheids is likely caused by the supernova kick to the primaries and the dynamical disruption of the binary as described in the paragraph below. The majority of the binary RanCeps which were initially secondaries have a compact companion: ≈30% of their companions are neutron stars, ≈50% white dwarfs, and only ≈20% gaseous stars, which are typically also Cepheids or red giants with a slightly larger mass than the secondary.

Apart from stellar evolution in binaries, the Cepheids which formed as secondaries are more influenced by stellar dynamics, usually in the form of the binary-binary or similar interaction between stellar subsystems, which gives the former secondaries enough speed to escape their birth clusters and possibly releases them from their companions. The ejected stars travel at larger speeds and thus they can be found at larger distances from the clusters. We take the cut-off distance to be 300 pc because it is the distance travelled by a star at a speed corresponding to the velocity dispersion of a 500 M_⊙ cluster (2 km s⁻¹) for the 150 Myr life-span of Cepheids at the Solar metallicity. While 34% of the Cepheids which were initially secondaries are located farther away than 300 pc from their birth clusters, this fraction is only 11% for the Cepheids which were initially primaries. This behaviour is expected from N-body dynamics where binary-binary or binary-single star interactions preferentially release the least massive stars from the subsystem, i.e. the secondaries (Valtonen 1974; Heggie & Hut 2003; Tanikawa et al. 2012).

Fig. 8 Eccentricity as a function of period for BiCeps in standard models (magenta) and control models (black). The thick lines show the mean eccentricity for given period, and the coloured areas represent the interquartile range. The figure quantitatively shows the boost of eccentricity for Cepheids which are subjected to cluster dynamics. The mean eccentricity at time t0 (identical for both models) is shown by the yellow line.

4.1.6 Binaries where both the primary and secondary are Cepheids

Some Cepheids have as a companion another Cepheid. To calculate their fraction p_C–C among all Cepheids, we divide the time interval of (0, 300 Myr) to bins of 30 kyr, which are sufficiently short to resolve the movement of the Cepheid in the instability strip for the most massive Cepheids (they have the fastest evolution), and we count the Cepheids, N_C,i, in each bin i as well as the Cepheids, N_C–C,i, in a binary with another Cepheid. The definition p_C–C = ∑_i N_C–C,i/ ∑_i N_C,i, where i runs over all time bins, provides the expected observed fraction of Cepheids residing in a binary of two Cepheids when averaged over Cepheids of all possible ages forming within a galaxy of a time-independent star formation rate. Both standard and control models provide p_C–C = 0.013 with median orbital periods log₁₀(P[days]) ≈ 4.6, which is slightly shorter than the orbital period distribution of all Cepheids (Table 2). The small stock of statistics available does not allow for the metallicity dependence of this quantity to be properly studied. We note that Karczmarek et al. (2022) found lower values of p_C–C ≈ 0.001 to 0.01 for most of their binary models.

4.1.7 Evolution of eccentricity

Star cluster dynamics is responsible for substantial changes in eccentricity. The lines in Fig. 6 connect the states in the orbital period-eccentricity plane (log₁₀(P) − e) of the binaries containing at least one ProCep at the ZAMS with their states at the TAMS. The comparison is done at the TAMS because the time-span covers the majority of the stellar life-time (≈90%) so that cluster dynamics has caused most of its evolution of orbital parameters while stellar evolution has still relatively minor influence. For the sake of clarity, the Figure shows the evolution only of the binaries having 3 < log₁₀(P[days]) < 7 and e < 0.25 at the ZAMS.

In the control models (left panel), only a small fraction of binaries undergo evolution of their orbital parameters between the ZAMS and TAMS. Most of the binaries that evolve their orbital parameters increase their orbital periods at constant eccentricity, which is due to the stellar mass loss from the primary in the binaries where the Cepheid is the secondary. In rare cases, some binaries increase their eccentricity, which is due to the SN explosion and kick in the primaries where the Cepheid is the secondary.

In contrast, cluster dynamics causes more significant changes of orbital parameters (right panel of Figure 6), particularly an increase of eccentricity. For example, in the control models practically all BiCeps with e > 0.75 had already e > 0.75 at the ZAMS (this is 11% of all BiCeps), while in the standard models cluster dynamics increases the fraction of BiCeps with e > 0.75 to 23%. Moreover, cluster dynamics also dissolves about half of the initially highly eccentric binaries so that only approximately one in four of the highly eccentric BiCeps had e > 0.75 already at the ZAMS (6% of all future BiCeps have e > 0.75 at the ZAMS), the rest being produced by encounters in the cluster.

The ratio between the number of ProCeps in binaries at t_TAMS and t₀ at a given area in the log₁₀(P) − e plane is shown in Fig. 7⁵. Between these two time events, the interaction of the binaries with other cluster members reduces the population of low eccentric binaries (e < 0.25) with 3 < log₁₀(P[days]) < 7 to 27% (upper right panel), which is a substantially larger depletion than due to stellar evolution only (population reduced to 70%; upper left panel). Cluster environment also forms binaries of longer orbital periods (log₁₀(P[days]) > 7) than present at t₀.

When the ProCeps become Cepheids, cluster environment leaves its imprints in the eccentricity of Cepheids with log₁₀(P[days]) ≳ 3.5 (Fig. 8). For example, Cepheids with log₁₀(P[days]) ≈ 6 have median eccentricity of 0.38 ± 0.25 (the error represents the interquartile range) in the control models, while they have median eccentricity of 0.73 ± 0.16 in the standard models. Future observations of eccentricity distribution among Cepheids with wide companions might help constrain the environment in which they were born.

Interestingly, albeit the depopulation of the low eccentricity binaries with 3 < log₁₀(P[days]) < 7 is striking and it practically leads to a desert (upper right panel of Fig. 7), this fact is not so apparent from the eccentricity distribution of all ProCeps as indicated by the dotted histogram in the left panel of Fig. 9. This is because the population of low eccentric binaries is dominated by the binaries of short orbital periods (log₁₀(P[days]) ≲ 2.5), whose dynamical evolution is less significant because of their compactness and thus smaller cross sections for interaction (Hut & Bahcall 1983). As the ProCeps evolve beyond the TAMS, their radii increase, and they interact with their companions. In this process, practically all ProCeps with log₁₀(P[days]) ≲ 2.5 interact with their companion so that they are either prevented from becoming Cepheids, or they are single stars when becoming Cepheids; in any case the states with log₁₀(P[days]) ≲ 2.5 are erased as the stars evolve from the TAMS to t_He,end (lower row of Fig. 7).

Since the close binaries are predominantly the ones which have low eccentricity at the TAMS, the eccentricity distribution changes between the t_tams and t_He,end for the standard models substantially. While the eccentricity distribution becomes almost flat in the standard models, it is only marginally modified in the control models between t₀ and the Cepheid stage (cf. the dashed lines in the left panel of Fig. 9), implying that binary evolution only is not able to change the eccentricity distribution to a larger extent. Thus, the desert for log₁₀(P[days]) ≳ 3 and e ≲ 0.5 and the flattening of eccentricity distribution is the result of cluster dynamics.

Fig. 9 The eccentricity distribution (left panel) and mass ratio distribution (right panel) of ProCeps and Cepheids in binaries. The solid, dotted and dashed line shows the distribution at t0, tTAMS, and tHe,end, respectively. The plot is for Cepheids formed in a population of star clusters with Rg = 8 kpc, Z = 0.014, and of the adopted ECMF of slope β = 2 Eq. (4).

4.1.8 Evolution of the mass ratio

Cluster evolution does not noticeably influence the mass ratio distribution between the BiCep and its companion (right panel of Fig. 9). The median mass ratio log₁₀(m₂/m₁) decreases from −0.16 for ProCeps at t₀ to −0.32 for Cepheids in the models including cluster dynamics (for models with R_g = 8 kpc, Z = 0.014) and to −0.29 in control models (for models with Z = 0.014).

4.2 Cepheids in triples and quadruples

We adopt the definition of the fraction of Cepheids in triples p_tr as the ratio of all Cepheids with at least two companions to the total number of Cepheids. As with binaries, we include only the systems for which the outer body has the orbital period around the inner binary shorter than 10¹⁰ days. The fraction of Cepheids in triples as a function of the initial mass of their birth cluster M_ecl is shown in the right panel of Fig. 4 separately for Cepheids in clusters (p_tr,in; solid lines) and in the field (p_tr,out; dashed lines). The plots provide a lower limit of Cepheids in triples because present models do not include primordial triples; all triples form dynamically. The triple fraction is larger in clusters (p_tr,in = 0.31) than in the field (p_tr,out = 0.11), and p_tr,{in,out} decreases with increasing M_ecl with only a slight dependence on R_g or Z. This is similar to the trend seen in binaries (the binary fraction is 1.7× higher in cluster than in the field), albeit it is more pronounced in triples (the triple fraction is 3× higher in cluster than in the field). The triple fraction of Cepheids in clusters and field combined is p_tr ≈ 0.18.

Outer orbits (i.e. the orbit of the outer star around the inner binary) have nearly thermal eccentricity distribution (left panel of Fig. 10, where f(e) ∝ e; Heggie 1975), which is consistent with their formation due to capture. This distribution stands in contrast to the flat eccentricity distribution of binary Cepheids (left panel of Fig. 9). Another difference between the distribution of binaries and outer orbits in triples is in the mass ratio distribution, which is more flat and bimodal with peaks at log₁₀(m₃/m₁) ≈ −0.9 and ≈ − 0.1 (right panel of Fig. 10), while the distribution of binaries increases with log₁₀(m₂/m₁) (right panel of Fig. 9). Triples in clusters tend to have outer companions of more equal masses than triples in the field (c.f. the cyan and black lines in the right panel of Fig. 10). The mass ratio distribution for triples is calculated as the ratio of the outer body of mass m₃ relatively to the Cepheid, which is of mass m₁. In less frequent cases, where a single Cepheid is being orbited by a binary, m₃ is taken to be the mass of the binary.

The period distribution of outer orbits is shown in Fig. 11 separately for Cepheids in different environments: for all Cepheids (yellow line), all Cepheids in field (black line), Cepheids in more massive clusters (M_ecl > 1000 M_⊙; cyan line) and Cepheids in lower mass clusters (M_ecl < 500 M_⊙; magenta line). In contrast to BiCeps (Fig. 3), the distribution does not depend sensitively on the environment. Outer orbits also have much longer orbital periods with median log₁₀(P[days]) of 8.5 and interquartile range of Q₇₅ − Q₂₅ = 2.2 (Table 2), nevertheless some outer orbits have orbital periods comparable to that of BiCeps (≈20% of outer orbits have log₁₀(P[days]) ≲ 6.7). Figures 10 and 11 are drawn for Cepheids formed by a population of star clusters with R_g = 8 kpc, Z = 0.014, and of the adopted ECMF of slope β = 2 (Eq. (4)).

The fractions of Cepheids in quadruples in clusters p_quad,in and in the field p_quad,out follow similar trends as binaries and triples: lower mass clusters have a significantly larger fraction of quadruples (p_quad,in ≈ 0.10 and p_quad,out ≈ 0.05) than more massive clusters, which have p_quad,in, p_quad,out ≲ 0.01. The combined quadruple fraction of Cepheids in clusters and in the field is p_quad ≈ 0.04. Given the low number statistics of quadruples in our simulations, we do not detect the dependence of p_quad,in and p_quad,out on R_g or Z. For a particular cluster mass, the absolute number of quadruples is lower by a factor of 5 than the number of triples.

All triples in our simulations form dynamically. As we describe in sect. 4.2 of Paper I, in lower mass clusters ProCeps mass segregate to the cluster centre before the cluster density is diluted as the result of gas expulsion. Their elevated density and relatively low velocity near the cluster centre enables binary formation where one or both of the inner components of the binary is another binary, so a triple or quadruple forms dynamically (Binney & Tremaine 2008; Moeckel & Bate 2010). This is the same mechanism we discuss in the context of binary formation at the end of Sect. 4.1.2.

An example of dynamical formation of a stellar multiplet that includes several ProCeps is shown in the online material⁶. The movie shows a 100 M_⊙ cluster (model M7), which mass segregates its ProCeps before gas expulsion, so that they are retained near the cluster centre, where they survive the potential rebound due to gas expulsion (see Sect. 4.2 of Paper I for details). Located in a small volume at the cluster centre, all the ProCeps in the clusters rapidly (by 1.1 Myr) form a hierarchical sextuplet. The sextuplet ejects a binary containing one of the Cepheids at ≈3.6 Myr. During the ejection, the new quadruple (consisting of two hard binaries orbiting each other) captures three low mass stars, forming a hierarchical septuplet. All the captured stars have eccentric orbits. The three captured stars can be seen in the movie (upper right panel) after the ejection event.

The formation of hierarchies is preferred in lower mass clusters because they contain a small number of stars of mass similar to ProCeps and these stars mass segregate rapidly to the cluster centre, which shortens the timescale for multiple star formation (Binney & Tremaine 2008, their Sect. 7.1). The multiple star system does not need to be inherently stable, but its decay-time (which, for unstable systems, corresponds to ≈10–100 crossing times; Heggie & Hut 2003; Wang et al. 2019) is relatively long (even comparable to t_He,end in some cases) so that the Cepheid might be still in a multiple stellar system. This explains why the fraction of Cepheids in triples increases with decreasing cluster mass. Triples also need time to be formed dynamically, which explains why field Cepheids (most of which escape when the clusters are young) have a lower triple fraction than the ones in clusters. Another way how to transport triples to the field is to eject them; however, ejecting a triple out of the cluster is not easy as these systems can be disrupted in the encounter which would eject it.

Fig. 10 Eccentricity distribution (left panel) and mass ratio distribution (right panel) of the outer bodies in triples including Cepheids. The distribution for all triples (i.e. in star clusters and in the field combined) and the triples in the field and clusters only is shown by the yellow, black and cyan lines, respectively. We note that both distributions are remarkably different from the distributions of Cepheids in binaries, which are shown in Fig. 9.

Fig. 11 Cumulative distribution of orbital periods of the outer orbits in triples including Cepheids. We show separately the distribution for all Cepheids in triples in clusters and field combined (yellow lines), only in field (black lines), the Cepheids originating from more massive clusters (cyan), and the Cepheids originating from lower mass clusters (magenta). The distribution is normalised to the total number of triples N3,tot in each group.

Fig. 12 Cumulative distribution of the time of coalescence for the Ran-Ceps that experienced a coalescence normalised to the ttams of each individual Cepheid. Results for clusters with Rg = 8 kpc are plotted separately for the metallicity Z = 0.014 (red line), Z = 0.006 (green) and Z = 0.002 (blue). The approximate time when the stars are located within the Cepheid instability strip is indicated by the grey area.

4.3 Stellar mergers and mass transfer events

We define the number fraction of ExCeps (i.e. stars which are prevented from becoming Cepheids by a merger, or Roche-lobe mass transfer, or common envelope evolution with their companion) as p_EC ≡ 1 − N_RC/N_PC, where N_RC is the total number of RanCeps and N_PC is the total number of ProCeps at time 0. The value of p_EC is independent of the Galactic radius R_g (p_EC ≈ 0.30 for Z = 0.014; Table 1) with a possible weak decrease with decreasing metallicity (p_EC ≈ 0.27 for Z = 0.002). These fractions are in good agreement with the estimate p_EC = 0.25 ± 0.15 provided by Moe & Di Stefano (2017, their Sect. 6.1). Cluster dynamics increases p_EC by approximately 20% (cf. with the control models in Table 1). The vast majority of ExCeps ultimately merge with their companion; only ≲4% of ExCeps are caused solely by mass transfer without a merger event.

While some mergers with companions prevent ProCeps from becoming Cepheids, other mergers still produce Cepheids. The fraction of Cepheids which have experienced a merger (i.e. the total number of Cepheids which have merged divided by the total number of Cepheids) and which are located in clusters (p_C,coll,in) and in the field (p_C,coll,out) are listed separately in Table 1. The Cepheids located in clusters have merged more often than the Cepheids in the field (e.g. p_C,coll,in = 0.45 and p_C,coll,out = 0.36 for R_g = 8 kpc and Z = 0.014). The difference is probably caused by the hardening events, which are possible only within star clusters. We do not detect a significant dependence of p_C,coll,in and p_C,coll,out on R_g or Z, except for a possible slight increase for Z = 0.002.

At which evolutionary stage is the ProCep most likely to merge with its companion? Fig. 12 shows the cumulative histogram of all Cepheids which have merged with their companion before time t normalised to t_tams for each Cepheid. Approximately 12% of mergers occur still on the MS for Z = 0.014; the ratio decreases in lower metallicity environments, which is probably because of the lower mass threshold for Cepheids (i.e. they have smaller radius) and generally smaller stellar radius at a lower metallicity for given stellar mass on the MS (Tout et al. 1996), both decreasing the closest distance between stars for an interaction. The majority of mergers (60–70%) occur shortly after t_tams when the star substantially increases its radius. Approximately 20% of mergers occur after the Cepheid stage (i.e. after t_He,end), which has no influence on the Cepheid.

Some ProCeps in clusters have merged with wide companions (log₁₀(P[days]) > 4; upper row of Fig. 5) even though such cases are excluded from purely binary evolution because even moderately eccentric (e ≈ 0.5) massive primaries (m ≳ 15 M_⊙) can interact with their companions during their red supergiant phase only when log₁₀(P[days]) ≲ 3.8 (Bertelli et al. 2009; Moe & Di Stefano 2017). These cases are absent in our control models (lower row of Fig. 5). We suspect that these events originate from the Kozai-Lidov mechanism (Perets & Fabrycky 2009; Naoz & Fabrycky 2014), where the initial relatively mild eccentricity of the inner binary was boosted by the outer component in a triple. The outer component was captured during the cluster evolution before the star became Cepheid. This idea is supported by the fact that ≈93% of the Cepheids which have log₁₀(P[days]) > 4 at t₀ and which have merged with their companions have an outer companion at t_He,end (this companion might have caused the merger via the Kozai-Lidov mechanism), while the percentage of the Cepheids having log₁₀(P[days]) < 4 at t₀ and a distant companion is only 42%.

The Cepheids which have experienced stellar mergers or Roche-lobe mass transfer in their history appear in their host clusters at a different time (usually later) than what is expected for their mass (pulsation period) and the cluster age, which makes them peculiar objects. Fig. 13 shows the cumulative histogram of the ratio between the expected age t_exp(m) of the Cepheid calculated from its mass to its physical age t_phys (taken as the time since t₀). The age of a substantial number of Cepheids (≈50%) has been influenced by their interaction with their companion, with some of them (≈3%) having t_exp(m) shorter by more than a factor of 3 than their physical age. This opens the interesting possibility for the star clusters which host multiple Cepheids that one or several of the Cepheids could appear to be of different age than the other Cepheids within the same cluster even though they actually formed in a coeval star forming event. These Cepheids, whose age can be estimated from period-age relations (Anderson et al. 2016b; De Somma et al. 2021), would not match the age of the cluster estimated from isochrone fitting using the same set of stellar evolutionary models. Therefore, Cepheids in clusters might be imprecise indicators of the cluster age, often underestimating it. A merger which appears younger than the rest of the cluster and which is likely to evolve to a Cepheid in the future has already been observed in the Pleiades (star Alcyone of mass ≈6.5 M_⊙ Brandt & Huang 2015).

A small fraction (≈3%) of Cepheids appear to be older by more than 10% than expected from their mass (see Fig. 13), i.e. these stars are of too low mass for their age. Most of these cases had been eccentric binaries before their evolution out of the MS, whereupon they interacted with their companion in such a way that resulted in some mass loss and tidal circularisation; nevertheless the secondary still provided the primary with enough room to become a Cepheid. These binaries usually coalesce immediately after the Cepheid stage. Apart from this evolutionary path, some of the Cepheids with t_exp(m)/t_phys > 1.1 originate from mergers of stars with very high angular momentum, which lose more mass in the merger than what they gain from the secondary because of angular momentum conservation.

In total, from 15 to 30% of all Cepheids found in clusters have either coalesced with their secondary or gained significant amount of mass from their former primary so that their derived age is substantially different (by more than 40%) from the age of the cluster. The fraction of these cases is labelled p_mma in Table 1, and they are represented by the grey regions of Fig. 13. We point out that Roche-lobe mass transfer and stellar mergers, which are crucial in determining the outcome of binary interaction, and the subsequent evolution of the merger are treated by NBODY6 using a simplified treatment of the complex hydrodynamics of stellar mergers, and this may influence individual results. Nevertheless, the main conclusions should reflect the real behaviour. A more detailed discussion of the adopted simplifications and their possible influence can be found in Sec. 6.1.

Next, we studied the Cepheids which originate from stars of the initial mass below the lower mass threshold m_min,Ceph of Cepheids (AdCeps of type I), but which later increase their mass above the Cepheid mass threshold by a merger or mass transfer event with their companion. These stars are searched for in the data in the same way as RanCeps; an example of such an evolutionary track is depicted in the right panel of fig. A.1 in Paper I. Our models with cluster dynamics terminate at 300 Myr, but by this time, they form a comparable number of AdCeps-I as the control models. It is probable that the main runs would produce a similar number of AdCeps-I also after 300 Myr, so we take the value of the relative number of AdCeps-I, p_ACI ≡ N_ACI(t′ < 2Gyr)/N_PC(0) formed in the control runs as a proxy for what is expected to be found observationally. The number N_ACI(t′ < t) of AdCeps-I formed in the control runs by time t is shown by the dotted lines in Fig. 14. AdCeps-I start appearing slightly before the time, t_He,end, and 80% of all AdCeps-I occur before 270 Myr, 340 Myr, and 390 Myr for Z = 0.014, Z = 0.006, and Z = 0.002, respectively. A small number of AdCeps-I occur as late as 1 Gyr. In contrast, RanCeps occur earlier with no RanCeps occurring after t_He,end (solid lines). The total fractions of AdCeps-I are p_ACI = 0.17, 0.11, and 0.14 for Z = 0.014, Z = 0.006, and Z = 0.002, respectively (Table 1). This means that a stellar population producing N_RC Cepheids from ProCeps also produces p_ACIN_RC/(1 − p_EC) Cepheids originating from stars of mass below m_min,Ceph. Thus, according to the adopted binary evolution recipes, approximately one in five observed Cepheids originates from a star with initial mass below the lower mass threshold m_min,Ceph assumed for Cepheids, and it becomes a Cepheid only due to the mass gain from its companion in a binary.

A small fraction of Cepheids ≲0.3% in our models originate from stars of mass above m_max,Ceph (AdCeps type II). These objects typically lost some mass to the secondary, but were not prevented by the secondary from becoming Cepheids. However, we caution that such rare and complex objects would require a specific in-depth study before firmer conclusions can be drawn.

Fig. 13 Cumulative distribution of the ratio of the Cepheid age texp(m) as expected from its present mass to its physical age tphys (i.e. the time since its formation) for all Cepheids which are located in star clusters. Cluster metallicity is shown by colour. The grey regions indicate the cases where texp(m) differs by more than 40% from tphys. These cases, whose fraction among all Cepheids is denoted pmma, are likely to be detected by observers as Cepheids of discrepant ages.

Fig. 14 Cumulative number of stars which have become Cepheids by the time, t, for the control simulations normalized to the number of ProCeps which is taken at the beginning. Cepheids originating from ProCeps (solid lines) are plotted separately from the ones originating from stars of m < mmin,Ceph (dotted lines). The metallicity of the models is indicated by colour. The thin vertical solid lines show the time tHe,end, namely, the time of the least massive non-interacting Cepheid which can occur at a given metallicity; all the Cepheids occurring later than tHe,end gained some of their mass in interaction with their companion.

4.4 Role of the initial orbital period distribution

We estimated the uncertainty in the initial period distribution by comparing the control model of Z = 0.014, which has the orbital period distribution of Sana Eq. (1), with the control model for the same metallicity but which populates more the longer orbital periods Eq. (2). The latter model likely presents an upper limit on the width of the orbital period distribution for ProCeps (cf. Fig. 1).

ProCeps with the wider distribution of orbital periods Eq. (2) experience fewer stellar mergers with p_C,coll = 0.27 (p_C,coll = 0.40 for the control models with Sana’s orbital period distribution of Eq. (1); hereafter in this paragraph, we list the results of the original Sana distribution in parentheses), and the fraction of ExCeps is also lower p_EC = 0.17 (p_EC = 0.24). The less extensive merger history then decreases the fraction of Cepheids of mismatched age to p_mma = 0.14 (p_mma = 0.20), and it increases the fraction of binary Cepheids to p_BC = 0.59 (p_BC = 0.50). We do not detect much difference in p_cmp nor in p_C–C. The results are summarised in Table 1 (last row). In total, even though the models with the wider orbital period distribution were taken as the likely upper limit on how wide the distribution can be, they do not change the results dramatically; the wider distribution changes any of the quantity listed in Table 1 by at most 50%, but it is generally less for most of the quantities.

5 On the birth companion frequency among B stars

5.1 Comparison of present models to observations

The estimates of the fraction of Cepheids in triples and quadruples are likely to be lower bounds because our simulations start with subsystems of a maximum of two bodies (there are no primordial triples or higher order hierarchies in our models), while many young stars of early spectral types are found in more numerous subsystems, for instance, as in the Trapezium in the Orion Nebula Cluster (Preibisch et al. 1999). The dynamical encounters in clusters mainly disrupt some of the binaries, but they can also form other binaries, triples and higher order subsystems. We recall that our models produce Cepheids with a binary fraction of 42% for cluster and field Cepheids combined. However, based on Gaia DR2 data, Kervella et al. (2019b) provides a lower estimate on Cepheids in binaries as high as 80%, which implies that assuming even 100% of primordial binarity among ProCeps is insufficient to account for the observed binarity of Cepheids. Possible solutions are primordial triples or higher hierarchies of ProCeps or perhaps a different initial orbital period distribution of ProCeps with most ProCeps having orbital periods between 10² to 10⁴ days (avoiding very short orbits to prevent mergers which produce single Cepheids). Such an orbital distribution would be in tension with observations of Sana et al. (2012), suggesting that the more likely explanation is mid-B star formation in hierarchical systems.

The fraction of Cepheids in triples is also lower in our models (≈18% for Cepheids in clusters and field combined) than what is observed (≈35% according to Evans et al. 2005 and Kervella et al. 2019b). When comparing the companion frequency of Cepheids, which we calculate as f_comp = p_BC + p_tr + p_quad⁷, we obtain f_comp = 0.64 for the present models, which is almost two times lower than the observed value of ≈1.15 (Evans et al. 2005; Kervella et al. 2019b). The observed value is still a lower estimate because it misses some binaries and triples, and it neglects quadruples.

Recently, Evans et al. (2020) find that all the Cepheids with a companion closer than ≈2000 AU (log₁₀(P[days]) ≈ 7.1) in their sample have another inner companion, so the systems are hierarchical triples. This differs from the occurrence of close companions to all Cepheids, where only ≈30% of all Cepheids are in spectroscopic binaries (Evans et al. 2015; Shetye et al. 2024); namely, the likelihood of a Cepheid to have a close companion is dependent on the presence of another wider companion. Clearly, the origin of these triples can be twofold: either primordial or dynamically formed. We searched our models for similar triples, and found that only ≈5% of Cepheids with a companion of orbital period shorter than log₁₀(P[days]) = 7.1 have another inner companion. Thus, it is unlikely that these outer companions formed dynamically; they are the result of the star formation process.

The necessity of primordial triples or even higher order stable hierarchical systems to explain the observed companion frequency of Cepheids points in the same direction as earlier works, which are available only for late-type stars (Kroupa 1995b; Portegies Zwart et al. 2004; Kouwenhoven et al. 2007; Moeckel & Bate 2010; Geller et al. 2013), where the observed companion frequency cannot be matched by simulations assuming primordial binaries only even if the birth binary fraction is as high as 100%. Our findings support these results and extends them to more massive primaries, where the necessity of primordial hierarchies is important not only for the fraction of triples and more numerous systems in evolved clusters, but also for the fraction of binaries.

5.2 Estimate of the true birth companion frequency of B stars

The cluster dynamics and binary evolution can be viewed as mathematical operators transforming the birth companion frequency of stars around the spectral type mid-B to the companion frequency of Cepheids⁸. Cepheids enable detection of less massive and more distant companions than MS mid-B stars because of their narrower spectral lines (Evans et al. 2013, 2015). In our models, the evolutionary operator transforms the initial mid-B star binary fraction of 1.0 (i.e. companion frequency of 1) to the Cepheid binary fraction of 0.42. Given the observed binary fraction for Cepheids to be ≳0.8, and assuming that the operator is linear (which it is likely not), we arrive at an estimate for the birth companion frequency of mid-B stars of ≳2.0, namely, all mid-B stars form at least in stable triples on average (with some proportion of quadruples compensated by binaries). This estimate is substantially higher than the observed companion frequency for mid-B MS stars (1.3 ± 0.2) (Moe & Di Stefano 2017); see also Abt et al. (1990).

Setting aside the companions to Cepheids which elude observation, this is still likely a lower estimate. More precisely, the evolutionary transformation is not only from mid-B stars ↦ Cepheids, but from mid-B stars and some lower mass stars ↦ Cepheids because some of the lower mass binaries gain mass by a merger with their companion and become AdCeps-I (Sec. 4.3). The lower mass stars are of late-B and early-A type. Approximately 20% of all Cepheids originate from lower mass stars (AdCeps-I), whereas the remaining 80% of Cepheids originate from the expected mass range (RanCeps, Sec. 4.3). At the ZAMS, the companion frequency of these stars is then given by the weighted mean of both populations, $$f_{\text{comp }}=w_1{f}_{\text{comp },1}+w_2{f}_{\text{comp},2},$$(6)

where w₁ and f_comp,1 are, respectively, the number fraction and birth companion frequency for lower mass stars, and w₂ and f_comp,2 are, respectively, the number fraction and birth companion frequency for ProCeps. Following the results of Sec. 4.3, we set w₁ = 0.2 and w₂ = 0.8. The median mass at the ZAMS is 4 M_⊙, and 6.2 M_⊙ for the future AdCeps-I and RanCeps, respectively. From Offner et al. (2023, the right panel of their Fig. 1; see also Shatsky & Tokovinin 2002; Rizzuto et al. 2013; Moe & Di Stefano 2017; Moe & Kratter 2021 for the observational basis), we extrapolate the observed companion frequencies to be 1.29 and 1.55 for the 4M_⊙, and 6.2 M_⊙ star, respectively. Assuming that the true birth ratio of f_comp,2/f_comp,1 is the same as the observed one (i.e. ≈1.2), Eq. (6) provides f_comp,2 = 2.05 and f_comp,1 = 1.7 for mid-B and late-B stars, respectively. This would indicate a substantially higher companion frequency for late-B stars than currently assumed; a typical late-B star would form as a triple. On the other hand, assuming that the true birth companion frequency for late-B stars is the same as the observed one (i.e. f_comp,1 = 1.29), Eq. (6) provides f_comp,2 = 2.2, which is even more than the observational estimate for O stars (f_comp = 2.1; Sana et al. 2012). In any case, the correction with respect to AdCeps-I increases the estimated companion frequency for mid-B stars above 2.

Another piece of evidence that the birth companion frequency of mid-B stars is substantially larger than 1.3 stems from the comparison with the observed companion frequency of Cepheids, which is at least 1.15 (Sec. 5.1). Given that ≈25% of mid-B primaries do not evolve to Cepheids because of the interaction with their companion (Moe & Di Stefano 2017), the mid-B star companion frequency would decrease from 1.3 to 1.0 due to stellar evolution alone, i.e. below the value of 1.15 observed for Cepheids, which is itself a lower estimate.

The estimate of the birth companion frequency for mid-B stars must be taken with great caution. The main limitation is the assumption that the evolutionary operator is linear. Including stable hierarchical triples would open more possibilities for interactions between these systems, different interaction cross-sections, and also the long-term stability of perturbed hierarchical systems. New sets of simulations starting with primordial triples will have to be calculated to constrain the birth companion frequency for B stars more accurately.

The additional uncertainty is the onset of gas expulsion, t_d, which we take as 0.6 Myr. If gas expulsion starts later than this, clusters evolve for longer in their compact sizes, which means that soft binaries will be exposed to the dense cluster environment for longer, and thus more of them will be ionised. Thus, if t_d > 0.6 Myr (t_d < 0.6 Myr), a higher (lower) birth companion frequency will be needed to evolve to the same companion frequency of Cepheids than what was obtained for t_d = 0.6 Myr. The properties of multiple systems at birth that we can deduce from the observations thus intimately depend on the feedback physics and compact state of their birth embedded clusters as emphasised by Kroupa (1995a,b).

In the estimate of the birth companion frequency, we assume that all triples are intrinsically stable. This assumption does not consider primordial unstable triples (Goodwin & Kroupa 2005), which typically decay to a hard binary and a single star. Including unstable triples (and possibly also unstable hierarchies of higher order), would raise the mid-B star birth companion frequency above 2.

6 Discussion

6.1 Stellar evolution of Cepheids in N-body models

In this paper, we have relied on a simplified picture of the evolution of stars into the evolutionary phase where they may become observable as classical Cepheids. The broad strokes of this picture are well understood. However, several open questions remain, notably concerning the mass-luminosity (M-L) relation and its dependence on mixing processes, mass-loss, among others. N-body simulations integrate stellar evolution via fast lookup functions (Aarseth 1999; Hurley et al. 2000) since full stellar evolution calculations, not to mention detailed modelling of binary evolution, would be prohibitively expensive. The fast lookup functions (Hurley et al. 2002) in NBODY6 depend on a set of detailed 1D stellar model calculations and notably incorporate tides that ensures realistic orbital evolution of close-in binaries on their way to the Cepheid stage. However, the detailed modelling of stellar mergers is a very complex subject with many uncertainties and our results, especially for rare cases or the properties of merger products, should be interpreted very carefully. For details on binary evolution and the uncertainties involved in modelling the common envelope phase, we refer to Hurley et al. (2002). Since predictions will generally differ between stellar evolution codes, a couple of warnings are in order.

Firstly, stellar ages are always model-dependent, whereas the dynamical timescale of cluster evolution only indirectly depends on the stellar evolutionary timescale (e.g. mass transfer triggered by evolution of the primary). Since the ages of Cepheids depend strongly on input physics (Bono et al. 2005; Anderson et al. 2016a; De Somma et al. 2021), care must be taken when interpreting the results presented here in combination with the output of stellar evolution models, such as period-age relations for Cepheids. For example, stellar evolution models featuring rotation predict Cepheids to occur at ages that are approximately 60–100% older than models that do not include the rejuvenating effects of stellar rotation on the MS (Anderson et al. 2016a). Since such changes are relative within an age scale set for a given set of models, this does not affect the statement above concerning apparent ages due to mergers or mass transfer events. However, the dynamical time affecting cluster evolution and the stellar evolution time scale causing evolution-driven interactions could be systematically at odds.

Secondly, the mass range within which Cepheids occur depends on stellar evolution models. We here have considered Cepheids within the mass range where the stellar evolution incorporated within NBODY6 predicts blue loops to intersect with a theoretically predicted instability strip (Anderson et al. 2016a) that has been shown to agree closely with observed Cepheid populations (Groenewegen 2020; Espinoza-Arancibia et al. 2024). This ensures internal consistency in our treatment of Cepheids. However, this leads to more massive and younger (mixing effects aside, cf. immediately above) Cepheids than those predicted by other models, notably ones shown not to be subject to the long-standing Cepheid mass discrepancy (e.g. Caputo et al. 2005; Anderson et al. 2014). In particular, the minimum mass of Cepheids is correlated with the extent of the blue loops, which are particularly uncertain stellar evolutionary effect (e.g. Walmswell et al. 2015). However, differences in stellar model masses are of order 10 to 20% and thus clearly subdominant given the factor of ≈3 in mass covered by Cepheids (≈4 to 12 M_⊙). While caution is warranted near the extremes of the mass range, most trends reported here are expected to be rather insensitive to the problem of mass scale.

6.2 Initial binary population of B stars

As we mention in Sect. 2, we assume that the initial orbital period distribution of B stars is close to the observed one. This assumption is adopted mainly because of the lack of previous theoretical works on this topic. Sec. 4.1 shows that the properties of binaries evolve substantially due to cluster dynamics, with many longer orbital period binaries (log₁₀(P[days]) ≳ 4) being softened, which usually leads to disruption (upper row of Fig. 5). If the initial period distribution were more populated at longer orbital periods (log₁₀(P[days]) ∈ (≈5, ≈8); analogously to the shape of the initial period distribution for late type stars as found by Kroupa 1995b) than observation by Sana et al. (2012), most of the widest binaries would be disrupted early. This would produce more single Cepheids, and thus decrease the fraction of Cepheids in binaries, the fraction of Cepheids which have merged or those which are of different age than expected from their current mass.

On the other hand, the estimated initial companion frequency of 2 (Sec. 5.2) implies very different initial conditions than estimated here. We can view the present paper as the first iteration towards understanding the true birth companion frequency of B stars. The second iteration is to perform N-body models of star clusters with mid-B stars initially in triples and quadruples, and then to compare the outcome with available observations, and possibly adjust the initial conditions again to perform another iteration until convergence is reached.

The choice of initial conditions for multiple subsystems is more challenging than for binary stars because of the possible correlations between the orbital parameters of the inner and outer bodies (see, e.g. White & Ghez 2001 for circumstellar disc correlation in early type pre-MS stars). The presence of initial triples or quadruples likely increases the probability of collisions between stars (e.g. by the Kozai-Lidov mechanism) and interaction between the subsystems because the subsystems would mass segregate to the cluster core early, where they would exchange some of their members, while ejecting the others as single stars or hard binaries (Perets & Šubr 2012).

The lack of primordial triples in our simulations implies that our results regarding the fractions of Cepheids in binaries, triples and quadruples underestimate their real values in absolute numbers. Despite this, the dependence of these quantities on cluster mass and cluster or field environment (as dealt with in Sects. 4.1 and 4.2) is instructive and likely follows the real trends.

The present work neglects Cepheids during their first crossing through the instability strip because of the short duration of the event. Although the more numerous low mass stars (which do not have second and third instability strip crossing) are expected to cross the instability strip once, and spend there a longer time than the stars that evolve from ProCeps, observational evidence suggests that these objects are rare. Among more than 3600 Milky Way classical Cepheids (Pietrukowicz et al. 2021), there are only four whose rapid period changes indicate a first instability strip crossing (Turner et al. 2006; Fadeyev 2014; Anderson et al. 2016b) and the pulsation periods of these stars do not imply a particularly low mass. Since identifying period changes requires very long temporal baselines (and since first crossings are more easily identified due to faster period changes), there is only limited room for a large population of low-mass stars that would only enter the instability strip during a first crossing. That said, period changes have only been measured for a couple hundred Milky Way Cepheids (Csörnyei et al. 2022), so it is not excluded that individual such cases could be found once sufficiently long temporal baselines become available for all known Milky Way Cepheids. A similar paucity of Cepheids on their first crossing is observed in the Large Magellanic Cloud (Rodríguez-Segovia et al. 2022), which has a lower metallicity (Z = 0.006) than Milky Way (Z = 0.014). The small population of observed Cepheids on the first crossing implies that they cannot noticeably influence the companion frequency of Cepheids and the estimate on the birth multiplicity of B stars.

6.3 The dependence on galactocentric radius and metallicity

The dissolution rate of star clusters is a strong function of the cluster orbital radius, R_g (e.g. Baumgardt & Makino 2003). It means that the cluster density and the time interval for which binaries are exposed to cluster environment depend also on R_g. However, we do not detect any clear trend on galactic radius among the properties of binaries, triples or mergers studied in this work (see Table 1). This result is in contrast to the radial dependence of the fraction of all Cepheids in clusters, which increases by a factor of 1.2 as R_g increases from 4 kpc to 12 kpc (Paper I). The probable explanation for the independence of binary properties of R_g is that most of the binary dynamical evolution (i.e. due to stellar encounters in the cluster environment) occurs during the earliest stage of cluster evolution (before gas expulsion at t_d ≈ 0.6 Myr or shortly afterwards), which is not influenced by the external field of the galaxy. The explanation is in line with numerical works (e.g. Marks & Kroupa 2012) and observations (e.g. Bouvier et al. 1997; Patience et al. 1998; Kraus et al. 2011; Geller & Mathieu 2012; Moe & Di Stefano 2017), which indicate that the statistical properties of binary stars in clusters stay largely constant from the age of several Myr to the age of several Gyr.

There is a very weak dependence of the binary properties on metallicity (Table 1). As the metallicity decreases, Cepheids reach larger radii (Anderson et al. 2016b), which slightly increases the merger rate with their companions (higher p_{C,coll,{in,out}} in Table 1). The lower minimum Cepheid mass, m_min,Ceph, at lower Z implies that a higher fraction of Cepheids originate as the secondary, while the primary had become a compact object when the secondary becomes a Cepheid. Thus, the fraction of Cepheids with a compact companion p_cmp increases with decreasing metallicity (Table 1). We note that the fraction of Cepheids with mismatched age, p_mma, decreases with decreasing Z substantially in control models according to the Table, but this is artificial because the simulations terminate at 300 Myr, while the mergers that evolve into Cepheids live longer for these metallicities. The control models, which terminate at 3 Gyr show no trend of p_mma with metallicity (Table 1).

7 Summary

We study the properties of stellar companions to Cepheids by means of direct N-body simulations. These models start when all stars (including mid-B stars that evolve into Cepheids) form in open star clusters as binaries (but not more numerous subsystems) on the ZAMS. The binaries are influenced not only by binary star evolution (which is modelled according to the recipes of Hurley et al. 2002), but also by changes of their orbital parameters due to interactions with other stars in the cluster, including dynamical exchanges, binary hardenings, and ionisations. We calculate the properties of binaries and multiple systems of these stars at the time when they evolve to Cepheids, and discuss the implications for the birth companion frequency of mid-B stars. The influence of the dynamical cluster environment is explored here for the first time; the only other population synthesis of Cepheids, which is due to Neilson et al. (2015) and Karczmarek et al. (2022) consider only stellar evolution of isolated binaries. The main results of the modelling are summarised below.

• During binary evolution, many ProCeps interact with their companions by coalescence, Roche-lobe mass transfer or common envelope evolution. Approximately 30% of ProCeps are prevented from becoming Cepheids due to binary evolution. In contrast, mass gain in lower mass stars enable some of the stars to become Cepheids even though they would not have been able to do so had they evolved in isolation. Approximately one in five Cepheids actually originates from stars with mass below the threshold mass for isolated single Cepheids. Cepheids from these stars appear at a later time (some of them even at the age of ≳500 Myr) than expected for Cepheids originating from ProCeps.

• A significant fraction of Cepheids (≈45% of Cepheids in clusters and ≈36% of Cepheids in the field) have merged with their companion; however, this does not prevent them from becoming Cepheids. Most of the mergers (60–70%) occur immediately after the TAMS (Fig. 12). The larger fraction of mergers in clusters is attributed to dynamical hardening. As a result of mergers or mass transfer in binaries, 15–30% of all Cepheids in clusters occur at a different time (more different than by 40%) than expected from their stellar mass and the cluster age (Fig. 13).

• The percentage of Cepheids in a binary is substantially higher in clusters (≈60%) than in the field (≈35%). Cepheids also have different orbital period distribution in different environments as the result of cluster dynamics (Sect. 4.1.2 and Fig. 5), with ≈10% of all binary Cepheids having exchanged their companion since the time of their formation (Sect. 4.1.3).

• Approximately 3–5% of all Cepheids have a compact companion (Sect. 4.1.4). These binaries have a lower mass ratio (q ≈ 0.2) than all Cepheids in binaries (q ≈ 0.5).

• Binaries where both the primary and secondary are Cepheids at the same time are found in ≈1.3% of all Cepheids (Sect. 4.1.6).

• Dynamical evolution in clusters substantially increases the eccentricities of binaries which have initially log₁₀(P[days])≳3 and e ≲ 0.5, forming a desert in the corresponding area of the log₁₀(P) − e plane (upper right panel of Fig. 7). This transforms the eccentricity distribution from decreasing with e (at the ZAMS) to being flat (when the stars become Cepheids; left panel of Fig. 9). The mass ratio distribution between the secondary and primary is not significantly influenced by cluster dynamics.

• The percentage of Cepheids in dynamically formed triples is larger in the less massive clusters than in the more massive ones (right panel of Fig. 4). In total, 30% of Cepheids in clusters and 10% of Cepheids in the field are in triples (Sec. 4.2). A similar dependence is seen in quadruples with 5 and 3% incidence among Cepheids in clusters and in the field, respectively. Some of the triples can trigger coalescence of the inner binary by increasing its eccentricity via the Kozai-Lidov mechanism. A movie showing an example of dynamical formation and disruption of a multiple system of four Cepheids is attached in the accompanying online material.

• In addition to stellar evolution, cluster dynamics decreases the fraction of all Cepheids (i.e. in clusters and field combined) in binaries (p_BC = 0.42 instead of p_BC = 0.50 in control models; see Table 1), and it decreases the fraction of Cepheids with compact companions (p_cmp = 0.03 instead of p_cmp = 0.05). On the other hand, cluster dynamics causes more Cepheids of mismatched age (p_mma = 0.27 instead of p_mma = 0.20), and prevents more ProCeps from becoming Cepheids (p_EC = 0.30 instead of p_EC = 0.24).

• Binary properties are unaffected by the position of the cluster within its galaxy (Sec. 6.3). In lower metallicity environment, there are more Cepheids with compact companions and slightly more Cepheids which have merged with their companion, but the overall impact of metallicity is minor (Sec. 6.3).

• Comparing the companion frequency in present models (≈0.64) with the observed value (≳ 1.15), we suggest that mid-B stars form with companion frequency substantially larger than 1, and probably larger than 2 (Sec. 5.2). This implies that mid-B stars typically form in triples and even more numerous subsystems, most of which are stable on a timescale of at least ≈1 Myr.

• The estimated true companion frequency of mid-B stars thus significantly differs from the value assumed in our work (i.e. f_comp = 1). Starting our simulations with f_comp = 2 would likely result in different properties of binary, triple, and multiple Cepheids than reported here, with mergers, Roche-lobe mass transfer events, hardenings, exchanges, and ionizations being even more abundant. However, the dependence of binary properties on the cluster mass and field or cluster environment likely reflects real trends. The present simulations likely provide a lower estimate on these quantities.

Data availability

Movie is available at https://www.aanda.org and at https://zenodo.org/records/13889050

Acknowledgements

We thank the referee for the constructive and very detailed comments, which improved the quality and clarity of the paper. F.D. appreciates the ESO Visiting Scientist Programme, which made this project possible via his stay at ESO Garching. We would like to thank to Sverre Aarseth for continuously developing the code NBODY6, and to Sambaran Banerjee and Seungkyung Oh for their advice about specifying integration parameters for the simulations. F.D. thanks to Stefanie Walch for her encouragement to work on the project while being employed in her group, and acknowledges support of the Deutsche Forschungsgemainschaft (DFG) through SFB 956 “The conditions and impact of star formation” (sub-project C5). R.I.A. is funded by the SNSF Eccellenza Professorial Fellowship PCEFP2_194638 and acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement No. 947660). F.D. and P.K. acknowledge support through grant 20-21855S from the Grant Agency of the Czech Republic. P.K. acknowledges support through the DAAD-Easter-Europe Exchange grant at Bonn University and corresponding support from Charles University. We appreciate the support of the ESO IT team, which was vital for performing the presented numerical models. This research made use of the Matplotlib Python Package (Hunter 2007).

Appendix A Comparison of terminology with Paper I

References

All Tables

All Figures

Fig. 1 Cumulative distribution functions of orbital periods of binary stars according to the primary mass. The observed distributions for primaries in specific mass intervals are shown by bars (Moe & Di Stefano 2017): O-type stars (m > 16 M⊙; blue bars), mid B-type stars (5 M⊙ < m < 9 M⊙; green), late B- and A-type stars (2 M⊙ < m < 5 M⊙; yellow), and solar type stars (0.8 M⊙ < m < 1.2 M⊙; red). The theoretical birth distribution of Kroupa (1995a) for late type stars (m < 3 M⊙; dashed black line) is transformed by pre-MS eigenevolution to the initial distribution (solid black line). The initial distribution for O-stars according to Eq. (1) (Sana et al. 2012) and the wider distribution of Eq. (2) with A = 3, which is the widest estimate on the Cepheid orbital periods, is plotted by the magenta and cyan line, respectively.

In the text

Fig. 2 Schematic illustration of the considered mass changes due to binary evolution. The abscissa represents time; the ordinate represents the initial stellar mass, m. ProCeps, which have an initial mass range of (mmin,Ceph, mmax,Ceph) evolve towards Cepheids (without mass exchange: violet line or with mass exchange: grey lines), the state of which is marked by the hatched area. The beginning of core helium burning is indicated by the yellow dotted line. See Sect. 3 for full explanation of the evolutionary channels and the terminology used in this paper.

In the text

Fig. 3 Cumulative distribution of orbital periods of ProCeps (taken at t0; black solid line) and BiCeps (taken at tHe,end; dashed lines). The orbital period distribution of BiCeps is shown separately for Cepheids in the lower (magenta dashed line) and more massive clusters (cyan dashed line), in the whole cluster population (yellow dashed line), in the field originating from the whole cluster population (black dashed line), and in the control model (brown dashed line). The orbital period distribution of Cepheids with a compact companion is indicated by the black dotted line. The distributions are normalised to the total number of BiCeps in each group. There are almost no BiCeps with log10(P[days]) ≲ 2 because of coalescence with their companion before the onset of core helium burning. The vertical dotted line represents the upper limit on the initial orbital periods (log10(P[days]) ≲ 6.7). All BiCeps with an orbital period longer than that formed due to dynamical encounters in clusters.

In the text

Fig. 4 Fraction of Cepheids in binaries and triples as a function of the initial stellar mass of the embedded cluster where the Cepheids formed. Left panel: fraction of Cepheids in binaries of orbital period shorter than 1010 days. The fraction for Cepheids occurring in clusters and in the field is plotted separately by solid and dashed lines, respectively. Clusters at different galactocentric radii are shown by symbols, and clusters of different metallicity are shown by colours. Some of the Cepheids in binaries are orbited by an outer body, namely, they are the inner binary of a triple. Right panel: fraction of Cepheids in triples where the orbital period of the outer body is shorter than 1010 days. The meaning of the symbols is the same as in the left panel.

In the text

Fig. 5 Outcome of binary evolution in ProCeps as a function of their initial orbital period. Upper row: the panels from left to right show the future of ProCeps formed in a population of star clusters of the assumed ECMF of Eq. (4) for clusters of mass Mecl < 800 M⊙, clusters with Mecl ≳ 800 M⊙, and all clusters combined. BiCeps which harden, soften or exchange their companions are shown by the magenta, cyan and yellow regions, respectively. Cepheids which have coalesced with their companion are shown by green regions, and ExCeps are shown by red regions. These evolutionary outcomes are shown as stacked histograms. The initial orbital period distribution for all Cepheids is shown by the thick black histogram (the hatched regions represent the rest of BiCeps, which have not underwent any of the aforementioned evolutionary changes). We note that the colour-coded histograms are not additive; for instance an exchange that hardens at the same time is shown in both areas as if it were two separate binaries. The orbital period of a binary can change between t0 and the time when the star becomes a Cepheid (see e.g. Fig. 6). The vertical black dashed lines (at log10(P[days]) = 2.6) show the upper orbital period for physical collision with their companion during the Cepheid phase. The vertical grey bands represent the threshold orbital period Psep between the soft and hard binaries for clusters of the particular mass range. The plots are for a population of star clusters with Rg = 8 kpc and Z = 0.014. Lower row: same as the upper row but for the control models with the same metallicity.

In the text

Fig. 6 Evolution of the orbital period and eccentricity between the ZAMS (grey triangles) and TAMS (coloured bars) for the binaries which contain at least one ProCep. Although the plot shows all the binaries at the ZAMS, it shows the state at the TAMS only of the binaries which originate with (log10(P), e) from the interval (3, 7) × (0, 0.25) (the green rectangle). The lines connect the states at the ZAMS with those at the TAMS. The purpose of the dot and line colour is only to aid clarity. The black lines indicate the maximum allowed eccentricity of a binary of given orbital period to avoid a physical collision with its companion for m1 + m2 = 9 M⊙ and stellar radii of R1 = 5 R⊙ (on the MS) and R1 = 60 R⊙ (before the second passage through the instability strip). Left panel: control models. Only the minority of binaries evolve; these are the secondaries where the primary mass loss results in an increase of the orbital period. Right panel: models with cluster environment, which result in a pronounced increase of eccentricity. Both cases are calculated for Z = 0.014, Rg = 8 kpc and for the cluster mass range from Mecl = 50 to Mecl = 6400 M⊙. Note: the binaries which either form or are destroyed between the ZAMS and TAMS are not shown.

In the text

Fig. 7 Change in orbital period and eccentricity between t0 and tTAMS (upper row), and between t0 and the end of the core helium burning tHe,end (lower row). The colour-scale shows the ratio between the number of ProCeps (at tTAMS) or Cepheids (at tHe,end) and the number of ProCeps at t0, which are located at the nearby position in the log10(P) − e plane. Values larger (lower) than 1 indicate that the number of ProCeps or Cepheids increased (decreased) since t0. The dashed and dotted white lines show the period threshold for coalescence (at log10(P[days]) = 2.6) and the upper limit on the birth period distribution (log10(P[days]) = 6.7), respectively. The black line shows the maximum allowed eccentricity for given orbital period at t0 (assuming m1 + m2 = 9 M⊙ and stellar radius of R1 = 5 R⊙). The absence of binaries with log10(P[days]) ≲ 2.6 at tHe.end is due to stellar collisions. The depletion of binaries with log10(P[days]) ≳ 3 and e ≲ 0.5 (the white rectangle with an indication of eccentricity evolution) seen already at the tTAMS in the standard model is due to dynamical encounters with other cluster members (the depletion is absent in the control model). The binaries with log10(P[days]) ≳ 6.7 in the standard model (and absent in the control model) are mostly outer orbits in triples.

In the text

Fig. 8 Eccentricity as a function of period for BiCeps in standard models (magenta) and control models (black). The thick lines show the mean eccentricity for given period, and the coloured areas represent the interquartile range. The figure quantitatively shows the boost of eccentricity for Cepheids which are subjected to cluster dynamics. The mean eccentricity at time t0 (identical for both models) is shown by the yellow line.

In the text

	Fig. 9 The eccentricity distribution (left panel) and mass ratio distribution (right panel) of ProCeps and Cepheids in binaries. The solid, dotted and dashed line shows the distribution at t0, tTAMS, and tHe,end, respectively. The plot is for Cepheids formed in a population of star clusters with Rg = 8 kpc, Z = 0.014, and of the adopted ECMF of slope β = 2 Eq. (4).
In the text

Fig. 10 Eccentricity distribution (left panel) and mass ratio distribution (right panel) of the outer bodies in triples including Cepheids. The distribution for all triples (i.e. in star clusters and in the field combined) and the triples in the field and clusters only is shown by the yellow, black and cyan lines, respectively. We note that both distributions are remarkably different from the distributions of Cepheids in binaries, which are shown in Fig. 9.

In the text

Fig. 11 Cumulative distribution of orbital periods of the outer orbits in triples including Cepheids. We show separately the distribution for all Cepheids in triples in clusters and field combined (yellow lines), only in field (black lines), the Cepheids originating from more massive clusters (cyan), and the Cepheids originating from lower mass clusters (magenta). The distribution is normalised to the total number of triples N3,tot in each group.

In the text

Fig. 12 Cumulative distribution of the time of coalescence for the Ran-Ceps that experienced a coalescence normalised to the ttams of each individual Cepheid. Results for clusters with Rg = 8 kpc are plotted separately for the metallicity Z = 0.014 (red line), Z = 0.006 (green) and Z = 0.002 (blue). The approximate time when the stars are located within the Cepheid instability strip is indicated by the grey area.

In the text

Fig. 13 Cumulative distribution of the ratio of the Cepheid age texp(m) as expected from its present mass to its physical age tphys (i.e. the time since its formation) for all Cepheids which are located in star clusters. Cluster metallicity is shown by colour. The grey regions indicate the cases where texp(m) differs by more than 40% from tphys. These cases, whose fraction among all Cepheids is denoted pmma, are likely to be detected by observers as Cepheids of discrepant ages.

In the text

Fig. 14 Cumulative number of stars which have become Cepheids by the time, t, for the control simulations normalized to the number of ProCeps which is taken at the beginning. Cepheids originating from ProCeps (solid lines) are plotted separately from the ones originating from stars of m < mmin,Ceph (dotted lines). The metallicity of the models is indicated by colour. The thin vertical solid lines show the time tHe,end, namely, the time of the least massive non-interacting Cepheid which can occur at a given metallicity; all the Cepheids occurring later than tHe,end gained some of their mass in interaction with their companion.

In the text