© The Authors 2025

1. Introduction

The formation and evolution of globular clusters (GCs) have long been subjects of interest for astrophysics and cosmology (Krumholz et al. 2019), with numerous questions that are bound to be resolved thanks to recent observational advances. A major breakthrough came over the past two decades from Hubble Space Telescope (HST) photometry combined with ground-based high-resolution spectroscopy on very large telescopes (e.g., Bedin et al. 2004; Mucciarelli et al. 2009; Piotto et al. 2015; Dalessandro et al. 2016; Marino et al. 2019; Gilligan et al. 2019; Rain et al. 2019; Carretta & Bragaglia 2024). This synergy revealed that all old Galactic and extragalactic GC scrutinized so far in the local universe host multiple stellar populations (MPs) with unique chemical properties. While their so-called first population (1P) stars share the chemical properties of halo stars of similar metallicity, their second population (2P) stars exhibit diverse degrees of N, Na, Al, and He enhancements that are anticorrelated with O (and sometimes Mg) depletion (see Gratton et al. 2019 for a review and Carretta & Bragaglia 2023 for more recent insight and references). The proportion of 2P stars in individual GC ranges between 30% and 80% and increases with the GC mass (Carretta 2006; Carretta et al. 2010; Milone et al. 2017; Lagioia et al. 2019; Carretta & Bragaglia 2024). Similarly, the extent of the O-Na and Mg-Al anticorrelations increases with GC mass, with Mg depletion being observed only in the most metal-poor and/or most massive GC (e.g., Bragaglia et al. 2017; Pancino et al. 2017; Mucciarelli et al. 2018). Helium variations between 1P and 2P are usually very small (ΔY≃0.01), with the strongest ones (ΔY≃0.1) derived in the most massive GC (e.g., Milone et al. 2018; Cadelano et al. 2023).

The new photometric window opened by the James Webb Space Telescope (JWST) on the coolest low-mass main sequence stars in GC has provided evidence that the same MP characteristics are found all over the color-magnitude diagram (Nardiello et al. 2022; Milone et al. 2023; Ziliotto et al. 2023; Scalco et al. 2024). Similar MP signatures have also been detected in extra-galactic massive star clusters with ages as low as ∼1.5 Gyr (e.g., Mucciarelli et al. 2008; Saracino et al. 2020; Martocchia et al. 2020; Kapse et al. 2022; Salgado et al. 2022).

There is a general agreement on the fact that 2P stars formed out of proto-GC gas mixed to various degrees (so-called dilution) with hydrogen-burning yields processed at high temperature and released in important quantities (so-called mass-budget problem) by a short-lived type of 1P so-called polluters (Prantzos et al. 2017, and references therein). Some debate persists however regarding the identity of the polluters (e.g., Charbonnel 2016; Bastian & Lardo 2018; Milone & Marino 2022, and references therein). However, we know that 2P stars inherited only H-burning products, with very slight amount of He and no products of more advanced nucleosynthesis, and they are seen only in massive and compact star clusters, young and old, in addition to the mass budget issue (e.g., Prantzos & Charbonnel 2006; Schaerer & Charbonnel 2011). This indicates that we are dealing with a very specific type of stellar polluters, which formed concurrently with their host proto-GC in very specific conditions and at any metallicity.

We followed the approach from Denissenkov & Hartwick (2014), who demonstrated that the dilution of the central H-burning products of a 10⁴ M_⊙ super massive star (SMS) close to the zero age main sequence (ZAMS) with original proto-GC can reproduce the ubiquitous O-Na anticorrelation and the scarce Mg-Al anticorrelation. Gieles et al. (2018) proposed a scenario where SMS with masses above 10⁴ M_⊙ can form at any redshift through runaway collisions in gas-rich and extremely compact cluster environments, independently of metallicity¹. The continuous rejuvenation of the SMS through collisions and gas accretion in Gieles et al. (2018) model can provide enough H-burning ejecta to fulfil the requisite mass budget to account for the ratios between 1P and 2P stars.

However, the required stellar density to form SMS in their runaway model is about four orders of magnitude higher than that inferred in the proto-GCs that have been identified in young and confirmed bound systems in gravitationally lensed galaxies at redshifts up to 10.6–11 (Adamo et al. 2024). In addition, Gieles et al. (2018) assumed total mass conservation regarding the products of the collisions. This is the hypothesis we are aiming to test in this paper.

On the other hand, Prantzos et al. (2017) showed that protocluster pollution by early main sequence stars with masses between ∼10³ and 10⁴ M_⊙ (which they also call SMS, and with actual values depending on metallicity) can also explain the abundance properties of MPs. According to the inertial inflow (Padoan et al. 2020) GC formation model by Gieles et al. (2025), stars in this mass range could form in proto-GCs (with masses of ≳10⁵ M_⊙) via gas accretion as a natural extension of the IMF beyond very massive stars (VMS, with masses between ∼10² and 10³ M_⊙), whose existence has been established in extragalactic super star clusters (e.g., Crowther et al. 2016; Smith et al. 2016, 2023; Meštrić et al. 2023; Shenar et al. 2023; Upadhyaya et al. 2024). Here, we follow the nomenclature of Gieles et al. (2025), who called such objects accreting extremely massive stars (aEMS; 10³⁻⁴ M_⊙). In their model, aEMS can reach an equilibrium between gas accretion and wind mass loss. This conveyor belt behavior is expected to rejuvenate the aEMS by supplying fresh H, making the star stay in the MS phase as long as the conveyor belt phase, thereby allowing it to live longer than what would traditionally be expected. Consequently, it releases H-burning yields with low He contamination and in sufficient quantities regarding the mass budget issue. However, stellar collisions, which cannot be avoided in dense stellar systems, might also interfere during the growth of VMS and aEMS, potentially leading to an intermediate picture between the runaway SMS and the aEMS scenarios to explain the origin of MPs in GCs.

The interest in collisions in star clusters also follows a long history. A number of theoretical and numerical studies have demonstrated that under specific conditions, stellar collisions and mergers can potentially result in the formation of VMSs, SMSs, and eventually IMBHs (Sanders 1970; Quinlan & Shapiro 1989; Portegies Zwart & McMillan 2002; Gürkan et al. 2004; Freitag et al. 2006; Giersz et al. 2015; Mapelli 2016; González Prieto et al. 2024; Vergara et al. 2025). However, the collision path and its end product have been demonstrated to be influenced by several factors (e.g., Glebbeek et al. 2009; González Prieto et al. 2024). The initial conditions of the cluster (mainly its mass and radius) represent the primary parameter that might drive its early and long-term dynamical evolution. These can be constrained by the type of observations described above or derived from galactic simulations in a cosmological context with the requisite resolution to allow the formation of bound clusters to be followed in a robust manner (e.g., Li & Gnedin 2020; Rodriguez et al. 2023; Bruel et al. 2024; De Lucia et al. 2024; Pascale et al. 2025). Other parameters are important, such as the stellar density and initial mass function (IMF), mass segregation, and the initial fraction of binaries.

The model predictions are also affected by the treatment of the evolution of single and binary stars (mass-radius relation, lifetime,...) as well as the assumptions about the resulting merger objects (mass, radius, rejuvenation, remnant), which are crucial, yet highly uncertain ingredients. Smoothed particle hydrodynamical simulations have been employed to model the collision of two stars, initially focusing on low-mass stars in the context of exploring the formation of blue stragglers in star clusters (e.g., Lombardi et al. 2002), but also extending to more massive stars (e.g., Gaburov et al. 2008), which are the potential progenitors of neutron stars and black holes. Further evolving a merger product is numerically challenging, as it involves mapping back the still not fully relaxed stellar mass distribution from the three-dimensional (3D) hydrodynamic simulation to a one-dimensional (1D) stellar evolution code (e.g., Sills et al. 1997, 2001; Glebbeek et al. 2013; Schneider et al. 2016, 2019, 2020; Wu et al. 2020). All the aforementioned studies look at stellar collisions with mass ratio in the range of ≃1−10 and do not consider stellar masses in the VMS+EMS+SMS regime. Therefore, the hydro and N-body simulations of massive and compact star clusters lack necessary and tractable input for mass-loss induced by collisions. Furthermore, while some of these models follow the dynamical evolution of dense star clusters over extended periods, they do not address the early phases of the evolution of their constituent stars before they arrive on the ZAMS. Given that pre-main sequence stars have larger radii and different internal structures and chemical compositions than their main sequence descendants, this could have implications for the collision rates, the collision-induced mass loss, and the merger product, as well as the maximum mass a star could reach in a cluster.

