© The Authors 2026

1. Introduction

In dense and massive stellar clusters, the constituents of the system can be segregated by mass. As the central density rises until the collapse of the cluster (Bettwieser & Sugimoto 1984), an increase in close gravitational interactions between stars occurs, possibly leading to the coalescence of these stars (Hills & Day 1976; Portegies Zwart et al. 1999). The dynamical merger of stars is considered to be the scenario behind creating peculiar objects, such as blue stragglers (see Mathieu & Pols 2025, and references therein), as well as very massive stellar objects (VMSOs). The stellar mass spectrum (following the definition by Gieles et al. 2025) given by

• very massive stars (VMSs), with 10² to < 10³ M_⊙;

• extremely massive stars (EMSs), with 10³ to < 10⁴ M_⊙;

• supermassive stars (SMSs), with masses equal or above 10⁴ M_⊙.

Despite the lack of direct observational evidence of the existence of EMSs and SMSs, considering their presence across the cosmic history is relevant as they are formidable candidates for: intermediate mass black hole (IMBH) progenitors (Portegies Zwart et al. 2004; Reinoso et al. 2018; Volonteri 2010; Devecchi & Volonteri 2009); chemical polluters for proto-globular clusters, evidenced by peculiar chemical abundances in their stellar populations (Denissenkov & Hartwick 2013; Gieles et al. 2018, 2025; Ramírez-Galeano et al. 2025); and explaining the enhanced nitrogen abundances observed in extremely compact star-forming galaxies observed in the early universe (see Charbonnel et al. 2023; Marques-Chaves et al. 2024; Schaerer et al. 2024), as well as the observed features of the compact emitters so called as “little red dots” (e.g., Begelman & Dexter 2025; Nandal & Loeb 2025; Coughlin & Begelman 2024).

Nevertheless, a complete self-consistent theory for the dynamical formation of VMSOs requires understanding the role of mass accretion and loss processes experienced in dense stellar environments. The main relevant processes are: stellar mergers contributing to the mass accretion and loss (e.g., Glebbeek et al. 2009, 2013; Gaburov et al. 2010; Ramírez-Galeano et al. 2025; Gieles et al. 2018); as well as the contribution of mass by gas accretion (e.g., Gieles et al. 2025, 2018) and stellar winds driving mass loss (see Ramírez-Galeano et al. 2025; Gieles et al. 2025; Glebbeek et al. 2009; Higgins et al. 2025).

In particular, collision-induced mass loss has been investigated mainly in a hydrodynamical framework for lower mass stars (i.e., from 1 to ∼50 M_⊙) in the context of triple systems (see Gaburov et al. 2010) and involved in head-on collisions (see Glebbeek et al. 2013); it has been reported that a maximum loss of 10%–25% for the system mass results from these interactions. Those results depend on the configurations involved and the structure of each star. Ramírez-Galeano et al. (2025) have provided a semi-analytical prescription to determine the mass loss induced by collisions for VMSOs created by high cold-gas accretion rates, whereas no hydrodynamical framework has been applied to obtain results for this range of stellar masses.

This work aims to study the resulting mass loss from stellar mergers involving EMSs and less massive stars. Approaching this problem in a complete and self-consistent manner is beyond our capabilities, as this would require 3D hydrodynamic simulations with radiation transport. Due to the wide range of temporal and spatial scales, this task becomes very challenging, if not impossible. Given such limitations, we adopted a 1D framework (adapted from the one used in Fragos et al. 2019) built over the Modules for Experiment in Stellar Astrophysics (MESA) stellar evolution software instrument (Paxton et al. 2010, 2013, 2015, 2018, 2019; Jermyn et al. 2023) and taking advantage of its Riemann hydrodynamic solver to model the response of the EMS structure during stellar circular inspirals. We present how the information gathered from circular inspirals can be used to determine the interaction outcomes across initial orbital configurations. The general methodology on the 1D hydrodynamic and stellar evolution framework used is discussed in Sect. 2, followed by the description of the stellar models considered and the set-up for our simulations in Sect. 3. The results are presented in Sect. 4, along with the final discussion and conclusions in Sect. 5.

2. Methods used to simulate stellar mergers

To model stellar mergers involving an EMS and a less massive star, it is essential to use a framework that captures both the orbital evolution of the system and the structural response of the EMS. As the less massive star becomes engulfed within the EMS envelope, the numerical framework must resolve the inspiral dynamics, account for the work done by hydrodynamical and gravitational drag forces, and follow the resulting changes in the EMS’s internal structure and potential mass loss. For this purpose, we employed the StellarInspiral1D¹ code (Fragos et al. 2019), which couples the orbital evolution with stellar physics to track a star’s reaction throughout the inspiral with another, lower mass star. This code is built as an external module to MESA (Paxton et al. 2010, 2013, 2015, 2018, 2019; Jermyn et al. 2023), extending the physics to track stellar inspirals and common envelope evolution events. Furthermore, we used a modified version of MESA, as some of the relevant physical prescriptions for this project were not available in MESA (particularly in release 23.05.1 used in this work) and cannot be added as an external module (or hook). The general description of the StellarInspiral1D code is described in Sect. 2.1 and the physics we incorporated in MESA is described in Sect. 2.2.

2.1. The StellarInspiral1D code

The StellarInspiral1D framework considers the interaction between two stellar objects by modeling one as a single-star MESA model and the other as a non-evolving spherical object. This approach was considered for cases where an extended star interacts with either a non-extended star or compact object such as a neutron star or a black hole (see Fragos et al. 2019, as an example). In the case of a collision between an EMS and a less massive star, the latter hereafter referred as the companion, we would expect the EMS to have a radius three to four orders of magnitude larger than that of the companion (see Ramírez-Galeano et al. 2025). Thus, the previous assumption also holds for the cases regarding this work. We assume the mass of the companion is negligible with respect to the mass of the extended star, as this framework does not account for the presence of the companion star in the gravitational potential of the system (see Bronner et al. 2024, for a possible approach on how to take into account the gravitational field of the companion).

For simplicity, the companion star was modeled as a non-evolving object and it is necessary to set its mass and radius as fixed quantities, M_comp and R_comp, respectively. The StellarInspiral1D code requires an initial configuration for the system, allowing us to freely choose the eccentricity and semimajor axis of the two-body system. The companion’s trajectory in the r-θ plane is computed in the reference frame fixed at the center of mass of the EMS. The equation of motion that describes the orbit between both stellar objects (see Ginat et al. 2020) is given by

$$\begin{array}{c}\hfill \mu \ddot{\mathbf{r}}=-\frac{G(M_{\text{EMS}}(r)+M_{\mathrm{comp}})\mu }{r^3}\widehat{\mathbf{r}}-{\mathbf{F}}_{\text{drag}}(\mathbf{r},\mathbf{v})+\text{Post-Newt.}\end{array}$$(1)

Here, r is the separation between center of mass the EMS and the one of the companion ($\mathbf{r}=r\widehat{\mathbf{r}}$), M_EMS(r) is the enclosed mass of the EMS within separation r, and μ = M_EMS(r)M_comp/(M_EMS(r)+M_comp) is the reduced mass. The term F_drag corresponds to the drag force that is relevant when the companion star is engulfed in the EMS envelope, while the final term accounts for post-Newtonian corrections of the order 2.5.

The hydrodynamical and gravitational drag force, following Ostriker (1999) and Binney & Tremaine (2008), is computed as