The present study aims to provide predictions for collision-induced mass loss using an analytical approach based on the mass-dependent structure of aEMSs that allowed us to cover a broad range of stellar masses and radii for both the aEMS targets and the colliding stars in the very early phases of cluster formation. In Sect. 2, we describe how we built the aEMS models that are required for our analysis; we discuss how the uncertainties in the transport of energy in their convective layers affect the predictions for their global properties and compare with models from the literature. The nucleosynthesis behavior of aEMSs is presented in Appendix A. In Sect. 3, we set the analytical prescriptions we then used to quantify collision-induced mass loss and gain. We identified the main processes at work, considering inspiraling and target stars over a broad mass range. We provide grids with model predictions that can be used for N-body and/or hydro simulations of star clusters. As a first example of application, in Sect. 4, we detail our use of a Monte Carlo approach to simulate multiple collisions on an aEMS target in a dense static proto-cluster during its first 5 Myrs and we derived the amount of mass it can lose and gain considering collisions, gas accretion, and mass loss via stellar winds. We summarize our results and the sources of uncertainties of our predictions in Sect. 5.

2. Mass-dependent structure of aEMS models

2.1. Preparatory aEMS model calculations

Our purpose is to create mass-dependent structures of aEMS to evaluate collision-induced mass loss. We used the 1D stellar evolution code Modules for Experiments in Stellar Astrophysics (MESA, Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023)². Since in most GCs 2P stars are hardly enriched in He with respect to 1P stars, aEMSs are expected to acquire mass and eject the nucleosynthesis products of hot H-burning very early on in the main sequence. Accordingly, we assumed a fast gas accretion path, starting from a 0.7 M_⊙ seed, until the models reach the desired mass range. We opted for the accreting formulation of Haemmerlé et al. (2019, their Eqs. (1), (2), and (3)), which is based on the Churchwell-Henning (CH) empirical relation (Churchwell 1999; Henning et al. 2000) and scaled it up by a factor of 10. In practice, this implies a mass-dependent mass-accretion rate that increases from ∼10⁻⁴ M_⊙yr⁻¹ (with the advantage of avoiding numerical instabilities at the beginning of the calculations) to 6.15×10⁻² (respectively ∼10⁻¹) M_⊙ yr⁻¹ when the model reaches 10³ (respectively 10⁴) M_⊙, with a mean value of a ∼4×10⁻² M_⊙ yr⁻¹. These values are of the right order of magnitude compare to the aEMS model of Gieles et al. (2025). Mass loss through stellar winds is not included in these preparatory calculations, as the numerical objective was to ensure a fast total mass growth rate of the models while limiting the He-enrichment when the models enter the EMS regime. This is described in Appendix A where we provide more information about the input microphysics of the models and discuss their nucleosynthesis behavior. We also explore the effects of metallicity on the internal structure and nucleosynthesis of aEMS by computing models for three different metallicity values corresponding to well-characterized Galactic GCs ([Fe/H] = −0.72, −1.14, and −2.00 dex, for NGC 104 (47Tuc), NGC 2808, and NGC 6397, respectively).

Since the mass-dependent structure of aEMS is what matters for our collision study, it is important to understand which assumptions made in the preparatory calculations affect the stellar properties. The impact of the gas accretion rate has been extensively studied in the literature (e.g., Hosokawa et al. 2013; Herrington et al. 2023; Haemmerlé et al. 2018; Nandal et al. 2024) and is recalled in brief below. However (and as discussed in Sect. 2.2) for the first time, the treatment of super-adiabatic convection has an important impact on the aEMS properties even for high accretion rates.

2.2. Implications of the treatment of super-adiabatic convection in radiation-dominated stellar layers on the aEMS global properties

The efficiency of energy transport by convection in massive stars has long been debated (e.g., Langer 1997; Maeder 2009). More recently, 3D radiative hydrodynamic simulations have begun to provide insights into the intricate physical processes occurring in this regime. These include the effects of hydrodynamic instabilities and of waves excited by near-sonic turbulent convection (e.g., Jiang et al. 2015; Jiang 2023; Schultz et al. 2020). Such phenomena cannot be adequately addressed within the confines of 1D stellar evolution models. This leads to very uncertain stellar radii in proximity to or above the Eddington limit, where radiative forces significantly influence convective flows and the overall energy transport mechanism (Maeder et al. 2008; Yusof et al. 2013; Köhler et al. 2015; Sanyal et al. 2015; Gräfener 2021; Agrawal et al. 2022a, b).

The challenges described above also apply to aEMSs. These radiation-pressure-dominated objects exhibit super-Eddington luminosities with values between approximately 10⁸ and 10¹⁰ L_⊙, as shown in Fig. 1. There we present key properties of the accreting models with the three values of [Fe/H] that we computed with two different formalisms for the transport of energy by convection in addition to the classical MLT (Mihalas 1978; Paxton et al. 2011). The so-called MLT++ formalism was developed by Paxton et al. (2013) to artificially reduce the super-adiabaticity implied by MLT in some radiation-dominated convective regions of massive stars nearing the Eddington limit. It uses nonlocal modifications to the energy transport equations to better handle some numerical stability issues that arise near the Eddington limit. A modified version of MLT++, here referred to as MLT++L.I. (Jermyn et al. 2023), also artificially enhances the transport of energy by convection when the Eddington factor is high, but incorporates a fully implicit and local formulation. It uses automatic differentiation to improve the stability and performance of the calculations across a wide range of masses and metallicities. As clearly stated by Jermyn et al. (2023), both MLT++ and MLT++L.I. are meant to serve as “stellar-engineering methods to circumvent complex evolutionary stages, rather than a specific physical model that accounts for how convection is modified near the Eddington limit”. Nevertheless, the comparison between the corresponding models that we present below allows for the investigation of the range of uncertainty on the stellar radius and internal structure of aEMSs due to our lack of knowledge about the transport of energy by convection in this regime. This is of particular importance, as it will impact the outputs of the collision-induced processes that we study.

Fig. 1. Evolution of the key properties of the accreting models starting from 0.7 M⊙ and reaching 2×104 M⊙ (the circles, squares and triangles indicate the points at which the accreting models reach a mass of 1000 M⊙ for reference). The models are computed with the same accretion rate and no mass loss, for the reasons discussed in the text.

Panel (b) in Fig. 1 illustrates the monotonic increase in the radius of the MLT model until the aEMS reaches the imposed maximum mass of 2×10⁴ M_⊙. The gas accretion rate is such that the advected entropy cannot be quickly radiated away, and the Kelvin-Helmholtz contraction timescale of the stellar interior is always too small compared to the increase in stellar radius induced by accretion, independently on the metallicity. Consequently, for the three [Fe/H] values we consider the MLT models are predicted to climb the Hayashi track as supergiant protostars that are highly inflated and red (see the blue tracks in the Hertzsprung-Russell diagram (HRD) in panel (a) of Fig. 1). These findings obtained with the MLT framework are in excellent agreement with models with high accretion rates typically above 10⁻² M_⊙yr⁻¹ (as in our models) from the literature (Hosokawa et al. 2013; Herrington et al. 2023; Haemmerlé et al. 2018; Nandal et al. 2024). The treatment of convection in these studies is not specified, which may be indicative of the similarity of their assumptions to those of our MLT models.

The radius of the models computed with the two other convection treatments undergoes inflation, similarly to what is seen for the MLT model during the first ∼30 Kyr (MLT++) and ∼50 Kyr (MLT++ L.I.), reaching up to a radius of ∼2400 R_⊙ and ∼3000 R_⊙, respectively (these numbers are for [Fe/H] = −2 dex and slightly depend on metallicity). However, and as expected, there are substantial differences in the radius and effective temperature (T_eff) of the models, depending on the convection treatments in the regime where radiative pressure is significant. For the [Fe/H] = −2 dex case, this happens for log(L/L_⊙) above ∼5.8, which is reached when the model has a mass of ∼40 M_⊙. Due to the imposed increase in the efficiency of the energy transport by convection, the MLT++ and MLT++L.I. models manage to evacuate the excess energy resulting from the gravitational settling of the accreted material after approximately one thermal time, causing them to deflate to smaller radii despite the relatively high accretion rates we assume. Whatever the metallicity, the MLT++ model is the most compact one, with the highest T_eff, due to the artificial enhancement of the energy transport from the core to the surface. The MLT++L.I. formalism maintains an entropy gradient that is sufficient to drive convection effectively, but not so steep as to cause excessive compactness or elevated surface temperatures. This results in a model that is intermediate in compactness and T_eff between MLT and MLT++.

The uncertainty on the stellar radius in the aEMS regime due to our lack of knowledge on the energy transport efficiency by convection, it will result in some uncertainty in the amount of mass that can be lost by the collision of a star with the aEMS, as discussed in Sects. 3 and 4.

2.3. Implications of the treatment of super-adiabatic convection in radiation-dominated stellar layers on the internal structure and central conditions of aEMS

The Kippenhahn diagrams in radius and mass shown in Fig. 2 for [Fe/H] = −2 dex illustrate the evolution of the internal structure depending on the convection treatment. Upon entering the EMS regime (around 25×10³ years) with continuous expansion of the stellar radius, the MLT model permanently exhibits two convective regions separated by an extended (in radius) radiative zone. The outer convective region is a consequence of the high opacity in most external layers of this model, which has the lowest T_eff and the largest radius. The convective core is driven by H-burning at high temperatures. A large part of the interior of this model is not thermally relaxed, with a positive outward entropy gradient created by gas accretion that prevents the expansion of the convective core. More specifically, since the Brunt-Väisälä N² frequency is directly proportional to the entropy gradient ds/dr, the expansion in radius of the convective core is quenched. While the external and internal convective regions have similar extensions in radius coordinates, the outer region corresponds to only about ∼6% of the star's total mass, whereas the inner region corresponds to about ∼44%. The MLT++ L.I. model initially resembles the MLT model, featuring two convective zones until approximately 6×10⁴ years when the mass is around 2×10³ M_⊙. Subsequently, when the star contracts, the effective temperature and the temperature of the outer layers increase, leading to a decrease in the opacity and convection in the outer layers. However, as the evolution proceeds, the convective core expands both in radius and mass, and the star becomes fully convective once it has reached a mass of ∼10⁴ M_⊙ around 135 kyrs. In the extreme MLT++ case, the star is even more compact and becomes fully convective after 54 kyrs when its mass has reached ∼1.7×10³ M_⊙. In all cases, more than 30% of the mass of the aEMS is convective.

Fig. 2. Kippenhahn diagrams in radius and mass coordinates (left and right respectively) showing the evolution of the internal structure of the accreting models with [Fe/H] = −2 dex as in Fig. 1 for MLT, MLT++L.I., and MLT++ (from top to bottom). The x-axis is the time since the beginning of accretion on the 0.7 M⊙ seed. The color maps show the temperature profiles. Areas highlighted with gray dots denote convective regions.

Concerning the central conditions, the models with the three different convection treatments at a given metallicity follow similar tracks in the T_c−ρ_c plot (Fig. 1, panel (c)), with the more compact MLT++ model reaching slightly higher central temperatures than the other two. This shows that while the convection treatment affects the macroscopic properties of the models, it does not significantly influence the conditions at their center nor the nucleosynthesis at a given metallicity. We refer to Appendix A where we show that the aEMS mass range (10³−10⁴ M_⊙) provides the relevant nucleosynthesis to explain MSP properties (see also Gieles et al. 2025).

2.4. Other sources of uncertainties

Beyond the uncertainties regarding convection, an additional source for the transport of energy and entropy in the radiative layers of the aEMS could be turbulence induced by fast accretion shocks and/or collisions. The 3D radiative-hydrodynamical simulations of the second gravitational collapse following the dissociation of H₂ molecules (Larson 1969) and leading to the birth of a low-mass protostar (1 M_⊙) have shown that strong turbulence develops inside their radiative interior and is maintained as long as accretion precedes onto the protostar (Ahmad et al. 2023; see also the 2D simulations from Bhandare et al. 2020 for protostars up to 20 M_⊙). Whether or not such a strong turbulence could be triggered and maintained in the aEMS remains to be investigated, along with its efficiency in transporting entropy that could potentially moderate the swelling of aEMS. We assume that this phenomena has no impact on our conclusions regarding collision-induced mass loss discussed in the next sections, where we probe the impact of the convection treatment that is the main driver of the swelling efficiency. We note, however, that the presence and efficiency of turbulence in the radiative layers could transport the H-burning products from the core to the external layers, allowing their ejection in the wind, and/or during the peeling of the aEMS due to stellar collisions.

3. The impact of a single collision – Analytical approach

Here, we focus on collision-induced mass loss when a star (hereafter inspiraling star) spirals and deposits orbital energy into the envelope of an aEMS (hereafter target star). We studied the effects of this energy deposition inside the aEMS analytically (i.e., we do not follow the hydrodynamical reaction of the aEMS to the collisions), following the physical criteria described in Sect. 3.1.

In Sect. 3.2, we explain how we used the internal profiles of the relevant physical quantities from the aEMS internal structures of the fast gas accreting models at different masses to infer whether a collision will result in mass loss from the aEMS. We considered three cases for the treatment of convection, but focus the discussion on the models with [Fe/H] = −2.00 dex.

3.1. Analytical prescriptions

We assumed as an initial condition that the inspiraling star is in a circular orbit just below the surface of the target aEMS. Then, as a result of the drag that it feels from the surrounding gas, it starts inspiraling within the envelope of the aEMS. Given the gas properties of an aEMS's envelope and the physical size of a typical star, the dominant drag is expected to be the hydrodynamic one, which is the only drag mechanism we consider. We expect the inspiraling star of mass m_insp and radius r_insp to continue its inward journey as long as two conditions are met: (1) its virial temperature (Siess & Livio 1999; Privitera et al. 2016),

$$T_{\mathrm{virial}}=\frac{\mu m_{\mathrm{H}}G{m}_{\mathrm{insp}}}{k_{\mathrm{B}}{r}_{\mathrm{insp}}},$$(1)

remains higher than the local temperature T_aEMS(r) within the aEMS; then (2) its radius remains smaller than its Roche lobe radius, which we compute as (Eggleton 1983):

$$R_L=r\frac{0.49\hspace{0.5em}{q}^{2/3}}{0.6\hspace{0.5em}{q}^{2/3}+ln(1+q^{1/3})},\hspace{0.5em}q=m_{\mathrm{insp}}/M_{\mathrm{aEMS}}(r),$$(2)

with r being the orbital separation of the system, namely, the distance between the center of the target aEMS and the position of the inspiraling in its envelope (r≤R_aEMS, as we assume that the inspiraling star is in an orbit below the surface of the target aEMS). If either of these conditions is not met at a given depth within the aEMS, we assume that the inspiraling star has merged with the aEMS, transferring all of its mass to the latter. As for the mass-radius relations, we consider the model predictions for the three treatments of convection for both the inspiraling and the target stars. However, for the inspiraling stars, we did not apply the factor of 10 to the CH relation.This is in agreement with the maximum accretion rates for stars with final masses of ∼10² M_⊙ from the inertial flow simulations of Padoan et al. (2020).

The orbital energy that is lost is released locally in the form of heat within the aEMS's envelope. Since we performed our calculations in 1D, we made the implicit assumption that this energy is instantaneously distributed within a spherical shell of the aEMS, around the location of the inspiraling star. If the energy released by the inspiraling star during its orbit is greater than the binding energy of the aEMS's layers above the current position of the orbit, as well as the energy transport timescale from that point to the surface (as reported by our stellar models) is longer than the inspiraling time scale, it can cause the aEMS to lose mass. The orbital energy of the inspiraling star at a given distance r from the center of the aEMS is computed as

$$E_{\mathrm{orb}}(r)=-\frac{{\mathit{GM}}_{\mathrm{aEMS}}(r)m_{\mathrm{insp}}}{2r},$$(3)

where M_aEMS(r) is the enclosed mass of the aEMS at r. The amount of orbital energy released is quantified as $E_{\mathrm{orb}}^{\mathrm{rel}}=E_{\mathrm{orb}}^{\mathrm{ini}}-E_{\mathrm{orb}}$, with $E_{\mathrm{orb}}^{\mathrm{ini}}$ the orbital energy at the surface of the aEMS. The binding energy of the aEMS, denoted as E_bind, is determined as the addition of the potential energy,