$$\begin{array}{c}\hfill {\mathbf{F}}_{\mathrm{drag}}(\mathbf{r},\mathbf{v})=-\frac{\pi \rho R_{\mathrm{acc}}^2{\mathbf{v}}^2}{\mu }\widehat{\mathbf{v}}\times \{\begin{array}{cc}ln\left(\frac{1+\mathcal{M}}{1-\mathcal{M}}{e}^{-2\mathcal{M}}\right),\hfill & \mathcal{M}>1,\hfill \\ ln({\mathrm{\Lambda }}^2-\frac{{\mathrm{\Lambda }}^2}{{\mathcal{M}}^2}),\hfill & \mathcal{M}<1,\hfill \end{array}\end{array}$$(2)

where v is the velocity of the companion in the considered reference frame, ρ is the local EMS density at r, $\mathcal{M}=|\mathbf{v}|/c_{\text{sound}}(r)$ is the Mach number, defined in terms of the local sound speed, and R_acc is the companion’s accretion radius defined in terms of the geometric and gravitational cross-section radii as

$$\begin{array}{c}\hfill R_{\text{acc}}=\text{max}\{R_{\mathrm{comp}},\frac{2G{M}_{\mathrm{comp}}}{v^2+c_{\text{sound}}^2}\}.\end{array}$$(3)

The Coulomb logarithm, Λ, is expressed as a function of the R_acc and the convective core radius of the supermassive star (see Binney & Tremaine 2008; Ginat et al. 2020).

The companion’s trajectory, described by Eq. (1), is solved explicitly within the MESA solver’s time step. The code uses a modern implementation of the Dormand-Prince 8(5,3) solver² based on the algorithm by Hairer et al. (1993). The MESA solver time step in our simulations is restricted to a maximum of 0.1 yr to prevent abrupt changes in the system’s orbital configuration.

The work done by the hydrodynamical drag force, F_drag  •  Δr, is deposited as an additional heat source in the EMS’s envelope at the companion’s position within one accretion radius (i.e., from r − R_acc to r + R_acc). While we tracked the change in the orbital energy and angular momentum as a consequence of F_drag, no transfer of orbital angular momentum to the supermassive star was considered in this study.

The StellarInspiral1D is adapted to use MESA’s implicit Riemann solver for Euler equations of hydrodynamics for an inviscid fluid (see Paxton et al. 2018). The advantage of using this solver lies on the ability to resolve contact discontinuities, such as shocks, without requiring explicit artificial viscosity (Toro et al. 1994), contrary to the scheme described in Paxton et al. (2015). Artificial viscosity can damp the EMS envelope’s kinetic energy and compromise the potential mass loss induced by the stellar collision. We used the adaptive mesh refinement (AMR) functionality that MESA provides for the usage of the Riemann solver (Paxton et al. 2015) instead of its default mesh refinement scheme. We adopted the same settings for MESA’s AMR as in Farag et al. (2022), but we targeted 3000 cells for the resolution of the EMS stellar model.

2.2. Extending the MESA code

In the previous section, we discussed how the StellarInspiral1D code extends MESA to consider the interaction between an extended envelope and a compact star. This is possible by taking advantage of the modular nature of MESA. For this study we identified the need for some physical prescriptions to enhance the physical accuracy of simulations such as those considered by the StellarInspiral1D code. Such modifications involve the re-implementation of the flux-limited radiation transport prescription by Levermore & Pomraning (1981), as discussed in Appendix A; limiting convection velocities, as discussed in Appendix B; and extrapolation of opacity tables for low temperature and density regions, as in Appendix C. The modifications described here were performed based on the release 23.05.1 of MESA³.

3. Setting up the stellar merger simulations

In this section, we describe the final structural adjustments applied to the EMS models after their creation, in preparation for the hydrodynamic modeling of the stellar interaction with the companion. We also outline the initial and termination conditions for the simulations.

3.1. Building models of extremely massive stars

To obtain a proxy for the internal structure of the EMSs, we followed the same process as Ramírez-Galeano et al. (2025), building these kinds of stellar models trough high mass accretion rates by considering the ten-fold of the rate prescription for star formation by Haemmerlé et al. (2019), postulated in terms of the works by Churchwell (1999), Henning et al. (2000), and Behrend & Maeder (2001):

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{acc}}=\frac{f}{1-f}·{\dot{M}}_{\mathrm{out}}\times 10,\{\begin{array}{cc}f=1/3\hfill & \text{if}M\le 5{\mathrm{M}}_{\odot },\hfill \\ f=1/11\hfill & \text{otherwise.}\hfill \end{array}\end{array}$$(4)

Here, Ṁ_out is the mass-loss rate driven by the protostar’s bolometric luminosity:

$$\begin{array}{c}\hfill log\left(\frac{{\dot{M}}_{\mathrm{out}}}{{\mathrm{M}}_{\odot }}\right)=-5.28+0.752log\left(\frac{L}{{\mathrm{L}}_{\odot }}\right)+0.0278{log}^2\left(\frac{L}{{\mathrm{L}}_{\odot }}\right).\end{array}$$(5)

The EMS model was assembled starting from a 0.7 M_⊙ pre-main sequence star. The mass was increased via cold gas accretion following Eq. (4) until it reached the desired mass. We considered EMSs with masses of 1000, 3000, and 5000 M_⊙. All our EMS models are computed for a metallicity of [Fe/H] = − 2.10 corresponding to a hydrogen fraction of X = 0.7296 and a metal fraction of Z = 10⁻⁴. In addition, all EMS models are non-rotating.

The accretion rate given by Eq. (4) is used to reach the desired masses for the EMS before important central H-burning occurs (for justification and details see Ramírez-Galeano et al. 2025). We adopted this variable accretion rate over a constant one for numerical reasons rather than considering it as an accurate physical assumption for the actual formation of EMS in massive and compact star clusters. Effectively, we require a prescription that provides a sufficiently low accretion rate when the mass of our model is still small, to avoid convergence issues, while quickly ramping up to higher values as the mass of the model increases. The accreted material has the same composition as the one in the EMS surface.

In this work, we employed the mixing length theory (MLT) of convection by Mihalas (1978) with α_MLT = 1.93 (Choi et al. 2016; Fragos et al. 2023), adopting the same scheme as in the “MLT” models by Ramírez-Galeano et al. (2025). In the latter study, the possible uncertainties on the stellar structure are investigated by other sets of models where the stellar superadiabaticity is artificially reduced given the MLT++ prescriptions by Paxton et al. (2013), Jermyn et al. (2023). The MLT++ schemes represent a numerical engineering approach for MESA to deal with the evolution of massive stars. Such schemes are not compatible with a hydrodynamic solution of the stellar structure. They artificially enhance the convection flux, rendering this type of simulation numerically unstable. The adjustments are not limited by any timescale linked to the hydrodynamic solver and, therefore, we did not produce any models with such schemes. For a more in-depth overview on the EMS models considered here, refer to Ramírez-Galeano et al. (2025) “MLT” models.

3.2. Relaxing the stellar structure

Once the EMS model reaches the desired mass, following Sect. 3.1, the cold gas accretion is halted. In this initial phase, the EMS structure is computed using the hydrostatic MESA solver. Before initiating the collision simulations, the EMS structure must be relaxed to eliminate numerical features inherited from the accretion phase that could otherwise affect the merger dynamics. Additionally, the numerical solver must be switched from the hydrostatic to the hydrodynamic Riemann solver. The relaxation process is performed in two distinct steps:

1. Post-accretion relaxation: accretion is halted. At this stage, the Newton hydrodynamic solver, without artificial viscosity, is used instead of the hydrostatic one. For numerical reasons, we first switched to the Newton hydro solver instead of the Riemann one. As the EMS model is not assumed to be in hydrostatic equilibrium at this point, the EMS is relaxed numerically to remove dynamical perturbations in its structure induced by the change of numerical solver. The relaxation time is not sufficient for the EMS’s thermal structure to adjust or for the star to reach hydrostatic equilibrium again. Dynamical and analytical studies indicate that the typical time between collisions in dense stellar environments, such as protoglobular clusters or nuclear star clusters, is ≲10³ yr (e.g., Gieles et al. 2018; Ramírez-Galeano et al. 2025; Rantala et al. 2024, 2025); this is shorter than the EMS timescale needed to reach hydrostatic equilibrium (i.e., ∼10⁴ yr) for the models considered in this work. Therefore, thermal perturbations remain un-relaxed between successive collisions and the star is expected to remain bloated. To keep the EMS structure bloated, the model was evolved for only 100 years, approximately one-hundredth of its thermal timescale.

2. Structural relaxation: the Riemann solver is activated to describe the stellar structure. Flux-limited radiation transport is enabled to improve the treatment of radiative losses, as well as allowing for the boundary condition to be set by reducing the outer optical depth from 2/3 to 2/3 × 10⁻¹ (see Appendix A for more context). Additionally, convection velocities are limited to ≤0.1 c_sound. To relax any perturbations induced by those changes, the duration of this relaxation step is 50 years (i.e. ≳5 times the EMS’s dynamical timescale depending on its mass).

A major assumption on the “structural relaxation” step is related to the surface boundary condition of the EMS, where τ_surf = 2/3  ×  10⁻¹ is set as the outermost optical depth for the simulated model (the default value of τ_surf in MESA is 2/3). The temperature and pressure for the outer non-simulated portion of the star (i.e., τ < τ_surf) is set by the surface boundary condition based on the Eddington gray atmosphere model (see Weiss et al. 2004, as an example for the derivation of this model). The Eddington gray atmosphere model assumes that the optically thin layers of the star are in hydrostatic equilibrium (i.e., that the bulk acceleration of the gas is zero). As a result, any natural motion or expansion of the stellar surface is artificially suppressed by this assumption. To avoid such artificial damping the outer stellar layer must reach the point where the surface pressure vanishes (i.e., τ → 0). However, for numerical reasons, the selected value of τ_surf ∼ 0.067 is the lowest optical depth tested where no convergence issues occur at the stellar surface. We note that our results involving the motion of the EMS surface are likely to be affected by the aforementioned limitations (see Sect. 4.1.1).

Throughout both relaxation phases and the merger simulations, convection across the front of supersonic shock waves is suppressed. This transport mechanism is disabled when the change in Mach number between adjacent cells exceeds 0.03, as in MESA’s default settings, to ensure the shock face remains radiative. The structural differences in the EMS before and after the relaxation steps are minimal, but performing these steps is crucial to ensure numerical stability during the subsequent merger simulations.

3.3. Initial and termination conditions for the stellar inspiral

After both relaxation phases, we set the initial conditions for the collisions using the StellarInspiral1D framework. We consider companion masses of 10, 15, 20, 30, 40, 50, 60, and 70 M_⊙. The corresponding radii for these companions are obtained from gas accretion models constructed using the unmodified scheme of Haemmerlé et al. (2019), specifically by applying Eq. (4), divided by 10, and adopting the same composition as the EMS models. An overview on the relevant parameters for the companion are presented in Appendix D.

Following similar assumptions to those in Ramírez-Galeano et al. (2025), all configurations begin in an initial circular orbit with a separation set to 0.8 R_EMS(t = 0). This separation was chosen to avoid the companion exit the EMS during the inspiral as the outer layers pulsate due to the orbital energy injected, as this behavior tends to compromise the convergence of the solver.

The companion’s trajectory is followed until one of the following conditions is satisfied:

• The companion fills its Roche lobe, computed using the formula by Eggleton (1983) in terms of M_comp and M_EMS(r);

• The temperature of the EMS at the companion’s position exceeds the companion’s virial temperature, T_vir, comp.

The virial temperature is defined as the characteristic temperature at which the companion’s gravitational potential energy and thermal kinetic energy are in equilibrium, establishing a relationship between its mass, radius, and internal temperature (see Appendix D as well as Ramírez-Galeano et al. 2025 for more context). When either of these conditions is met, the companion is considered to dissolve within the EMS’s interior. This event is hereafter referred to as the merger event. Once a merger condition is reached, the companion star is removed from the simulation, and the orbital evolution is halted. No further energy is deposited into the EMS envelope. For numerical simplicity, the companion’s mass is not added to the EMS envelope.

All the simulations were tracked until they either reached 5000 yr of evolution, consumed approximately 600 CPU hours, or terminated due to solver failure. The results of these simulations are discussed in the following section.

4. Results

This section presents the impact of stellar mergers on the structure and mass loss of EMSs. To establish a baseline for interpreting the results, we begin with a detailed analysis of two representative cases: stellar mergers involving a 1000 M_⊙ EMS with 10 and 70 M_⊙ companions (Sect. 4.1). These examples highlight the key physical processes driving mass ejection and structural changes. Section 4.2 presents the resulting mass-loss estimates from all our simulations. Finally, in Sect. 4.3, given the results of the hydrodynamic simulations methodology, we provide mass-loss estimates for a broader set of initial orbital configurations, as well as a criterion to determine if a two-star interaction will end up leading to a merger or a scattering event, aiming to present results that can be implemented in dynamical simulations of stellar clusters.

4.1. Stellar mergers between low- and high-mass companions: Study of two cases

We first analyzed the two cases where companions with 10 and 70 M_⊙ masses merge with a 1000 M_⊙ star. Both configurations are particularly informative with respect to our understanding the EMS’s evolution through a merger.

4.1.1. Evolution of a 1000 M_⊙ EMS merging with a 10 M_⊙ companion

Let us focus on the inspiral evolution of the merger between a 1000 M_⊙ EMS and a 10 M_⊙ companion star. Figure 1 shows the system’s evolution across time. The top panel of Fig. 1a presents the Kippenhahn diagram showing the evolution of the EMS’s binding energy (showing the value sign(E_bind)×log(max(1, |E_bind/[erg]|))) across its radial coordinate and time. The bottom panel of Fig. 1a shows a similar diagram with the logarithm of the bulk velocity relative to the local escape speed, log(v/v_esc). In both plots, the position of the companion is shown with the solid black line.

Fig. 1. Inspiral evolution of a 1000 M⊙ EMS with a 10 M⊙ companion. Panel a: Structural and dynamical evolution of the 1000 M⊙ EMS during inspiral with a 10 M⊙ companion. Top: EMS Kippenhahn diagram, in terms of the radius coordinate with respect to time, with the signed logarithm of the binding energy encoded by the background color; Bottom: same but for log(v/vesc). In both panels, dotted cyan regions mark convection zones. The solid black line represents the position of the companion, the black dot marks the merger event. Panel b: Zoom-in on the merger event. Top: EMS radius with v/vesc color map, the value of the radius at t = 0 is shown by the dashed red line; Middle: Unbound (Munb) and escaping (Mesc) masses. Bottom: Surface pressure, its initial value is shown by the dashed red line. In all three plots the time when the merger event is met is marked by the black vertical line.

As can be seen in Fig. 1a, most of the EMS remains gravitationally bound, as indicated by the persistent negative binding energy. However, the EMS reacts to the injected energy by pulsating with a period of approximately 6 years, close to its dynamical timescale of 5.8 years (see Clayton et al. 2017, as a similar behavior is experienced in 1D common envelope simulations). These pulsations are evident by tracking the evolution of the stellar surface across time. As the companion approaches the merger condition, the amplitude of these oscillations increases. After this point, orbital energy injection ceases, and the subsequent evolution is governed by the EMS’s hydrodynamic response. The EMS expands, reaching a radius of approximately twice its initial one.

In the bottom panel of Fig. 1a, we see that after the merger condition is reached, the outer envelope accelerates outward, increasing the amplitude of the stellar pulsations. However, after a few decades, the surface motion decelerates, reducing the amplitude of the oscillations. To investigate the mass loss induced by the stellar merger, we define two diagnostic quantities:

• The unbound mass, M_unb, corresponding to the outermost region of the EMS with positive binding energy, satisfying the condition

  $$\begin{array}{c}\hfill M_{\mathrm{unb}}=M_{\mathrm{EMS}}-M_{\mathrm{low},\mathrm{unb}},\text{where}{\displaystyle {\int }_{M_{\mathrm{low},\mathrm{unb}}}^{M_{\mathrm{EMS}}}}{\epsilon }_{\mathrm{bind}}dm\ge 0.\end{array}$$(6)

  Here ε_bind = ε_grav + ε_int + ε_kin is the specific binding energy, which is the sum of the specific gravitational, internal, and kinetic energies, respectively, the injected orbital energy is already considered in the those terms. M_low, unb ≤ M_EMS is the innermost mass coordinate where the stellar surface remains unbound.

• The escaping mass, M_esc, given by material in the outer EMS layers that satisfy

  $$\begin{array}{cc}\hfill M_{\mathrm{esc}}& =M_{\mathrm{EMS}}-M_{\mathrm{low},\mathrm{esc}},\hfill \\ \hfill & \text{where}{\displaystyle {\int }_{M_{\mathrm{low},\mathrm{esc}}}^{M_{\mathrm{EMS}}}}{\epsilon }_{\mathrm{bind}}dm\ge 0\text{and}{\epsilon }_{\mathrm{kin}}\ge {\epsilon }_{\mathrm{int}}.\hfill \end{array}$$(7)

  The last condition ensures that M_esc only considers the portion of the star that is unbound by its bulk motion, namely, one whose velocity reaches or exceeds the local escape velocity. As in the definition of M_unb, the injected orbital energy had already been considered in the binding energy of the star.

It is important to note that M_esc is the mass that is already escaping the system’s gravitational potential, as its kinetic energy is greater that its gravitational one. In contrast, M_unb quantifies the material that can potentially be ejected, as the binding energy of such material can be positive due to its internal energy while its velocity has not exceeded yet the local escape velocity. As the binding energy of the EMS evolves due to the time dependence of the orbital injected energy, both M_esc and M_unb are time-dependent. The main goal of this case study is to determine the values of M_esc and M_unb.

The time evolution of the EMS radius, M_esc, M_unb, and EMS surface pressure are illustrated in Fig. 1b. The top panel shows the EMS radius with the color map encoding the fraction v/v_esc, the middle panel tracks the evolution of M_esc and M_unb, and the bottom panel presents the surface pressure evolution. As the companion approaches the merger condition, both M_unb and M_esc become significant, reaching values between ∼1 and ∼60 M_⊙ for M_unb, and around 0.1 M_⊙ for M_esc. The maximum value reached by the unbound mass is a consequence of the binding energy of the EMS outer layers and the injected orbital energy by the companion. The EMS’s outermost ∼70 M_⊙ have a binding energy that is about two orders of magnitude lower than that of the deeper stellar layers. Moreover, the energy injected by the companion is almost an order of magnitude larger than the binding energy of this outer material. In contrast, the value of the escaping mass remains sensitive to how the EMS responds by expanding its envelope and converting internal energy into gravitational and kinetic energy.

A comparison between the top and middle panels of Fig. 1b allows us to understand the evolution of M_unb and M_esc. The phases during which M_unb > 0 and M_esc > 0 correspond to the periods of high-amplitude oscillations. However, there are also intervals during which both masses drop to zero and the EMS radius contracts. This indicates that some material that had initially reached or exceeded the escape velocity is decelerated and becomes bound again. This behavior is a consequence of the boundary conditions imposed by the Eddington gray atmosphere, whereby the surface pressure is computed assuming a hydrostatic atmosphere. In a more realistic treatment, the surface pressure would be expected to vanish once the local escape velocity is reached. As shown in the bottom panel, however, the surface pressure decreases by only about one order of magnitude from its initial value. The EMS pulsations are damped within a few decades, not only as a consequence of the hydrostatic atmosphere but it is also expected that radial pulsations to be damped as a consequence of the backward differences applied by the implicit numerical integration scheme used by MESA (see, e.g., Heger et al. 1997; Appenzeller 1970 for context on the damping of radial pulsations by numerical integration schemes similar to the one used in this work). In general, the EMS transforms its internal energy, which contains the thermal energy injected by the companion, into gravitational and kinetic energies as it pulsates. Then, the kinetic energy gets artificially damped, losing a fraction of the energy that can potentially drive mass loss. In this simulation, the EMS loses a maximum of ∼0.1 M_⊙ per oscillation period of approximately 6 years, where such an estimate relies on the physical assumptions imposed in our models. We expect that in nature, the EMS would relax its structure and stop pulsating on a timescale comparable to its thermal timescale following the orbital energy injection.

To estimate an upper limit for the total escaping mass driven by EMS pulsations, we estimated the time it takes for the EMS to thermally relax after the stellar merger. Prior to the perturbation by the companion, the EMS luminosity does not vary significantly over a timescale comparable to the inspiral duration. As the lost orbital energy is deposited in the envelope, the difference between the instantaneous EMS luminosity, L_EMS(t), and the initial luminosity, L_EMS(t = 0), represents the rate of energy loss due to the injection of orbital energy. We define the luminosity fraction, f_L, as the average fractional increase in EMS luminosity over the initial value, integrated from the beginning of the simulation to the time the merger condition is reached,

$$\begin{array}{c}\hfill f_L\equiv \frac{1}{L_{\mathrm{EMS}}(t=0)t_{\mathrm{merger}}}{\displaystyle {\int }_0^{t_{\mathrm{merger}}}}\{L_{\mathrm{EMS}}(t\prime )-L_{\mathrm{EMS}}(t\prime =0)\}dt\prime ,\end{array}$$(8)

where t_merger is the time taken to reach the merger condition from the beginning of the simulation. We integrated only up to t_merger rather than the full simulation duration because different configurations have varying termination times. Additionally, the boundary conditions of the EMS atmosphere may result in an underestimated f_L after the merger.

The time it takes for the EMS to relax thermally is given by

$$\begin{array}{c}\hfill {\tau }_{\mathrm{th},\mathrm{insp}}=\frac{-\mathrm{\Delta }{E}_{\mathrm{orb}}(t_{\mathrm{merger}})}{f_L{L}_{\mathrm{EMS}}(t=0)}·\end{array}$$(9)

The expected duration of EMS pulsation induced mass loss is therefore of order τ_th, insp − t_merger, as a portion of the deposited energy is radiated away before the merger condition is reached. For the case involving the 1000 M_⊙ EMS and 10 M_⊙ companion, this time difference is approximately 420 yr − 137 yr = 283 yr. With a mass-loss rate of ∼0.1 M_⊙ per oscillation cycle of ∼6 yr over 283 years, the total escaping mass is estimated to be M_esc ∼ 5 M_⊙. As our EMS structures are based on the models of Ramírez-Galeano et al. (2025), we compared this result to their mass-loss predictions, where no mass loss was expected for the configuration considered here. This hints that the assumptions in their semi-analytical work might underestimate the amount of collision induced mass loss, as discussed in more details in Sect. 4.2.

4.1.2. Evolution of a 1000 M_⊙ EMS merging with a 70 M_⊙ companion

The resulting evolution of the stellar merger between a 1000 M_⊙ EMS and a 70 M_⊙ companion is presented as a representative case for stellar inspirals involving a high-mass companion. Figure 2 shows the evolution for the corresponding system similar to Fig. 1. In the top panel of Fig. 2a, we observe that the binding energy of the EMS is predominantly negative; however, as the companion approaches the merger condition, a portion of the outer envelope becomes unbound. This behavior is more noticeable for this inspiral with respect to the one involving the 10 M_⊙ companion. As the amount of deposited energy is seven times that seen for the least massive companion and the inspiral time is reduced to less than one half for the 70 M_⊙ companion, the resulting energy rate injection ∼15 times larger than in the 10 M_⊙ companion case.

Fig. 2. Inspiral evolution of a 1000 M⊙ EMS with a 70 M⊙ companion until the simulation crashes. Panel a: structural and dynamical evolution of the 1000 M⊙ EMS during inspiral with a 70 M⊙ companion. Top: EMS Kippenhahn diagram, in terms of the radius coordinate with respect to time, with the signed logarithm of the binding energy encoded by the background color; Bottom: Same but for log(v/vesc). In both panels dotted cyan regions marks convection zones. The solid black line represents the position of the companion, the black dot marks the merger event. Panel b: Zoom-in on the merger event. Top: EMS radius with v/vesc color map, the initial radius value is shown by the dashed red line; Middle: Unbound (Munb) and escaping (Mesc) masses. Bottom: Surface pressure, its initial value is shown by the dashed red line. In all three plots the time when the merger event is met is marked by the black vertical line.

As in the inspiral involving the 10 M_⊙ companion, the EMS pulsates in response to the injected energy. The mean period of the oscillation for this configuration is ∼8 yr. In the bottom panel of Fig. 2a, we present the evolution of v/v_esc fraction for the new configuration. As the companion plunges deeper in the EMS envelope, the magnitude of the bulk velocity, and the amplitude of the oscillations increase. Once the merger condition is reached the EMS undergoes a rapid expansion of its outer layers increasing by nearly one order of magnitude in radius. The simulation is unable to evolve further due to numerical convergence issues, but it is likely that if further pulsations occur the amount of mass lost will increase, as in the 10 M_⊙ merger case.

The time evolution of the EMS radius, the resulting unbound and escaping masses, as well as the surface pressure for this simulation are presented in Fig. 2b. The top panel shows the EMS radius and v/v_esc increase as the companion approaches the merger event. When the merger condition is reached, both M_unb and M_esc continue to increase, reaching a final value of ≳10² M_⊙.

In the bottom panel of Fig. 2b, the surface pressure is shown to decrease by three orders of magnitude with respect to its initial value as a consequence of the energy injected by the 70 M_⊙ companion. The low, but not null, value of the surface pressure reached as the EMS expands is one order of magnitude less than the lowest value reached by the inspiral involving the 10 M_⊙ companion. The motion of the outer envelope has not been decelerated completely by the time the simulation stops. Such behavior is interpreted as a consequence of the higher energy injection rate by this high mass companion with respect to the 10 M_⊙ one.

If we consider that the EMS will continue to eject mass (similarly to the case involving the 10 M_⊙ companion), then it is expected that a maximum of ∼100 M_⊙ will be ejected per each 8 yr period. The inspiral thermal timescale for this configuration is τ_th, insp = 592 yr, while t_merger = 60 yr. A total pulsation mass loss can be estimated with this information, but for this configuration its value of ∼6000 M_⊙ exceeds the mass of the EMS. This is a consequence of considering that in each pulsation, 100 M_⊙ will be ejected as in the first one but this result must be interpreted as an upper value estimate. As a first-order approximation, the total escaping mass should not exceed the maximum amount of mass that can be energetically unbound by the injected orbital energy. The theoretical estimate of this mass is obtained by considering that all deposited −ΔE_orb is distributed just at the surface of the unperturbed EMS such that the binding energy of the external layers is exactly zero; for this configuration, this mass is ∼300 M_⊙. We compared our mass-loss estimates with calculations from other works in the following section.

4.2. Resulting mass loss across EMS and companion mass configurations

The unbound and escaping masses are quantities that evolve over time and due to the assumed atmospheric boundary condition in our models and assuming the treatment of other physical processes to be accurate for our simulations, they might artificially decrease. Thus, analyzing their values at the end of the simulations would lead to an underestimation of both quantities. As stated in the previous section, the EMS structure reacts to the injected orbital energy by pulsating, converting the thermal energy into kinetic energy. During this process, the actual mass lost from the EMS corresponds to the mass that escapes the gravitational potential of the star; meanwhile, the energetically unbound mass that has not yet reached the escape velocity is considered to be the portion of the EMS that can be lost in future pulsations. In most of our inspiral simulations, the outer envelope of the EMS exceeds its local escape velocity in the first post-merger pulsation but any subsequent mass loss after this event is damped quickly. In order to obtain a global estimate for the mass loss produced by the stellar mergers, we introduced two quantities:

• the maximum unbound mass, across time, post-merger (M_unb, max) during the inspiral evolution as a proxy for the maximum amount of mass that can be lost;

• the escaping mass produced by the first pulsation post-merger (M_esc, puls), as a proxy of the maximum amount of mass that can be lost from each pulsation of the EMS produced by the injected energy.

In addition, we define the predicted total mass loss produced by the EMS pulsations post-merger as

$$\begin{array}{c}\hfill M_{\mathrm{puls},\mathrm{tot}}=min\{M_{\mathrm{unb},\max },\frac{M_{\mathrm{esc},\mathrm{puls}}}{{\tau }_{\mathrm{puls}}}({\tau }_{\mathrm{th},\mathrm{insp}}-t_{\mathrm{merger}})\},\end{array}$$(10)

where τ_puls is the dominant period of the EMS radius oscillations before the merger event, calculated in terms of the inverse frequency-spectrum from the fast Fourier transform of the EMS radius evolution in such a time domain. As mentioned before, the value of τ_pul is of the order of the EMS dynamical timescale. The previous definition assumes that if the motion of the outer envelope is not damped, then the EMS would produce an escaping mass of M_esc, puls per oscillation period, τ_pul. This would last until the EMS manages to radiate away the total injected energy after the merger condition is reached. In addition, we considered M_unb, max as the maximum mass that can be lost by the EMS and, thus, we set limits on M_pul, tot of M_pul, tot > M_unb, max.

The values of M_unb, max, M_esc, puls, and M_puls, tot are shown in Fig. 3 as a function of companion mass for all the EMS and companion-mass configurations considered in this work. The discussion of this plot is continued in the following sections.

Fig. 3. Maximum unbound mass (Munb, max, purple triangles pointing down), maximum escaping mass (Mesc, max, mint triangles pointing up), and predicted total pulsation mass loss (Mpul, tot, turquoise diamonds) as a function of companion mass for the three EMS models considered. The dashed gray line indicates the one-to-one relation, a predicted mass loss below this line means that the EMS gains more mass by the stellar merger than what it loses as a result of it and vice versa. The dotted blue line shows the maximum unbound mass estimation if there is no loss of injected energy during the inspiral, namely, α = 1 at all times. The dashed black line shows the extrapolated mass-loss prescription from Glebbeek et al. (2013). Analytical predictions from Ramírez-Galeano et al. (2025) are shown with the solid cherry-hue lines.

4.2.1. Comparison with analytical mass loss estimates

To compare our results with the predictions by Ramírez-Galeano et al. (2025), we need to extract the fraction of deposited orbital energy that is used to drive mass loss, defined here as the parameter α. The value of α is determined by the amount of −ΔE_orb that has not been lost yet (e.g., has not been radiated away) by the EMS at t_merger, namely,

$$\begin{array}{c}\hfill \alpha =\frac{\text{total injected energy}-\text{injected energy lost}}{\text{total injected energy}}·\end{array}$$(11)

In the framework by Ramírez-Galeano et al. (2025) the value of α is not calculated explicitly but is determined by the following assumptions:

• The orbital energy deposited in regions where the local thermal timescale τ_{thermal, local}(r) = E_int(r)/L_surf is shorter than the local inspiral timescale ${\tau }_{\mathrm{insp}}=r/\dot{r}$ is promptly radiated away and excluded from the total orbital energy injected ΔE_orb;

• Only the orbital energy deposited above the innermost coordinate, r, where α |E_orb(r)| > |E_bind(r)| is used to unbound the mass above r.

Within such a framework, it is assumed that all the predicted unbound mass will end up escaping. In our simulations, the value of α until merger can be estimated in terms of the amount of −ΔE_orb and by assuming the injected energy can only be lost by being radiated away as

$$\begin{array}{c}\hfill \alpha =\frac{-\mathrm{\Delta }{E}_{\mathrm{orb}}(t_{\mathrm{merger}})-f_L{L}_{\mathrm{EMS}}(t=0)t_{\mathrm{merger}}}{-\mathrm{\Delta }{E}_{\mathrm{orb}}(t_{\mathrm{merger}})}·\end{array}$$(12)

An overview on the comparison on the values of α from the framework by Ramírez-Galeano et al. (2025) and our results is presented in Appendix E.

The mass-loss estimates from the cited analytical framework are shown in Fig. 3 by the solid cherry-hued line. It is observed that our simulations predict higher unbound mass compared to the analytical estimates. This is based on the assumption in Ramírez-Galeano et al. (2025) that no mass can be unbound if the injected orbital energy is lower than the local EMS binding energy (i.e., α = 0). Consequently, no mass loss is expected from mergers involving low-mass companions (with the minimum companion mass to induce mass loss increasing with the mass of the EMS). From our simulations, however, a negligible value of α is produced by systems where t_merger ≫ τ_th, insp as expected from Eq. (12), and if one substitutes the definition of τ_th, insp (Eq. (9)) it is obtained that

$$\begin{array}{c}\hfill \alpha =1-\frac{t_{\mathrm{merger}}}{{\tau }_{\mathrm{th},\mathrm{insp}}}·\end{array}$$(13)

An analytical estimate for the maximum amount of mass that can be energetically unbound is shown in Fig. 3 with the dotted blue line. This estimate is determined by considering the initial EMS profile, following the same arguments as in Ramírez-Galeano et al. (2025), but assuming that α = 1; namely, none of the orbital injected energy is lost during the inspiral and all of it is distributed to set the EMS binding energy of the surface to zero. The resulting M_unb, max from our simulations is overall consistent with this analytical estimate for high-mass companions and significant discrepancies arise for low-mass ones. Such differences are mainly a consequence of not taking into account orbital energy losses during the inspiral.

It is relevant to point out that the resulting values of α from our simulations are dependent on the choice of initial conditions. All of our numerical experiments start from a circular orbit. It is expected that for eccentric configurations, the ratio t_merger/τ_th, insp changes for the same EMS and companion masses, resulting in different values for α. If α gets closer to unity, the maximum unbounded mass estimate serves as a better representation of the possible mass loss from these kinds of stellar mergers.

Regarding the total mass-loss estimation from the pulsation behavior, for the 1000 M_⊙ EMS the resulting M_esc, puls is consistent with the maximum unbound mass. We observed that the match between M_esc, puls and M_unb, max shows a growing discrepancy for the configurations with more massive EMS models. We interpret this to be a consequence of the assumed hydrostatic atmosphere condition, as in the case of more massive EMS models, we would require a higher injected energy to expand the star, and this expansion will end up damped by the hydrostatic atmosphere, leading to an underestimated M_esc, puls and, therefore, an underestimated M_pul, tot.

4.2.2. Comparison with mass-loss estimates from other hydrodynamic simulations

We compared our results to the mass-loss prescription from Glebbeek et al. (2013) based on their hydrodynamic simulations for head-on collisions of lower-mass stars (i.e., from 5 to 40 M_⊙). Since the dynamical configuration, as well as the physical solution of the stellar merger, is different for that work with respect to ours we perform a comparison of both to raise awareness on the impact of the considered physical assumptions on the resulting mass-loss estimates from both studies.

The predicted mass loss reported by Glebbeek et al. (2013) is shown in Fig. 3 as a black dashed line. Such results predict a maximum mass loss of 10% of the total mass of the system, but this high level of mass loss is experienced only when the value of the mass of both stars in the collision are similar (i.e., m_star1 ≃ m_star2). The configurations considered in our study consist of extremely low mass ratios, M_comp/M_EMS ≲ 0.1, where Glebbeek et al. (2013) predicts a mass loss of < 2% of the system mass, M_EMS + M_comp. From our results, the value of M_unb, max represents ≲20% of the system mass, similarly for M_pul, tot in the case of the 1000 M_⊙ EMS. We attributed the high percentage of unbound mass from our simulations and assumed mass loss to the structure of the EMS. The EMS envelope is mostly radiation pressure dominated as the EMS luminosity approaches the Eddington limit; in addition, the EMS is marginally bound as the value of the adiabatic exponent is close to 4/3 across the star. There are no further details available on the structure of the stars considered by Glebbeek et al. (2013) but any differences in terms of the resulting mass loss from both studies might arise from the marginally stable EMS models considered here.

4.3. Possible merger outcomes across eccentricities

Within the limitations of the presented numerical experiments, the final mass loss from a stellar merger depends primarily on the total orbital energy injected, the EMS binding energy distribution, the inspiral duration, and the timescale over which the EMS can radiate the injected energy.

The StellarInspiral1D framework allows for simulations to be performed with initially eccentric orbits. However, tests with such configurations show that the solver timestep is strongly modulated by the net injected energy. In eccentric orbits, the loss of orbital energy increases near periastron. This reduces the timestep of the calculations by several orders of magnitude earlier in the inspiral evolution, with respect to the initial circular orbit configuration. Following these mergers would require 10–100 times more computational time than the models considered in this work, namely, 6000–60 000 cpu hours.

To estimate the possible outcomes of stellar interactions across orbital configurations, we adopted an alternative method, where Eq. (1) is solved assuming a static EMS structure during the companion’s trajectory. This allowed us to determine whether a given orbital configuration leads to a stellar merger and to estimate the orbital energy lost in the interaction that can be used to drive mass loss.

We characterized the initial configuration using eccentricity (e) and semilatus rectum (ℓ), the latter being always positive unlike the semimajor axis (a). The semimajor axis is related to the eccentricity and semilatus rectum as a = ℓ/(1 − e²). Only conics intersecting the EMS radius are considered for the initial conditions. The companion’s trajectory is integrated until completing one orbit or exiting the EMS. Outcomes are classified by the final orbital parameters as

• A.

  the companion never exits the EMS;

• B.

  the companion exits the EMS in the first orbit with e_final < 1;

• C.

  the companion exits the EMS in the first orbit with e_final > 1;

• D.

  the merger condition is reached during the first orbit.

For outcomes A, B, and D, the post-interaction orbit is bound and a merger will occur in the absence of external perturbations. In case B a merger is expected, although in a dense stellar cluster the companion could, in principle, interact with other stars after leaving the EMS envelope and prevent the subsequent merger. The total injected energy for these cases is

$$\begin{array}{c}\hfill -\mathrm{\Delta }{E}_{\mathrm{orb},\mathrm{merger}}=\frac{G{M}_{\mathrm{EMS}}(r_{\mathrm{merger}})M_{\mathrm{comp}}}{2r_{\mathrm{merger}}}-\frac{G{M}_{\mathrm{EMS}}{M}_{\mathrm{comp}}}{2a_{\mathrm{initial}}}·\end{array}$$(14)

Here r_merger is the separation at which the merger condition is met, M_EMS(r_merger) is the enclosed mass at that radius, and a_initial is the initial semimajor axis.

For outcome C, the orbit remains unbound and no merger occurs; this is classified as an orbital scattering. The injected energy is then

$$\begin{array}{c}\hfill -\mathrm{\Delta }{E}_{\mathrm{orb},\mathrm{scattering}}=\frac{G{M}_{\mathrm{EMS}}{M}_{\mathrm{comp}}}{2a_{\mathrm{final}}}-\frac{G{M}_{\mathrm{EMS}}{M}_{\mathrm{comp}}}{2a_{\mathrm{initial}}}·\end{array}$$(15)

To illustrate the previous methodology, Fig. 4 shows, for a 1000 M_⊙ EMS and a 70 M_⊙ companion, the collision outcomes and deposited orbital energy across initial orbital configurations. The EMS profile after relaxation is used for trajectory integration. The semilatus rectum ranges from 10 to 10⁵ R_⊙, the eccentricity from 0 to 2, and the EMS radius is indicated by the dashed gray line. Initially bound orbits (e < 1) remain bound, as drag reduces orbital energy and angular momentum, decreasing the semilatus rectum and eccentricity. For initially unbound orbits (e > 1), mergers require a semilatus rectum ∼100 R_⊙ (∼1/100 of the EMS radius). For e ≈ 1, all allowed initial configurations yield bound orbits, but only A and D guarantee a merger; B requires no external interaction to remain bound.

Fig. 4. Collision outcomes and deposited orbital energy for the interaction between a 1000 M⊙ EMS and a 70 M⊙ companion across initial orbital configurations. The EMS structural profile after the relaxation steps is used for trajectory integration. The initial semilatus rectum ranges from 10 to 105 R⊙ and the eccentricity from 0 to 2. The EMS radius is indicated by the dashed gray line. Left plot: Classification of collision outcomes (A–D; see Sect. 4.3) as a function of the initial semilatus rectum and eccentricity. Right plot: Deposited orbital energy for mergers and scatterings in the absence of external interactions.

It seen from the right panel of Fig. 4, for the mergers without external interactions (cases A, B, and D), the deposited energy is nearly independent of the initial semimajor axis. Thus, mass loss should mainly depend on the rate of orbital energy injection and the EMS response timescale. In C, the injected energy decreases steeply with increasing semilatus rectum; without a merger, mass loss is still possible, but at much lower levels given the amount of injected energy. If in a different scenario, we considered case B as scattering events as only cases A and D ensure a direct merger, then the injected energy for case B will be the same as in case C.

Even though a static profile method expands the accessible orbital parameter space, it ignores the EMS’s reaction to energy injection by successive collisions. As described in Sect. 4.1, the EMS structure will be affected, depending on the companion star’s mass, possibly leading to a more extended and less bound EMS. To illustrate the possible collision outcomes for a successive stellar merger of a 1000 M_⊙ EMS with two companions of the same mass, Fig. 5 shows the collision outcomes for the configuration of a 1000 M_⊙ EMS and 70 M_⊙ companion but using the EMS profile at the end of the hydrodynamic simulation from the merger with the same mass configuration (see Sect. 4.1.2 for the evolution EMS structure during a first merger with also a 70 M_⊙ star). The EMS radius has expanded to ∼5 × 10⁴ R_⊙ (dashed gray line). With this larger structure, initially bound orbits with higher semilatus rectum can be enclosed. For initially unbound orbits, the outcome changes little: direct mergers still require ℓ ≲ 100 R_⊙. Thus, using the initial EMS profile underestimates mergers for extended bound orbits but yields similar results for unbound orbits.

Fig. 5. Same as Fig. 4, but using the EMS profile at the end of the hydrodynamic collision simulation for the interaction between a 1000 M⊙ EMS and a 70 M⊙ to determine the possible outcomes of a second collision. The EMS has expanded to a radius of ∼5 × 104 R⊙, shown by the dashed gray line in both panels. Left plot: Collision outcome classification, right plot: Deposited orbital energy for the same configurations.

With the previous methodology, we can estimate the orbital energy lost during a two star collision. We can also estimate the maximum amount of mass loss by assuming all the lost orbital energy is used to unbind the mass in the EMS outer layers. Such prediction corresponds to value of outer mass where the orbital energy injected is equal to the EMS binding energy integrated from the surface inwards. The mass-loss estimate for mergers depends on the mass of both stars as well on the initial eccentricity and semilatus rectum; whereas for scattering events, the mass loss will depend on the final orbital parameters as well (see Eq. (15)).

To illustrate the predicted mass loss for direct mergers across stellar masses, we considered a fixed initial eccentricity of e_initial = 1.01 to reduce the dimension of the parameter space required to determine such an estimate. The selection of initial eccentricity is representative of nearly parabolic encounters in stellar cluster simulations (e.g., see Rantala et al. 2024, 2025, where ≳90% of stellar merges, given by the sticky sphere approximation, involve eccentricities from 0.9 to 1.1). As we only sought to determine the mass loss for direct stellar mergers (i.e., case A and D), this implies that the initial semilatus rectum is small enough to lead to such a collision outcome. We also determined such a threshold value for the semilatus rectum to produce a direct merger with e_initial = 1.01.

Figure 6 shows the estimates of the mass loss and semilatus rectum threshold for mergers, in units of EMS mass and radius, respectively, as a function of EMS masses (10³–10⁴ M_⊙) and companion masses (5–500 M_⊙). In the left panel, the mass-loss fraction reaches ∼50% for high mass ratios (M_comp/M_EMS → 1), and is negligible for small mass ratios (M_comp/M_EMS → 0). In the right panel, the maximum ℓ_initial for a merger with e_initial = 1.01 reaches at most ∼10% of R_EMS for high mass ratios; in other cases, it depends on the orbital energy loss needed to trap the companion and on its gravitational cross-section, both scaling with the masses involved. From the discussion around the data in Fig. 5, we see that the EMS response can modify the semilatus rectum criterion, but for a high eccentricity (e.g., 1.01,), such variations are assumed to be minimal. For the implementation of the resulting mass loss from this work in dynamical simulations with multiple stellar interactions, we recommend using the semilatus rectum criterion rather than the EMS radius to determine stellar mergers, instead of applying the sticky-sphere approximation for collisions in terms of the EMS radius. For absolute values on the mass loss and semilatus rectum threshold, we refer to Appendix F.

Fig. 6. Relative estimates of merger-induced mass loss and the semilatus rectum threshold between stellar mergers and scatterings, for orbits with initial eccentricity of 1.01. Considering EMS masses from 103 to 104 M⊙ and companion masses from 5 to 500 M⊙. Left plot: Ratio of the maximum merger-induced mass loss to the EMS mass. Right plot: Ratio of the maximum initial semilatus rectum for a stellar merger to the EMS radius.

5. Conclusions

With the aim to understand one of the processes involved in the formation of EMSs and SMSs, this work investigates the role of mass loss induced by stellar mergers. To accomplish that task, we solved the evolution of circular stellar inspirals between massive and extremely massive stars by the StellarInspiral1D framework adapted to the MESA stellar evolution code. The simulations involved 1000, 3000, and 5000 M_⊙ EMSs and inspiraling companions of < 100 M_⊙. We find that the thermal energy deposited by the stellar companions traveling inside the EMS envelope, resulting from the hydrodynamic and gravitational drag force, perturbs the equilibrium structure of the EMS and leads to pulsations. Such pulsations increase in amplitude with the energy that has been injected, while the stellar surface can be accelerated until it exceeds the local escape velocity, at which point mass loss is produced. The ejected mass by the pulsations becomes relevant as the companion star gets closer to the point where it merges with the EMS structure, as there the amount of orbital energy lost is at its maximum.

When the companion stars merge with the EMS structure (they reach the merger condition), the deposition of the lost orbital energy ceases. After that point, the stellar pulsations are usually damped over a few dynamical timescales, namely, on the order of ∼10 yr.

To estimate the total amount of mass that could be lost through the pulsations induced by the inspiral, we extracted the amount of mass lost in the first expansion post-merger. Such a level of mass loss is on the order of ≳10⁻² to 10² M_⊙. The oscillation period was obtained by performing a fast Fourier transform (FFT) of the evolution of the EMS radius over time; we also note that this value does not vary significantly during the evolution of the simulations. The pulsation period does vary across configurations; nevertheless, its value is always of the same order of magnitude as the EMS dynamical timescale.

With the resulting pulsation mass loss and period, we conclude that inspirals involving low-mass companions produce a mass-loss rate from ∼10⁻³ to ≲10⁻¹ M_⊙/yr; the mergers involving high-mass companions produce a mass-loss rate from ≳10⁻¹ to ∼10 M_⊙/yr.

To estimate the duration of the EMS pulsations post-merger, we estimated the timescale for the star to radiate away the total injected orbital energy, as in Eq. (9). We consider that the deposited orbital energy is radiated away at a rate of average increment of the EMS luminosity with respect to its initial value. EMS luminosity increases, on average, from ≳1% for inspirals with lower mass companions (i.e., ∼10 M_⊙), to ≲200% for high-mass ones.

By scaling the resulting mass-loss rate from the pulsations across the time it takes for the EMS to radiate away the injected energy, after the merger, we estimated the total mass loss from the stellar inspirals. The scaled mass loss is consistent with the maximum mass that is energetically unbound across the inspiral evolution. The match between the two values, total ejected and unbound mass, is better for the least massive EMS. For this EMS mass, the resulting damping for the pulsating motion, related to the physical assumptions and the numerical integration scheme, is at a lower level. We conclude that most of the total injected energy is used to drive the ejection of mass (see Sect. 4.2). This last point is explained by the fact that the EMS models considered here have an average adiabatic exponent of ≃4/3, meaning that the stars are almost dynamically unstable and the perturbation by the stellar inspiral can efficiently lead to mass ejection.

In terms of the system’s mass (M_EMS + M_comp), the total pulsation mass loss varies from ∼1% to ≲30% across configurations. Such values are consistent with the predictions by Ramírez-Galeano et al. (2025), with the exception of cases the study would predict no mass loss; however, the hydrodynamical simulations here did indeed produce it.

Given the computational time constraints imposed by the StellarInspiral1D framework, we did not perform simulations for eccentric mergers in this study. Nevertheless, we provide a semi-analytical prescription for such cases. By solving the companion’s equation of motion, Eq. (1), on a static EMS profile, we discriminated between the orbital parameters that lead to stellar mergers and orbital scatterings. Given the initial and final semimajor axis for each interaction, we computed the total orbital energy lost through the work done by the drag force. By assuming all the deposited energy is used to drive mass loss, we can estimate the latter quantity. With this methodology, we find that the amount of orbital energy deposited by stellar mergers does not vary significantly across orbital configurations, implying a similar resulting mass loss for those systems. However, the lost orbital energy is more sensitive to the initial orbital configuration when the companion is scattered by the EMS. Finally, initially bound orbits remain bound after the interaction, whereas initially unbound orbits require a semilatus rectum of ≲0.1 R_{EMS, initial} across eccentricities to produce a merger instead of scattering the companion’s trajectory (further details are discussed in Sect. 4.3).

Our results illustrate the importance of collision-induced mass loss as a fundamental process regulating the formation and growth of EMSs. In the future, an open task is to assess how the cumulative effect of this process influences the formation and long-term growth of EMSs and SMSs. In particular, its role in setting the relevance of EMSs across cosmic time, whether as progenitors of intermediate-mass black holes or linked to other astrophysical scenarios, remains to be addressed in future studies.

Data availability

The data generated for this study are publicly available via Zenodo and can be accessed in the associated repository: https://zenodo.org/records/18075985.

Acknowledgments

This work was supported by the Swiss National Science Foundation (PI Fragos, project number CRSII5_213497). EF acknowledges support by the Yale Center for Astronomy and Astrophysics Prize Fellowship. JRG acknowledges N. Lahén and A. Rantala for their insightful comments on this manuscript, as well as the MESA community and developers for their for the valuable input to aid the progress of this work.

References

Appendix A: Flux-limited radiation transport and surface boundary condition

Having a consistent model for radiative transport valid in the extreme cases of optically thin and thick regimes is necessary to model the extended envelope during the stellar mergers. Then for this work, the flux-limited diffusion (FLD) theory developed by Levermore & Pomraning (1981) has been re-implemented in MESA⁴. This feature was previously available in the code but had been removed in more recent versions. The FLD theory provides a closure relation for the moment equations of the radiation intensity. As well, it provides a description for radiation transport on the transition between the diffusive regime (optically thick, τ ≫ 1) and the free-streaming regime (optically thin, τ ≪ 1). The zeroth and first order moments of the radiation intensity, the radiative energy density (E_rad) and the radiative flux (F_rad) respectively, are related as

$$\begin{array}{c}\hfill F_{\mathrm{rad}}=\frac{1}{3}\frac{c}{\kappa \rho }\lambda \mathrm{\nabla }{E}_{\mathrm{rad}},\end{array}$$(A.1)

where λ is the flux-limiter derived by Levermore & Pomraning (1981), which parametrizes the angular dependence of the radiation beam’s intensity and is a function of the flux-ratio R_fl ≡ |∇E_rad|/κρE_rad. If τ ≪ 1 then R_fl → ∞ and the flux-limiter should asymptotically approach zero, contrary to the case of τ ≫ 1 where R_fl → 0 and λ → 1.

Unfortunately, MESA does not solve the radiation transport moment equations to obtain E_rad and implement FLD theory in a self consistent manner. This implementation of FLD requires a closure assumption on the value of E_rad to compute R_fl and then λ. The assumption is to consider E_rad = aT⁴ everywhere in the star (similarly to Bersten et al. 2011). However, such an assumption is only valid in the optically thick regime, leading to an underestimation of R_fl where τ ≪ 1 and therefore an overestimation of λ in such a regime. Even if the implemented FLD prescription is still too diffusive at τ < 1, it is the best option to describe radiative transport there within MESA’s limitations.

Appendix B: Limiting convective velocities

With the motivation of avoiding supersonic convection transport, the capability to limit convection velocities has been re-implemented in MESA⁵. For this work, we choose the mixing length theory (MLT) of convection by Mihalas (1978) and Kurucz (1979) with α_MLT = 1.93 (Choi et al. 2016; Fragos et al. 2023). Convective velocities are limited in terms of the local sound speed (c_sound) to be not grater than 0.1c_sound. The considered MLT prescription is only valid for non-supersonic convection velocities where the turbulent pressure is negligible with respect to the gas pressure (see Weiss et al. 2004, for more context on subsonic convection).

Appendix C: Extrapolation of the opacity tables

In order to describe the interaction between radiation and matter, MESA’s opacity module (kap) contain several options for gray opacity tables that are blended in the ℛ-T plane, where T is the local stellar temperature and ℛ ≡ (ρ/g cm⁻³)/(T/10⁶ K)³. In a simple manner, there are four sets of the opacity data that are blended:

• the low temperature tables (usually to describe κ from the temperature of formation of molecules and dust grains, (T ∼ 500 K) until a T value below the hydrogen recombination peak, T ∼ 6000 K),

• the high temperature tables (for κ from the hydrogen recombination T to T ∼ 10⁸ K),

• the Compton scattering prescription by Poutanen (2017, to compute κ for T > 10⁸ K), and

• the thermal conductive opacity based on the study by Cassisi et al. (2007).

During the development of this work we noticed that MESA was blending the Compton scattering opacity with the lower end of the low temperature tables, leading to an underestimation of κ of almost 30 orders of magnitude for regions with T < 3000 K or ℛ < 10^−7.56. Given the previous nonphysical behavior, we fixed the blending of the Compton opacity to strictly be used only if T > 10⁸ K or ℛ > 10, that correspond the T and ℛ limits of the tables used respectively.

For this study we consider the low temperature tables by Ferguson et al. (2005) and the high temperature tables by Iglesias & Rogers (1993, 1996). Both tables are only defined for 0.1 < ℛ < 10 and T > 500 K. To avoid not having a defined opacity value for the low-temperature and low-density EMS surface as it expands, we extrapolate only the low and high temperature tables by considering the nearest-neighbor value, in the ℛ-T plane, for the regions where ℛ < 0.1 or T < 500K.

The resulting opacity table used for the stellar models in this work is presented in Fig. A.1. The default MESA opacity tables are extrapolated below their minimum values for temperature and ℛ. The EMS models in this work cross these lower boundary values of the default tables only once the stellar surface expands due to the energy injected by the stellar merger.

Fig. A.1. Values for the stellar opacities in the log(T) vs. log(ρ) plane for Z = 10−4. The color map denotes the value of the logarithm of the opacity in [cm2 g−1]. The limits for the domains of each opacity table/prescription are shown in white. The solid lines represent the boundaries in temperature, while the dashed lines represent the boundaries in ℛ values. The intermediate regions between the low T and high T tables, as wells as between the high T and the Compton opacities, are the blending zones between such domains. The low and high T tables are being extrapolated when ℛ < 10−1 or T < 500K by considering the nearest neighbor value in terms of the temperature and ℛ.

Appendix D: Companion stars: Mass, radius, and virial temperature

The companion star models were built with MESA trough cold gas accretion following directly the prescription by Haemmerlé et al. (2019), in the same manner as in Ramírez-Galeano et al. (2025). We extract the mass radius relationship from those stellar models as it is shown in Fig. D.1. As well, the companion’s virial temperature is computed as

$$\begin{array}{c}\hfill T_{\mathrm{virial},\mathrm{comp}}=\frac{G\mu m_{\mathrm{H}}{m}_{\mathrm{comp}}}{k_{\mathrm{B}}{R}_{\mathrm{comp}}},\end{array}$$(D.1)

where μ is the companion’s mean molecular weight, m_H is the mass of the hydrogen atom and k_B is the Boltzmann constant. The values for m_comp, R_comp and T_{virial, comp} for this work are shown in Table D.1.

Fig. D.1. Mass–radius relation for companion stars used in the collision simulations. The stellar radii are plotted as a function of companion mass, this relation is used to to determine the companion’s minimum cross-section as well as its virial temperature.

Appendix E: On the remaining fraction of the deposited energy at the end of the inspiral

To estimate the amount of deposited orbital energy remaining at merger, the values of the luminosity change fraction f_L are calculated following Eq. 8. The resulting values for this parameter are shown in Fig. E.1 for all EMS mass and companion mass configurations.

Then, the fraction of orbital energy remaining at merger with respect to the total injected energy, α(t_merger), as defined in Eq. 12, is shown in Fig. E.2. The solid-dotted lines correspond to the values of α obtained from our simulations, while the dashed lines represent the analytical estimates from Ramírez-Galeano et al. (2025), based on their assumptions for the fraction of deposited orbital energy available to drive mass loss. The different colors indicate the three EMS masses considered in this work. A direct comparison shows that the analytical prescription of Ramírez-Galeano et al. (2025) systematically underestimates α, and therefore underestimates the amount of deposited orbital energy available to drive mass loss.

Our simulations also reveal that collisions with low-mass companions yield significantly smaller values of α. By rewriting Eq. 12 as in Eq. 13, we identify two main factors responsible for this trend:

• the smaller gravitational cross-section of low-mass companions leads to considerably longer merger timescales, t_merger, compared to those for high-mass companions;

• due to their lower mass, the total orbital energy deposited by such companions is expected to be smaller than that injected by high-mass ones, which in turn requires less time for the EMS to radiate away this energy.

These two effects combined reduce the fraction of deposited orbital energy remaining at merger for low-mass companions.

Fig. E.1. Fractional change in luminosity, fL, averaged from the inspiral start until the merger condition is reached, as function of companion mass. The different colors correspond to the three EMS masses adopted in this study, i.e., 1000, 2000 and 3000 M⊙.

Fig. E.2. Fraction of the deposited orbital energy remaining at merger, α(tmerger), as a function of companion mass. Solid-dotted lines show simulation results, while dashed lines show the analytical estimates given the assumptions by Ramírez-Galeano et al. (2025). Colors indicate the three EMS masses considered, i.e., 1000, 2000 and 3000 M⊙.

Appendix F: Mass-loss in absolute values

As in Fig. F.1, we show in Fig. 6 the maximum mass-loss estimate and the latus rectum threshold for stellar mergers in absolute values, instead of the relative ones to the EMS mass and radius respectively. For further context for this figure, see Sect. 4.3.

Fig. F.1. Absolute estimates of merger-induced mass loss and the latus rectum threshold between stellar mergers and scatterings, for orbits with initial eccentricity of 1.01. Considering EMS masses from 103 to 104 M⊙ and companion masses from 5 to 500 M⊙. Top plot: maximum merger-induced mass loss in M⊙. Bottom plot: maximum initial latus rectum for a stellar merger in R⊙.

All Tables

All Figures

Fig. 1. Inspiral evolution of a 1000 M⊙ EMS with a 10 M⊙ companion. Panel a: Structural and dynamical evolution of the 1000 M⊙ EMS during inspiral with a 10 M⊙ companion. Top: EMS Kippenhahn diagram, in terms of the radius coordinate with respect to time, with the signed logarithm of the binding energy encoded by the background color; Bottom: same but for log(v/vesc). In both panels, dotted cyan regions mark convection zones. The solid black line represents the position of the companion, the black dot marks the merger event. Panel b: Zoom-in on the merger event. Top: EMS radius with v/vesc color map, the value of the radius at t = 0 is shown by the dashed red line; Middle: Unbound (Munb) and escaping (Mesc) masses. Bottom: Surface pressure, its initial value is shown by the dashed red line. In all three plots the time when the merger event is met is marked by the black vertical line.

In the text

Fig. 2. Inspiral evolution of a 1000 M⊙ EMS with a 70 M⊙ companion until the simulation crashes. Panel a: structural and dynamical evolution of the 1000 M⊙ EMS during inspiral with a 70 M⊙ companion. Top: EMS Kippenhahn diagram, in terms of the radius coordinate with respect to time, with the signed logarithm of the binding energy encoded by the background color; Bottom: Same but for log(v/vesc). In both panels dotted cyan regions marks convection zones. The solid black line represents the position of the companion, the black dot marks the merger event. Panel b: Zoom-in on the merger event. Top: EMS radius with v/vesc color map, the initial radius value is shown by the dashed red line; Middle: Unbound (Munb) and escaping (Mesc) masses. Bottom: Surface pressure, its initial value is shown by the dashed red line. In all three plots the time when the merger event is met is marked by the black vertical line.

In the text

Fig. 3. Maximum unbound mass (Munb, max, purple triangles pointing down), maximum escaping mass (Mesc, max, mint triangles pointing up), and predicted total pulsation mass loss (Mpul, tot, turquoise diamonds) as a function of companion mass for the three EMS models considered. The dashed gray line indicates the one-to-one relation, a predicted mass loss below this line means that the EMS gains more mass by the stellar merger than what it loses as a result of it and vice versa. The dotted blue line shows the maximum unbound mass estimation if there is no loss of injected energy during the inspiral, namely, α = 1 at all times. The dashed black line shows the extrapolated mass-loss prescription from Glebbeek et al. (2013). Analytical predictions from Ramírez-Galeano et al. (2025) are shown with the solid cherry-hue lines.

In the text

Fig. 4. Collision outcomes and deposited orbital energy for the interaction between a 1000 M⊙ EMS and a 70 M⊙ companion across initial orbital configurations. The EMS structural profile after the relaxation steps is used for trajectory integration. The initial semilatus rectum ranges from 10 to 105 R⊙ and the eccentricity from 0 to 2. The EMS radius is indicated by the dashed gray line. Left plot: Classification of collision outcomes (A–D; see Sect. 4.3) as a function of the initial semilatus rectum and eccentricity. Right plot: Deposited orbital energy for mergers and scatterings in the absence of external interactions.

In the text

Fig. 5. Same as Fig. 4, but using the EMS profile at the end of the hydrodynamic collision simulation for the interaction between a 1000 M⊙ EMS and a 70 M⊙ to determine the possible outcomes of a second collision. The EMS has expanded to a radius of ∼5 × 104 R⊙, shown by the dashed gray line in both panels. Left plot: Collision outcome classification, right plot: Deposited orbital energy for the same configurations.

In the text

Fig. 6. Relative estimates of merger-induced mass loss and the semilatus rectum threshold between stellar mergers and scatterings, for orbits with initial eccentricity of 1.01. Considering EMS masses from 103 to 104 M⊙ and companion masses from 5 to 500 M⊙. Left plot: Ratio of the maximum merger-induced mass loss to the EMS mass. Right plot: Ratio of the maximum initial semilatus rectum for a stellar merger to the EMS radius.

In the text

Fig. A.1. Values for the stellar opacities in the log(T) vs. log(ρ) plane for Z = 10−4. The color map denotes the value of the logarithm of the opacity in [cm2 g−1]. The limits for the domains of each opacity table/prescription are shown in white. The solid lines represent the boundaries in temperature, while the dashed lines represent the boundaries in ℛ values. The intermediate regions between the low T and high T tables, as wells as between the high T and the Compton opacities, are the blending zones between such domains. The low and high T tables are being extrapolated when ℛ < 10−1 or T < 500K by considering the nearest neighbor value in terms of the temperature and ℛ.

In the text

	Fig. D.1. Mass–radius relation for companion stars used in the collision simulations. The stellar radii are plotted as a function of companion mass, this relation is used to to determine the companion’s minimum cross-section as well as its virial temperature.
In the text

	Fig. E.1. Fractional change in luminosity, fL, averaged from the inspiral start until the merger condition is reached, as function of companion mass. The different colors correspond to the three EMS masses adopted in this study, i.e., 1000, 2000 and 3000 M⊙.
In the text

	Fig. E.2. Fraction of the deposited orbital energy remaining at merger, α(tmerger), as a function of companion mass. Solid-dotted lines show simulation results, while dashed lines show the analytical estimates given the assumptions by Ramírez-Galeano et al. (2025). Colors indicate the three EMS masses considered, i.e., 1000, 2000 and 3000 M⊙.
In the text

	Fig. F.1. Absolute estimates of merger-induced mass loss and the latus rectum threshold between stellar mergers and scatterings, for orbits with initial eccentricity of 1.01. Considering EMS masses from 103 to 104 M⊙ and companion masses from 5 to 500 M⊙. Top plot: maximum merger-induced mass loss in M⊙. Bottom plot: maximum initial latus rectum for a stellar merger in R⊙.
In the text