$$U=-G{\displaystyle {{\displaystyle \int }}_0^{M_{\text{aEMS}}(r)}}\frac{m\prime }{r\prime }\mathit{dm}\prime ,$$(4)

and its internal energy, I, which is one of the structure outputs of the models.

As mentioned above, both the amount of orbital energy deposited and the corresponding timescale, ${\tau }_{\mathit{dr}}=r/\dot{r}$, are important. If τ_dr is longer than the thermal timescale τ_th, we assume that no mass loss would be expected from the aEMS because of the injected energy can be dissipated efficiently, even in the cases where the energy deposited by the inspiraling star exceeds the binding energy of the aEMS's envelope layers outside the orbit. To compute τ_dr, we ought to consider $\dot{r}$ as proportional to the rate of orbital energy decay attributed to the drag, which we express as

$$\dot{r}={\dot{E}}_d\times {\left[-\frac{1}{2}\frac{{\mathit{Gm}}_{\text{insp}}}{r}\left(\frac{{\mathit{dM}}_{\text{aEMS}}}{\mathit{dr}}-\frac{M_{\text{aEMS}}(r)}{r}\right)\right]}^{-1},$$(5)

where ${\dot{E}}_{\text{d}}\approx F_{\text{d},x}{v}_{\mathrm{cir}}$. F_d,x represents the drag force experienced by the inspiraling star along the direction of its orbital motion. This can be written as

$$F_{\text{d},x}=\pi r_{\text{insp}}^2\rho v_{\mathrm{cir}}^2,$$(6)

where v_cir denotes the circular velocity of the orbiting star at that radius, namely,

$$v_{\mathrm{cir}}={\left[\frac{G(M_{\text{aEMS}}(r)+m_{\text{insp}})}{r}\right]}^{1/2}.$$(7)

In Fig. 3 we compare the absolute value of the binding energy |E_bin| within the interior of a 10³ M_⊙ aEMS model computed with MLT to the absolute value of the orbital energy released $|E_{\mathrm{orb}}^{\mathrm{rel}}|$ by a 30.0 M_⊙ inspiraling star along its journey towards the aEMS's center (the model is for [Fe/H] = −2.00 dex). We indicate the depth at which the orbital energy released is higher than the binding energy of the outer layers of the aEMS, as well as the depth at which the merger can occur due to the condition of the local temperature (T_aEMS>T_vir) or because the inspiraling star fills its Roche lobe (R_L<r_insp). The grey area indicates where these three conditions are met and shows the amount of mass that might be removed from the aEMS. For the specific case shown in Fig. 3, we expect that the mass loss caused by the collision is about ∼14% of the mass of the aEMS, which corresponds to around 141.4 M_⊙. This inspiraling star can effectively remove more than four times its own mass from the aEMS. Subsequently, we examine the mass loss anticipated under the aforementioned assumptions for a diverse range of inspiraling and target stars, and the influence of the mass-radius relations derived from the distinct approaches to convection.

Fig. 3. Comparison between |Ebin| within an aEMS of 1000 M⊙ and the orbital energy $|E_{\mathrm{orb}}^{\mathrm{rel}}|$ of a 30 M⊙ inspiraling star (blue full and orange dotted lines respectively), for models computed with MLT and [Fe/H] = −2.00 dex, with the abscissa in mass and radius in the left and right panels, respectively. The red-dotted lines show the coordinates where $|E_{\mathrm{orb}}^{\mathrm{rel}}|$ is higher than |Ebin|. The green-dotted lines represent the points where the local temperature of the aEMS exceeds the virial temperature of the inspiraling star. The dotted-purple lines represent the points where the radius of the inspiraling star fills the Roche lobe.

3.2. Collision-induced mass loss and mass gain predictions for a large mass range of inspiraling and target stars

We applied the above analytical prescriptions and criteria to different combinations of masses for the inspiraling and the target stars assuming the value of [Fe/H] = −2.00 dex for their metallicity. We provide the predictions in tabular form for use by the community performing N-body and hydro-simulations of star clusters. For this reason, we consider inspiraling stars with masses m_insp between 0.1 M_⊙ and 500 M_⊙ and targets with masses M_aEMS³ ranging from 140 M_⊙ to 20 000 M_⊙. Our analytical analysis is applicable only if m_insp/M_aEMS<1. Actual calculations were made with mass steps for the inspiraling stars of 0.1 and 5.0 M_⊙ for masses between 0.1 and 10 M_⊙ and between 10 and 500 M_⊙ respectively. We derived the values of mass loss and net mass gain over the entire domain by linear interpolation.

The results are shown in Fig. 4 for the models with [Fe/H] = −2.00 dex and for the three different implementations for the transport of energy by convection. In the upper row, the colors scale with the amount of mass that can potentially be lost by the target of a given mass as a function of the mass of the inspiraling star. In the lower row, the colors show the corresponding net mass-change M_net=m_insp−M_loss for the target, with the blue and red regions corresponding respectively to a net mass loss (M_net<0) or mass gain (M_net>0).

Fig. 4. Color maps of the mass loss (upper row) and net mass change (lower row) due to a collision with inspiraling stars with minsp between 0.1 and 500 M⊙ as a function of the mass of the target stars with MaEMS between 140 and 20 000 M⊙. The three panels in each row correspond to the different convection treatments, as indicated. The hatched area delineates the region where our analytical approach is not applicable, namely, minsp/MaEMS≥1. In the bottom row the black lines indicate the change in the sign of Mnet.

Starting the analysis with the MLT assumption, which produces more extended stars, we observe several key patterns in the M_aEMS-m_insp plane. The white area in the upper left panel of Fig. 4 corresponds to cases where no mass loss is expected because the orbital energy released by the inspiraling stars is insufficient to overcome the binding energy of the target aEMS. Consequently, the minimum mass of an inspiraling star required to induce mass loss in the target aEMS increases with the target's mass. For an aEMS target of 1000 M_⊙, the minimum inspiraling mass needed to produce mass loss is approximately 25 M_⊙, resulting in a mass loss of about 126 M_⊙ (i.e., 12.6% of the mass of the target). For a target aEMS of 5000 M_⊙, the minimum inspiraling mass needed to cause mass loss is around 90 M_⊙, with a corresponding mass loss of approximately 375 M_⊙ (i.e., 7.5% of the mass of the target). There is a limited area in the middle part of the M_aEMS-m_insp plane where the amount of mass possibly lost by a target star exceeds 20% of its mass. Within this area, for example, an inspiraling star of 100 M_⊙ can induce a mass loss of 25% in a 1000 M_⊙ aEMS. In the lower left panel, we see a large (blue) area of the parameter space where the collisions are expected to induce a net mass loss from the targets (M_net<0). This region is surrounded by two (red) areas where the target star actually increases its mass after the collision (M_net>0).

In the cases of MLT++ L.I. and MLT++, the aEMSs are more compact and have higher binding energy, making mass loss through collisions more difficult. As in the MLT case, in the upper panel, the white area located at the bottom part of the M_aEMS-m_insp plane corresponds to cases where the collisions are not expected to induce mass loss, because the orbital energy of the inspiraling is insufficient to overcome the binding energy of the target aEMS. Furthermore, it is expected that no mass loss will occur in the white area located at the upper part of the plane. This is due to the Roche lobe criteria, which states that more compact aEMS stars result in smaller orbital distances in the orbit of the inspiraling star. Consequently, it is easier for the inspiraling stars to fill it according to Eq. (2). In this case, we assume that the mass is transferred entirely to the aEMS and that the collision does not induce mass loss from the aEMS.

In particular for the MLT++ L.I models, the radius of the aEMS model reaches a maximum when its total mass is around 1000 M_⊙, and it decreases later on (see panel (b) of Fig. 1). At this stage, the minimum inspiraling mass required to induce mass loss is approximately 30 M_⊙, which can unbind about 5% of the aEMS's mass. Subsequently, as the aEMS becomes increasingly compact, it becomes more challenging for collisions to induce significant mass loss. For instance, for an aEMS with a mass of 2000 M_⊙, the minimum inspiraling mass required to produce mass loss increases to approximately 70 M_⊙. For aEMSs with masses between about 2500 and 6000 M_⊙, inspiraling stars can induce a mass loss of up to 5%. However, for aEMSs with masses exceeding 6000 M_⊙, the mass loss cannot exceed 0.5%. As a result and as can be seen in the lower middle panel of Fig. 4, the area where the target has a lower mass after the collision (i.e., M_net<0) is strongly reduced compare to the MLT case.

In the MLT++ case, the aEMSs are even more compact, resulting in higher binding energy. This compactness makes it more difficult for inspiraling stars to induce mass loss. This considerably reduces the region in the M_aEMS-m_insp plane where collisions are expected to induce mass loss. Consequently, and as illustrated in the lower right panel of Fig. 4, the mass injected by the inspiraling star is always higher than the mass loss it induces from the target whose mass always increases after the merger.

This analysis underlines the necessity for a comprehensive treatment of the transport of energy by convection in aEMS that are radiation-pressure-dominated objects. However, our grid of models and the predictions for mass-loss and mass-gain induced by collisions over a large range of masses for both inspiraling and target stars provides for the first time an opportunity to explore how the stellar M-[Fe/H]-R-structure relation and its uncertainties impact the effects of collisions in the simulations of star cluster. Our grid can directly be used in N-body and/or hydro simulations of star clusters.

4. The impact of multiple collisions in a dense stellar cluster

In this section, we test our predictions for individual collisions in a simplified and illustrative way based on the case of a static proto-GC cluster. Our aim is to gain a deeper understanding of the circumstances under which collisions and mergers might lead to significant mass loss from the aEMS, as well as those where an aEMS in a typical proto-GC might enter into a conveyor belt regime where its mass gain is balanced by its mass loss. As in Sect. 3, here we focus on the models with [Fe/H] = −2.00 dex.

4.1. Model assumptions

We followed the evolution of the mass of an aEMS target of initial mass of 1000 M_⊙ that undergoes multiple collisions with individual stars in a dense star cluster. We considered a static proto-globular cluster with a total initial number of stars of 10⁵ in the mass range between 0.1 and 120 M_⊙. The masses of these stars were sampled according to the initial mass function (IMF) proposed by Kroupa (2001), leading to a total stellar mass in this mass range around M₀≃10⁵ M_⊙. We assumed a half-mass radius of 0.2 pc, leading to a stellar density ρ_h0≃1.5×10⁶ M_⊙/pc³ and a stellar mass surface density of Σ_h≃7×10⁵ M_⊙/pc². This falls within the range of values derived for the candidate proto-GC at high redshift (Sect. 1).

Without any preconceived hypothesis about its formation, we assumed that the aEMS is located at the center of the cluster, it has the structure of the aEMS as predicted in Sect. 3.2, and it is the only target with which the stars can collide (we discuss the consequences of this latest assumption in Sect. 5).

Following each collision, we modified the actual mass of the aEMS taking into account the expected collision-induced mass loss and gain predicted as described in Sect. 3.2 as well as the amount of mass the aEMS has lost by stellar winds and gained through accretion since the last collision (see also Sect. 4.3). For the mass loss due to the stellar winds, we followed the prescription from Vink (2018), namely, $log{\dot{M}}_{\mathrm{winds}}=-9.13+2.1log(M/M_{\odot })+0.74log(Z/Z_{\odot })[M_{\odot }{\mathrm{yr}}^{-1}]$. Contrary to Sect. 2.1, where the assumed mass-dependent high accretion rate to obtain the structure of the aEMSs, we now consider different values for the constant gas accretion rates during the collision phase, namely, Ṁ_acc,gas= 10⁻⁴, 10⁻³, 10⁻², 10⁻¹, and 1.0 M_⊙ yr⁻¹. As for the next collision, we consider the mass-radius relation and the structure predicted in Sect. 3.2 for the aEMS with the new mass calculated as described previously. We derived the collision rate and incorporated it into a Monte Carlo-based framework, as described in Sects. 4.2 and 4.3, respectively.

4.2. Collision rate in a static cluster

Following Hills & Day (1976), the collision rate of the target aEMS with cluster stars is given by

$$\dot{N}(M_{\text{aEMS}},m_{\text{insp}})=n_{\text{insp}}\gamma (M_{\text{aEMS}},m_{\text{insp}}),$$(8)

where n_insp is the average number of stars per unit of volume in the proto-globular cluster, and γ is the rate coefficient with a star of specific mass m_insp defined as

$$\gamma (M_{\text{aEMS}},m_{\text{insp}})=\frac{2{\mathrm{\Sigma }}_0}{\sqrt{\pi l^2}}\left(1+l^2{v}_{\infty }^2\right).$$(9)

Here, Σ₀ accounts for the geometric aspect of the effective cross-section,

$${\mathrm{\Sigma }}_0=\pi (R_{\text{aEMS}}+r_{\text{insp}}{)}^2,$$(10)

and v_∞ is the relative velocity required for two stars to close approach from infinity,

$$v_{\infty }^2=\frac{2G(M_{\text{aEMS}}+m_{\text{insp}})}{(R_{\text{aEMS}}+r_{\text{insp}})},$$(11)

with R_aEMS and r_insp their respective radii given from the mass-radius relation obtained from our accreting models. The parameter l is inversely proportional to the dispersion of relative velocities. Assuming a static proto-globular cluster where all stars share the same mean square speed (i.e., an unsegregated cluster), we set

$$l^2=\frac{3}{4\langle v_m^2\rangle },\text{with}\hspace{0.5em}\langle v_m^2\rangle \simeq 1.8\frac{3\mathit{GM}}{8R_h}.$$(12)

From there, we can introduce the quantity τ(M_aEMS,m_insp), which represents the average time in between collisions involving the aEMS with a star of mass m_insp. It is given by

$$\tau (M_{\text{aEMS}},m_{\text{insp}})=\frac{1}{\dot{N}(M_{\text{aEMS}},m_{\text{insp}})}.$$(13)

4.3. Monte Carlo-based simulations

The following algorithm was employed to manage collisions and their effect on the target aEMS over time. At a time t = 0 of a given simulation, the target aEMS star has a mass of M_aEMS. At each iteration, i, of the calculation, we followed the sequence detailed below:

1. A time step δt of sufficiently small magnitude is selected to minimize the probability of multiple collisions, assuming $\delta t=\frac{\tau (M_{\text{aEMS}},m_{+})}{10}$, where m₊ is the mass of the most massive protostar in the cluster.

2. A mass value m_k is randomly selected from the IMF.

3. The probability of the collision between the target and the star of mass m_k occurring during δt is computed as $p=\frac{\delta t}{\tau (M_{\text{aEMS}},m_k)}$.

4. A random number, u, between 0 and 1 is drawn.

   • (a)

     If u<p, we assume that the collision and the merger have occurred. The target aEMS is assigned a new mass, given by M_{aEMS, i}=M_aEMS,i−1+m_k−M_col,k−M_winds+M_acc,gas and the corresponding radius. M_col,k is the mass lost predicted in Sect. 3.2 for the collision obtained assuming the internal structure of the fast accreting model of mass M_{aEMS, i−1}. Also, M_winds and M_acc,gas are, respectively, the mass lost by the winds and the mass gained from gas accretion by the aEMS during the time step δt. In the last step, a protostar of a specific mass, m_k, is removed from the initial sample of the IMF for the next iteration.

   • (b)

     In the event where u>p, it is assumed that the aEMS is not subjected to a collision during the interval, δt, but still loses mass from winds and gains mass from gas accretion; the protostar of a specific mass, m_k, is kept in the IMF for the next iteration.

The simulation is terminated when it reaches 5 Myrs. This is when the first supernova (SN), stellar feedback, and binary heating effects are expected to impact the cluster properties (see e.g., Gieles et al. 2018, 2025), or when the mass of the aEMS falls below 150 M_⊙ (lower limit) or rises above 2×10⁴ M_⊙ (upper limit). When the aEMS reaches the lower limit, it is assumed that the star will lose the majority of its mass. Conversely, when the star reaches the upper limit, it reaches a state where it becomes extremely challenging to cause further mass loss through collisions. In the MLT case, the minimum inspiraling star mass required to induce mass loss on a 2×10⁴ M_⊙ target is approximately 280 M_⊙, which we did not consider as we have assumed that the highest mass in the IMF is 120 M_⊙. We discuss this limitation further in Sect. 4.4. Finally, given that this is an analytical approach, we do not account for the rejuvenation effects resulting from the collisions, nor the reaction of the star.

4.4. Results of the Monte Carlo simulations

We ran the Monte Carlo simulations 30 times for each combination of the following assumptions: transport of energy by convection (MLT, MLT++LI, and MLT++), and Ṁ_acc,gas over the ranges indicated above. We use the same initial mass distribution along the IMF for the stars for all the simulations. The results are presented in Fig. 5. The x-axis indicates the assumed values for Ṁ_acc,gas (in M_⊙/yr). In the upper panels, the y-axis shows the respective mean values of the contributions of gas accretion and collisions to the total mass gain. In the lower panels, it shows the respective mean values of the contributions of winds and collisions to the total mass loss.

Fig. 5. Top: Respective contributions of gas accretion and collisions to the total mass gain of the 1000 M⊙ target for different gas accretion rates in the Monte Carlo simulations (mean values). Bottom: Respective contributions of winds and collision-induced mass loss to the total mass lost for the same simulations. The time above the bars shows the mean time to the end of the simulations, and the mean number of collisions is displayed in the corresponding bar chart.

A first examination of Fig. 5 indicates that the total quantity of material lost by the aEMS target and the respective contributions of the collisions and of the winds, are primarily influenced by the structural characteristics of the aEMS resulting from the treatment of energy transport by convection. For a given treatment of convection, the second crucial parameter is the rate of gas accretion during the collision phase. When this quantity is varied across the three treatments of convection, it is found that there are configurations where the total mass lost by the aEMS target exceeds its initial mass.

In the MLT case where the aEMS target has the most extended structure, its effective cross-section is higher (see Eq. (10)), leading to the fastest rate of collisions for a given cluster configuration. Consequently, and even if the individual collisions are relatively efficient at inducing mass loss (see Fig. 4), the mass of the aEMS increases steeply thanks to the mergers. All the “MLT simulations” actually stop before ∼5 Myrs, because the target star has reached the maximum mass of 2×10⁴ M_⊙.

The higher the gas accretion rate, the faster this happens (from ∼10 to ∼68 Kyrs for gas accretion rates of 1.0 and 10⁻⁴ M_⊙.yr⁻¹ respectively), with the mass gain of the aEMS being dominated by the mergers except for the highest value of the gas accretion rate (i.e., 1.0 M_⊙.yr⁻¹), where the gas accretion starts to become dominant. Further calculations were conducted for Ṁ_acc,gas = 0.01 M_⊙. yr⁻¹, and an upper mass limit of 4×10⁴ M_⊙, which was reached in ∼0.5 Myr. Regarding the total mass loss, over these short timescales the collisions contribute more to the total ejecta than the stellar winds. Notably, configurations in which the total mass lost exceeds the initial mass of the target occur predominantly in the low and intermediate gas accretion regimes, with the exception of the case where Ṁ_acc,gas = 1.0 M_⊙. yr⁻¹.

In contrast, a more compact aEMS structure, obtained with MLT++L.I. and MLT++, results in a lower effective cross section for collisions. Furthermore, individual collisions are less efficient in inducing mass loss (Fig. 4).

For the two highest gas accretion rates considered (Ṁ_acc,gas = 10⁻¹ and 1.0 M_⊙/yr), the aEMS reaches the upper mass limit in about 19–200 Kyrs. The total mass lost is approximately 200–2300 M_⊙.

On the other hand, all the MLT++L.I. and MLT++ simulations with the three lowest gas accretion rates exhibit a distinctive “conveyor-belt” behavior. This is shown in Fig. 6 where we plot for the MLT++L.I. case the evolution of the mass of the target aEMS for the different Ṁ_acc,gas. During the first ∼1−1.5 Myrs, the mass of the aEMS target increases steadily, mainly due to the collisions and mergers with the stars, as shown in Fig. 5. After this initial mass-rising phase, the aEMS reaches a quasi-equilibrium state where the rates of accretion and mass loss balance each other. The equilibrium mass depends on the value of Ṁ_acc,gas, a lower value of Ṁ_acc,gas resulting in a lower value of the equilibrium mass (∼7500−8500 M_⊙ for the two lowest values of Ṁ_acc,gas, and ∼16 000 M_⊙ for the intermediate one, see Fig. 2 for the MLT+L.I. case). In these configurations, the corresponding evolution models are predicted to be fully convective. The chemical composition of the ejected material should then reflect that of the convective core, as expected to explain MPs in old GC and the peculiarities of the N-emitters. During the conveyor belt stage, the rate of mass gain (due to mergers and accretion of gas) and the rate of mass loss (due to collisions and winds) is approximately ∼4.45×10⁻³, 5.35×10⁻³ and 2.08×10⁻² M_⊙/yr, for Ṁ_acc,gas = 10⁻⁴, 10⁻³ and 10⁻² M_⊙/yr, respectively. The equilibrium is maintained until the end of the simulation.

Fig. 6. Trend of the evolution of the mass of an aEMS with an initial mass of 1000 M⊙ ([Fe/H] = −2 dex) for the different accretion rates for MLT++L.I. The shaded regions enclose the 30 MC realizations, and the solid line shows one realization.

This conveyor-belt behavior results in higher total mass loss compared to the results obtained with the MLT assumption. The maximum value for the total mass lost obtained over the range of parameters investigated, is approximately 4×10⁴ M_⊙, for Ṁ_acc,gas = 10⁻² M_⊙/yr within an interval of 5 Myrs. Importantly, the amount of processed material liberated by the aEMS in the conveyor belt cases depicted in Fig. 6 is one order of magnitude higher than its initial mass.

In the conveyor belt scenario, the aEMS is supposed to be continuously rejuvenated with fresh hydrogen that is both accreted and brought by the merging stars that have not yet burned H (a large fraction of them have a Kelvin-Helmholtz time longer than 10 Myrs). This continuous injection of unburnt material from gas accretion and stellar mergers is expected to effectively sustain a main sequence-like structure for a longer period, delaying the aEMS from evolving along the main sequence and entering advanced evolutionary phases (see also Appendix A.2). While our Monte Carlo simulations indicate that the conveyor belt regime is reachable within the simplistic framework of our static proto-star cluster model, exact numbers for the mass budget should be considered with caution. One of the main caveats is the assumption that the aEMS is the only target for collisions, with more realistic models showing that VMS can be formed by gas accretion and collisions in compact clusters (Gieles et al. 2025, and references therein). Additionally, recent N-body stellar collision studies have shown that the mass function of the colliding stars is more top-heavy than the Kroupa IMF due to mass segregation in the simulated star clusters (Rantala et al. 2024). The next obvious step is to incorporate our model predictions for collision-induced mass loss and gain into realistic high-resolution hydrodynamical and N-body simulations of the formation of massive clusters and their host stellar populations (including remnant BHs).

Another limitation of our analysis is the assumption of circular orbits for the inspiraling stars, while the collisions typically expected in a cluster occur along hyperbolic orbits. This assumption might lead to an underestimation of the velocities and energies involved in the collisions. Furthermore, in all three cases; MLT, MLT++, and MLT++ L.I, the timescale between collisions is shorter than the Kelvin-Helmholtz timescale. This implies that the aEMS does not have enough time to radiate away the energy gained from each collision before the next one occurs. In the future, it will be necessary to evaluate the dynamical response of the aEMS following a collision.

5. Summary and conclusions

Two main channels can contribute to the formation of EMS in non-metal free environments, namely, gas accretion and collisions, with expected co-lateral effects such as those our work aims to investigate for the first time. We proposed and applied analytical prescriptions and criteria to derive the mass loss induced by individual collisions when an inspiraling star spirals and deposits energy into the envelope of a target aEMS.

We considered the internal structure properties of the targets that we obtained for fast accreting stellar models computed with the public stellar evolution code MESA, using different treatments of super-adiabatic convection. While we cannot claim this is how aEMS form, our choice for the gas accretion rate to build the preparatory aEMS models was driven by the necessity to reach extreme mass before burning too much H, as required by the MSP abundance properties. The global properties we obtain for the aEMS as a function of stellar mass agree with those with similarly high accretion rates from previous, studies where the impact the gas accretion rate on the M-R relation has been covered. We show that the predicted M-R and corresponding structures also strongly depend on the treatment of super-adiabatic convection. The more efficient this mechanism is assumed to be, the more compact and the bluer the model is at a given mass for a given gas accretion rate. This highlights the critical necessity for a comprehensive understanding of the transport of energy by convection or other means, including turbulence in radiative layers with high entropy gradient caused by gas accretion, to reliably predict the evolution path of accreting aEMS in the HRD and the extent of their envelope and of their convective core. The question of what happens once the aEMS stops being rejuvenated either through gas accretion and/or collisions and evolves in more advanced phases is out of the scope of the paper, since MSP show only very modest He enrichment and no sign of He-burning yields. Our analytical method can, however, be applied in the future to evaluate collision-induced mass loss on the structure of more evolved stars.

For the assumed mass- and luminosity-dependent gas accretion rates that we apply to reach the mass of the aEMS targets (later used to investigate collision-induced mass loss), the central temperature and density of the models are only modestly affected by the treatment of super-adiabaticity. These quantities depend mostly on the metallicity of the initial chemical composition of the stars, which is assumed to be similar to that of the gas of the proto-GC, and on the actual mass of the proto-aEMS. Consequently, more metal-poor conditions, as well as more massive aEMS, lead to more extreme abundance variations of the isotopes involved in the CNO cycle and the NeNa and MgAl chains in the convective core of the models. This is consistent with the extensions of the C-N, O-N, Na-O, and Mg-Al anticorrelations among MPs in old GC that are observed to depend primarily on the total mass and metallicity of the host cluster. The gas-accretion timescales of our models being much shorter than the nuclear burning timescales, the He-enrichment resulting from this hot H-burning nucleosynthesis is very limited, although slightly higher in the most metal-poor models, in agreement with photometric derivations of He abundance variations among MPs in old GC.

The analytical approach we have proposed to quantity collision-induced mass loss enabled us to investigate a wide range of masses for both the inspiralings and the targets (0.1–500 M_⊙ and 140–20 000 M_⊙, respectively) with three different values of [Fe/H] covering the metallicity range of old Galactic GC. We probed the consequences of the different assumptions about superadiabaticity on the estimation of the mass loss induced by single collisions. In all cases, we found no disruptive collision configuration. As anticipated, the expected mass loss due to the deposition of orbital energy is higher in the case of the most extended targets that have less bound envelopes. It can reach ∼20−25% of the mass of the target star in a limited area of the M_aEMS-m_insp plane, with a minimum m_insp to induce mass loss that increases with M_aEMS. Inducing mass loss from more compact targets, such as those obtained when assuming highly efficient energy transport through convection, proves to be more challenging for two principal reasons. Firstly, because of their higher bounding energy, and, secondly, because the smaller orbital separation reduces the Roche lobe radius, which is rapidly filled by inspiraling stars across a significant portion of the M_aEMS-m_insp plane. These results appear to be qualitatively reasonable. However, the individual collision analysis was conducted using pre-computed target structures, without accounting for the hydrodynamic response of the target to the deposition of energy that might also induce the expansion of the star as well a internal mixing beyond classical convection.

To gain a comprehensive picture of the dynamics of the collisions and accurately quantify their impact in terms of mass loss, it is now necessary to use hydrodynamic simulations that can track the reaction and the structural readjustment of the target aEMS to single and multiple collisions. It should be noted, however, that hydrodynamic simulations of stellar inspirals, even in 1D (e.g., Fragos et al. 2019) are computationally challenging and expensive. Our analytical approach is thus very powerful as it covers a vast range of parameters and can provide guidance on the most relevant configurations for hydrodynamical investigation. Furthermore, this study provides predictions that can be used to investigate the impact of multiple collisions on a target aEMS put in simulations of a more realistic clustered environment. We have made them publicly available together with the mass-radius relations of our accreting star models computed with the three different assumptions for superadiabaticity.

Finally, we have studied collision-induced mass loss and mass gain from an aEMS and quantified its contribution to the mass budget for forming MPs in a proto-GC. We used the predictions for single collisions in Monte Carlo simulations of multiple collisions onto an initially 1000 M_⊙ aEMS target located at the center of a compact and dense star cluster. We assumed the target is surrounded by stars with masses between 0.1 and 120 M_⊙ populating the Kroupa IMF. We chose a number of stars and a cluster radius that resulted in a stellar mass surface density comparable to that determined for candidate proto-GC at a high redshift. We ran 30 Monte Carlo simulations for each combination of the following parameters, applying a super-adiabaticity treatment, while the value of the constant gas accretion rates (between 10⁻⁴ and 1.0 M_⊙/yr) during the collision phase. We could map the mass loss from collisions and winds versus the mass gain from mergers and accretion, as well as the timescales involved, over this range of parameters. The uncertainty about super-adiabaticity has again important consequences, together with the assumption about the gas accretion rate. Target stars with the largest radii are subject to the highest collision rate and reach the upper limit of 2×10⁴ M_⊙ that we have imposed in our simulations in a very short timescale (10⁴–6×10⁴ yr). The total mass lost by the target is limited, even though it can surpass the initial mass of the target in a specific area of the parameter domain we explored. On the other hand, a more compact target is submitted to a lower collision rate. In this case, when the gas accretion rate is low or moderate, the target attains a relatively high quasi-equilibrium mass through mergers (up to 16×10³ M_⊙ in 1 or 2 Myrs), whereby the rates of accretion and mass loss nearly balance each other out. Within the most favourable range of parameters, up to 10^4.9 M_⊙ were ejected in 5 Myrs by the conveyor target. Provided that the target is kept at an almost constant level of rejuvenation through the regular supply of fresh hydrogen, which prolongates its lifetime on the main sequence, this is sufficient to address the mass budget issue for MPs in old GC. A caveat of our Monte Carlo simulations is that we did not address how the 1000 M_⊙ target was formed (i.e., we just assumed very high initial gas accretion). Furthermore, we considered that collisions happen only with that specific object. However, all the member stars are potentially subject to collisions with their masses changing accordingly, which might modify the slope of the stellar mass function by favoring the formation of a VMS and of more than one EMS (see Rantala et al. 2024, and references therein). Finally, we considered a static proto-GC, while the very early dynamical evolution of a proto-star cluster and the corresponding collision rate shall be impacted by gas accretion, mass segregation, stellar feedback, and binary heating.

Very recent Monte Carlo star cluster simulations addressed the impact of these parameters as well as that of the initial conditions (e.g., number of stars, cluster compactness, stellar initial mass function, and primordial binary fraction) on the IMF and find collisional regimes that form VMS and EMS (González Prieto et al. 2024; Vergara et al. 2025). The growth of VMSs in massive star clusters has also been predicted in recent N-body simulations and star-by-star hydrodynamical simulations, with important consequences on the formation paths of IMBHs (Rantala et al. 2024; Lahén et al. 2025). However, these state-of-the-art studies do not account for the proto-stellar phase, which is characterized by a high proto-stellar radius and a loose internal structure, and they assume that colliding stars merge without mass loss, unless a CE event is triggered. The next important step is to include our model projections into sophisticated Monte Carlo and collisional N-body simulations, to explore novel rejuvenation pathways for colliding proto-stars and conveyor belt solutions for the mass budget, as well as the resulting stellar and BH mass function, which will help to elucidate the MP phenomenon in GC.

Data availability

The files needed to reproduce the aEMS MESA models the predictions for the mass-loss and mass-gain are available at this repository.

Acknowledgments

We thank the referee as well as N.Lahén and A.Rantala for their useful comments. The authors acknowledge support by the Swiss National Science Foundation (Project 200020-192039 PI CC and Project CRSII5_213497 PI TF). MG acknowledges grants PID2021-125485NB-C22, CEX2019-000918-M funded by MCIN/AEI/10.13039/501100011033 (State Agency for Research of the Spanish Ministry of Science and Innovation) and SGR-2021-01069 grant (AGAUR).

References

Appendix A: Hydrogen-burning nucleosynthesis, mass, and metallicity effects in aEMS models

A.1. Input physics

To probe the metallicity effects on the aEMS internal structure and nucleosynthesis, we computed models with the same assumptions as in Sect. 2.1 for three values of [Fe/H] (-0.72, -1.14, and -2.00 dex) that correspond to those of three well-characterized Galactic GCs, NGC 104 (47Tuc), NGC 2808, and NGC 6397 respectively. For each metallicity, we run models for the three different treatments for convection as described in Sect. 2.2. In all cases, we consider α−element abundance enhancement of [α/Fe] =+0.3 dex and adopt solar mixture from Asplund et al. (2009) for the other isotopes. The chemical composition of the accreted material is the same as that of the initial proto-stellar seed and it is deposited at the surface of the star. We employ the nuclear reaction network sagb_NeNa_MgAl from MESA, which includes 29 species from neutrons to ²⁷Al, and contains the pp-chains, CNO, NeNa, and MgAl cycles. The equation of state is the standard MESA combination of the SCVH (Saumon et al. 1995), OPAL (Rogers & Nayfonov 2002), HELM (Timmes & Swesty 2000), PC (Potekhin & Chabrier 2010), FreeEOS (Irwing 2004) and Skye (Jermyn et al. 2021) equations of state (Jermyn et al. 2023). We used eps_grav as the standard stellar structure energy equation form (see section 5 of Jermyn et al. 2023). The opacity tables are taken from Iglesias & Rogers (1993, 1996) for the solar mixture of Asplund et al. (2009). At high temperatures (T≥10⁸K), Compton scattering opacities are computed using the prescription of Poutanen (2017). The conductive opacities are taken from Cassisi et al. (2007) with the correction of Blouin et al. (2020).

A.2. Nucleosynthesis

Figure A.1 illustrates the relationship between the increasing central temperature and the stellar mass for all the convection and metallicity combinaisons. The curves initially overlap regardless of the [Fe/H] value and of the treatment of convection, until the stellar mass reaches approximately 10² M_⊙. This is the point where the CNO-cycle reaches equilibrium at the very beginning of central H-burning, causing the core to expand and the central density to decrease (see also Fig. 1). Later on, the more metal-poor models reach a higher central temperature at a given stellar mass, due to the impact of the metals on the stellar opacity. For each [Fe/H], the effect of the treatment of convection appears only when the stellar mass reaches ∼ 500-600 M_⊙, and is more pronounced at lower metallicity. As a result, the central abundance variations of the isotopes involved in the CNO cycle, NeNa, and MgAl chains, as the stellar mass increases in the aEMS model, are more extreme in the model with [Fe/H]=−2.0 dex than in that with −0.72 dex, as shown in Fig. A.2. Dashed gray lines mark the masses when the central temperature T_c reaches specific values at which the abundances of specific isotopes start to change significantly in our models. For both metallicities,¹²C starts to be heavily depleted while ¹⁴N is increasing around T_c∼30 MK. After the models have reached ∼10² M_⊙ and the central temperature is above ∼40 MK, the typical CNO equilibrium patterns build up (this is when core expansion occurs, as discussed above). ²³Na and ²⁶Al also start to increase (from ²²Ne, and ^25,26Mg respectively), while ¹⁶O decreases significantly. Above T_c∼60 MK, namely, for masses above ∼10³ M_⊙ (few 10³ M_⊙) for the [Fe/H]=-2 (-0.74 dex) models, ²³Na and ²⁶Al have increased by more than ∼ 1.5 dex and almost 2 dex respectively. The depletion of ²⁴Mg results in more Al production for T_c higher than ∼70 MK. For [Fe/H]=−2.0 (−0.74 dex) this is reached when the stellar mass is above ∼6×10³ M_⊙ (few 10⁵ M_⊙⁴).

Fig. A.1. Central temperature as a function of the stellar mass for three different metallicities, [Fe/H] =−0.72,−1.14, and −2.00 dex and for the three treatments for convection. The horizontal dashed-gray lines are located at 40MK and 70MK respectively, as a reference.

Fig. A.2. Evolution of the central abundances of relevant isotopes with respect to the mass of the star in the fast accreting models with [Fe/H]=-2 and -0.72 dex (top and bottom respectively; with MLT++L.I. in both cases).

Finally, when the stellar mass reaches 2×10⁴ M_⊙, the helium enrichment is higher in the most metal-poor model (ΔY≃0.07 compared to 0.04), and ²⁴Mg has decreased from 1.18×10⁻⁵ to 2.09×10⁻⁶). However, Mg depletion is not seen in most GC, but only in the most massive and/or metal-poor ones. This means that stars with masses above ∼10⁴ M_⊙ cannot have contributed in most cases (see Gieles et al. 2025 for more details about the Mg-Al anticorrelation). Our results show that accreting models meet the GC nucleosynthesis constraints within a broad mass range, with the most massive and/or the most metal-poor stars leading to the largest abundance variations. The substantial mass content of the aEMS large convective core provides an optimal reservoir for addressing the mass budget issue.

It is important to note, however, that the exact relationship between stellar mass and abundance variations described above should be interpreted with caution. In addition to the influence of metallicity, the actual masses where similar abundance variations are reached shall also depend on the gas accretion rate. For slower accretion, He abundance variations would be higher at a given mass. Because the central temperature increases with the mean molecular weight, which would then increase in the core, the same abundance variations would be reached in slightly less massive stars than indicated above (see Fig. 4 in Prantzos et al. 2017). On the other hand, Na is expected to be depleted if the central temperature is too high. This will be quantified in greater detail in a subsequent study. However, this will not affect the conclusion of the present paper regarding collision-induced mass loss.

All Figures

	Fig. 1. Evolution of the key properties of the accreting models starting from 0.7 M⊙ and reaching 2×104 M⊙ (the circles, squares and triangles indicate the points at which the accreting models reach a mass of 1000 M⊙ for reference). The models are computed with the same accretion rate and no mass loss, for the reasons discussed in the text.
In the text

Fig. 2. Kippenhahn diagrams in radius and mass coordinates (left and right respectively) showing the evolution of the internal structure of the accreting models with [Fe/H] = −2 dex as in Fig. 1 for MLT, MLT++L.I., and MLT++ (from top to bottom). The x-axis is the time since the beginning of accretion on the 0.7 M⊙ seed. The color maps show the temperature profiles. Areas highlighted with gray dots denote convective regions.

In the text

Fig. 3. Comparison between |Ebin| within an aEMS of 1000 M⊙ and the orbital energy $|E_{\mathrm{orb}}^{\mathrm{rel}}|$ of a 30 M⊙ inspiraling star (blue full and orange dotted lines respectively), for models computed with MLT and [Fe/H] = −2.00 dex, with the abscissa in mass and radius in the left and right panels, respectively. The red-dotted lines show the coordinates where $|E_{\mathrm{orb}}^{\mathrm{rel}}|$ is higher than |Ebin|. The green-dotted lines represent the points where the local temperature of the aEMS exceeds the virial temperature of the inspiraling star. The dotted-purple lines represent the points where the radius of the inspiraling star fills the Roche lobe.

In the text

Fig. 4. Color maps of the mass loss (upper row) and net mass change (lower row) due to a collision with inspiraling stars with minsp between 0.1 and 500 M⊙ as a function of the mass of the target stars with MaEMS between 140 and 20 000 M⊙. The three panels in each row correspond to the different convection treatments, as indicated. The hatched area delineates the region where our analytical approach is not applicable, namely, minsp/MaEMS≥1. In the bottom row the black lines indicate the change in the sign of Mnet.

In the text

Fig. 5. Top: Respective contributions of gas accretion and collisions to the total mass gain of the 1000 M⊙ target for different gas accretion rates in the Monte Carlo simulations (mean values). Bottom: Respective contributions of winds and collision-induced mass loss to the total mass lost for the same simulations. The time above the bars shows the mean time to the end of the simulations, and the mean number of collisions is displayed in the corresponding bar chart.

In the text

	Fig. 6. Trend of the evolution of the mass of an aEMS with an initial mass of 1000 M⊙ ([Fe/H] = −2 dex) for the different accretion rates for MLT++L.I. The shaded regions enclose the 30 MC realizations, and the solid line shows one realization.
In the text

	Fig. A.1. Central temperature as a function of the stellar mass for three different metallicities, [Fe/H] =−0.72,−1.14, and −2.00 dex and for the three treatments for convection. The horizontal dashed-gray lines are located at 40MK and 70MK respectively, as a reference.
In the text

	Fig. A.2. Evolution of the central abundances of relevant isotopes with respect to the mass of the star in the fast accreting models with [Fe/H]=-2 and -0.72 dex (top and bottom respectively; with MLT++L.I. in both cases).
In the text