© The Authors 2025

1. Introduction

Stellar collisions can occur in dense stellar environments such as globular clusters or nuclear stellar clusters, as shown in Fig. 1 depicting the parameter ranges of the relative velocity between stars and stellar density for collisions¹. As indicated by the magenta curve representing the critical density at which the collision timescale equals the lifetime of 1 M_⊙ main-sequence stars², stars can undergo multiple collisions in both globular and nuclear stellar clusters, including the Galactic nuclear cluster corresponding to the line for 4×10⁶ M_⊙ black hole.

Fig. 1. Ranges of the relative velocity between stars and stellar density for destructive (magenta-shaded) and non-destructive (gray-shaded) collisions, separated by the surface escape velocity of a 1 M⊙ main-sequence star. The velocity v and distance R from the black hole are related by $v/c=\sqrt{r_{\mathrm{g}}/R}$, with the black-shaded region (“Inside black hole”) indicating v/c>1. Here, rg is the gravitational radius. Black diagonal lines indicate nuclear clusters with the same mass as the black hole, assuming a Bahcall & Wolf (1976) profile extending inward (arrow) from the influence radius until 2 M⊙ is enclosed or the event horizon is reached. The red and blue crosses indicate the observationally constrained density at 0.25 pc from the Galactic center (Genzel et al. 2010) and the semimajor axis of the S2 star (Gillessen et al. 2009), respectively.

Collisions can be divided into two regimes based on the collision velocity relative to the escape velocity of two colliding stars (gray vertical line in Fig. 1)³. When the relative velocity between two colliding stars is low enough that the total energy is negative (gray-shaded region), the two stars would merge into a more massive star, which may appear as a blue straggler (e.g., Sandage 1953; Burbidge & Sandage 1958; Sills et al. 1995; Wang & Ryu 2024; Ryu et al. 2024a). These low-velocity collisions typically occur in globular clusters (cyan region in the lower-left corner, Baumgardt & Hilker 2018), where the stellar velocity dispersion is ≲15 km s⁻¹ (Cohen 1983; Pryor & Meylan 1993; Harris 2010). In contrast, in nuclear stellar clusters surrounding a supermassive black hole (diagonal black lines), the relative velocities between stars near the black hole can reach ≳O(10³) km s⁻¹. When two stars collide at such high velocities, exceeding their stellar binding energy, they can be fully or partially destroyed (magenta-shaed region) (e.g., Freitag & Benz 2005). This process, which occurs at a rate of 10⁻³−10⁻⁷ yr⁻¹ galaxy⁻¹ (Rose et al. 2020, 2023; Amaro Seoane 2023a; Balberg & Yassur 2023; Hu & Loeb 2024), can produce bright electromagnetic transients (Balberg et al. 2013; Balberg & Yassur 2023; Amaro Seoane 2023a, b; Ryu et al. 2024b; Dessart et al. 2024), making them promising candidates for nuclear transients for current and future surveys.

What do their light curves and spectra look like? This is one of the central questions to be addressed in identifying transients caused by destructive stellar collisions. Previous analytic studies (e.g., Balberg et al. 2013; Amaro Seoane 2023a) estimated transient luminosities and durations, assuming photons diffuse out in spherically expanding ejecta. More specifically, Balberg et al. (2013) made an analytic estimate of the luminosity for the first time, considering a scenario where the emitted energy from the region where the diffusion time equals the expansion time contributes to the luminosity, as done in supernova light curve analysis. On the other hand, Amaro Seoane (2023a) attempted to refine the estimate by introducing a factor for the conversion efficiency from kinetic energy to radiation energy and by speculating on the time evolution of the diffusion time. While these estimates already suggested potentially luminous transients, they remained approximate. Recently, Dessart et al. (2024) performed detailed 1D radiative-transfer calculations of light curves and spectra for high-velocity collisions between giant stars – these calculations were themselves based on the results from hydrodynamics simulations performed by Ryu et al. (2024b). The radiative-transfer simulations suggest that the luminosities of giant collisions are on the order of 10⁴²−10⁴⁴ erg s⁻¹ at 1 d after the onset of the collision, which is in the ballpark of tidal disruption events (e.g., Gezari 2021; Wevers & Ryu 2023) or supernovae. In addition, the collision debris reveal a similar spectral evolution as that of Type II supernovae from blue-supergiant star explosions (Dessart et al. 2024).

This paper investigates the observables of high-velocity collisions between realistic main-sequence stars at their terminal age, exploring a wide range of collision configurations, including stellar masses (M₁ = 1, 3, 10 M_⊙ and M₂ = 0.9M₁), mass ratios (q=M₂/M₁≃0.1−0.9), relative collisional velocities (2500−10 000 km s⁻¹), and impact parameters ranging from nearly head-on to off-axis collisions. Employing the same methodology as in Dessart et al. (2024), we first performed hydrodynamics simulations of collisions using the moving-mesh hydrodynamics simulation code AREPO (Springel 2010; Weinberger et al. 2020; Pakmor et al. 2016). In most cases, these eventually produce quasi-spherically expanding collision debris. We then conducted 1D, nonlocal thermodynamic equilibrium, time-dependent radiative-transfer calculations with the code CMFGEN (Hillier & Dessart 2012) using 1D, angle-averaged profiles built from the 3D hydrodynamics simulations at 0.5 d.

This paper is organized as follows. In Sect. 2, we describe our hydrodynamics simulations (Sect. 2.1), followed by a dynamical overview of stellar collisions (Sect. 2.2). In Sect. 3, we provide a description for the setup of our radiative transfer calculations (Sect. 3.1) and an overview of the photometric and spectroscopic properties of representative models (Sect. 3.2). Using the full set of nearly forty simulations, we discuss correlations between collision parameters and the modeled observables in Sect. 4. We present our conclusions in Sect. 5.

2. Hydrodynamics simulations

2.1. Numerical setup

Following Ryu et al. (2024b), we performed a suite of 3D hydrodynamics simulations of high-velocity collisions between main-sequence stars using the moving-mesh magnetohydrodynamics code AREPO (Springel 2010; Pakmor et al. 2016; Weinberger et al. 2020). We employed an equation of state for a mixture of ideal gas and radiation assuming local thermodynamic equilibrium, in which the total pressure is the sum of the radiation pressure and gas pressure, i.e., P=aT⁴/3+ρk_BT/[μm_p], where a is the radiation constant, T the temperature, ρ the density, k_B the Boltzmann constant, μ = 0.63 the mean molecular weight assuming a fully ionized plasma with Solar metallicity, defined in terms of the proton mass m_p. This assumption is valid because the collision product is highly optically thick until the end of simulations. For the initial conditions used in AREPO, we used stellar models evolved at solar metallicity (Z = 0.02) with MESA (version r24.03.1; Paxton et al. 2011, 2013; Jermyn et al. 2023). We adopted the prescriptions for overshoot, convection, semiconvection, and thermohaline mixing from Choi et al. (2016). We used the Ledoux criterion Ledoux (1947) to identify the boundary of the convective regions. We considered main-sequence stars with six different masses, ranging from 0.9 to 10 M_⊙ initially, and evolved them from their original uniform composition until the central hydrogen mass fraction reached 0.01, representing the terminal age of the main sequence. This choice is motivated by the fact that the stellar radius increases toward the terminal age, which will increase the collision probability⁴. At that time, the surface radius of our six models is 1.1 R_⊙ for 0.9 M_⊙ star, 1.2 R_⊙ for 1 M_⊙ star, 4.6 R_⊙ for 2.7 M_⊙ star, 4.9 R_⊙ for 3 M_⊙ star, 9.3 R_⊙ for 9 M_⊙ star, and 10 R_⊙ for 10 M_⊙ star. We show the density and the hydrogen mass fraction versus Lagrangian mass of each stellar model in Fig. 2.

Fig. 2. Internal structure of stellar models. Density (top) and H mass fraction (bottom) profiles of our main-sequence stellar models with a metallicity of Z = 0.02 and at the core H mass fraction of 0.01, as a function of enclosed mass (left) and radius (right).

Using the mapping routine of Ohlmann et al. (2017), we relaxed each 1D stellar model in 3D using N = 10⁶ cells. Then, we conducted collision experiments by placing a pair of such relaxed main-sequence stars 8 (R₁+R₂) apart and with a relative velocity V_rel (R₁ and R₂ are the surface radii of the colliding stars). Using this numerical setup, we simulated stellar collisions in which were varied the stellar masses M₁ and M₂, the relative velocity, and the impact parameter b – this parameter b will be quoted in units of (R₁+R₂). We considered three suites of (nearly) equal-mass collisions with q=M₂/M₁ = 0.9 (i.e., modelsets # 1, 2, and 3), as well as collisions with widely different mass ratios of 3 to 10 (i.e., modelsets # 4, 5, and 6). For each modelset, we considered initial configurations with b = 0.1, 0.25, and 0.5 for V_rel = 5000 km s⁻¹, and b = 0.1 and 0.5 for V_rel = 2500 km s⁻¹ and 10⁴ km s⁻¹. The complete list of models simulated with AREPO is given in Table A.1.. The model nomenclature gives the mass of each stellar component, the relative velocity in units of 1000 km s⁻¹ and the impact parameter b in units of (R₁+R₂). For example, model ms1/ms0p9/vrel2p5/b0p1 corresponds to a collision with parameters M₁ = 1 M_⊙, M₂ = 0.9 M_⊙, V_rel = 2500 km s⁻¹, and b = 0.1 (R₁+R₂).

We followed the evolution with AREPO up to 0.5 d after the collision. In all simulations, the fractional errors in total energy remain small (≲1%) until the end of the simulations. To ensure the robustness of our simulations, we conducted resolution tests for an equal-mass collision involving a 10 M_⊙ star with V_rel = 5000 km s⁻¹ and b = 0.1, varying N per star from 1.25 to 20×10⁵. We found that the overall evolution and key quantities, such as the conversion efficiency between kinetic energy and radiation energy or the structure of the debris, are essentially unchanged when N≳5×10⁵, which is two times lower than the resolution adopted for the entire grid.

2.2. Dynamical overview

We show the evolution of the density for collisions with different impact parameter b, first for three models from the set # 1 (M₁ = 1 M_⊙ and M₂ = 0.9 M_⊙) in Fig. 3 and and then for three models from the set # 5 (M₁ = 10 M_⊙ and M₂ = 1 M_⊙) in Fig. 4. The overall evolution is qualitatively similar to collisions between giants studied in Ryu et al. (2024b). At the onset of the collision, strong shocks form along the contact surface, as shown in the first columns of both figures, converting 0.1%−50% of the total kinetic energy into thermal energy, depending on the collision parameters. We note that the conversion efficiency for collisions with b = 0.5 is relatively low (0.1−4%) and independent of the stellar masses and V_rel. Although the shocks are weak in some cases, they lead to the formation of homologously, supersonically expanding debris after a few hours. While small structures exist within such ejecta, the overall shape is approximately spherical, which is consistent with the assumption for the debris made in Balberg et al. (2013) and Amaro Seoane (2023a).

Fig. 3. Gas density distribution at successive moments in main-sequence star collisions with M1 = 1 M⊙, M2 = 0.9 M⊙, Vrel = 5000 km s−1, but different impact parameters. These AREPO simulations of stellar collisions are characterized by b = 0.1, 0.25, and 0.5 (from top to bottom) and shown here at t = 0.0, 0.01, 0.05, and 0.5 d (from left to right).

Fig. 4. Same as Fig. 3, but for collisions with M1 = 10 M⊙ and M2 = 1 M⊙, Vrel = 5000 km s−1. In all cases, the collision unbinds only a fraction of the envelopes of both stars, leaving behind remnants.

The shocks convert kinetic energy into radiation energy. Figure 5 shows the radiation energy E_rad,ej and kinetic energy E_kin,ej within the ejecta in the two models considered in Figs. 3 (left) and 4 (right), as a function of time. At collision (t = 0), E_rad,ej peaks, with values ranging from 10⁴⁷−4×10⁴⁹ erg for modelset #1 and from 6×10⁴⁸−10⁵⁰ erg for modelset #5, depending on the collision parameters. It then drops by roughly two to three orders of magnitude in 0.5 d, which is the start time of our radiative-transfer calculations (see Sect. 3). Our hydrodynamics simulations do not account for radiative energy loss– they assume adiabatic conditions. While this approximation should hold well for the optically thick regions of the debris, the radiation energy in the outer regions – where the optical depth approaches unity – may be overestimated. However, this should not introduce significant errors, as the excess heat is radiated away in the first time step of the radiative transfer calculations. E_kin,ej ranges from 10⁴⁸−4×10⁵⁰ erg for both modelsets. In highly energetic collisions with small impact parameters, E_kin,ej remains constant, whereas in collisions with relatively low V_rel and large b, E_kin,ej initially tends to decrease gradually and then levels off within 0.5 d after the collision.

Fig. 5. Evolution of the ejecta kinetic and radiative energies in a sample of stellar-collision simulations. We show radiation energy Erad,ej (top) and kinetic energy Ekin,ej (bottom) for modelsets #1 (left, 1 M⊙+0.9 M⊙) and #5 (right, 10 M⊙+1 M⊙). The legend indicates the initial relative velocity (e.g., vrel5 means Vrel = 5000 km s−1) and impact parameter (e.g., b0p5 means b = 0.5). The diagonal black lines in the top panels indicate the power law of t−1, expected for an adiabatic expansion of a radiation-dominated plasma.

Even if the collision kinetic energy exceeds the binding energy of both stars, the collision does not generally destroy them, but instead leaves behind stellar remnants, composed of the more compact, denser He-rich core of the star prior to collision. The only cases where both stars are completely destroyed are when M₁ = 1 M_⊙, M₂ = 0.9 M_⊙, V_rel≳5000 km s⁻¹, and b = 0.1.

The results for the full sample of simulations is stored in Table A.1. and indicate debris or ejecta masses M_ej in the range from a few 0.01 M_⊙ to slightly over 10 M_⊙, ejecta kinetic energies in the range from a few 10⁴⁷ erg to several 10⁵¹ erg, and stored radiative energies at 0.5 d in the range from a few 10⁴⁷ erg to about 6×10⁴⁸ erg.

3. Radiative-transfer simulations

3.1. Numerical setup

Our main goal is to obtain accurate light curves and spectra of high-velocity collisions. As in Ryu et al. (2024b), the bolometric luminosities from the debris can be estimated based on the radiation energy and photon diffusion time derived from the 3D hydrodynamic simulations. However, such an estimate is crude. Moreover, spectral calculations are not possible solely using the hydrodynamic properties of the debris. Therefore, once the debris begin to expand in a quasi-spherical and homologous fashion, which occurs at ≃0.5 d, we stopped the 3D hydrodynamic simulations and performed radiative transfer calculations using the properties of the debris obtained from the 3D simulations as initial conditions to calculate light curves and spectra. The radiative-transfer calculations performed here on main-sequence star collisions follow the same procedure as described in Dessart et al. (2024), where collisions between giant stars were studied. Hence, we refer the reader to that work and in particular their Section 2 for a detailed presentation of the numerical approach. Here, we merely provide a concise summary.

We performed 1D nonlocal thermodynamic equilibrium (NLTE) time-dependent radiative-transfer calculations with the code CMFGEN (Hillier & Dessart 2012) using the AREPO simulation results for a whole range of stellar collisions. Since this radiative-transfer code does not treat hydrodynamics and assumes spherical symmetry, we started the simulations at a sufficiently late time after the onset of the collision, that is once the debris of the collision were coasting. Because of the relatively compact structure of the main-sequence stars used (see Fig. 2), this occurs after a few hours. Furthermore, we built a spherical average of the fluid quantities obtained in the AREPO simulation. The spherical average was calculated so that the ejecta mass in each radius bin was the same in 1D and 3D profiles. However, the kinetic energy and radiation energies (or radial velocity and temperature) were not strictly conserved. However, as shown in Table A.1., the difference between the two energy estimates is not significant, less than a factor of two in most cases. It becomes significant only for b = 0.5 and V_rel = 2500 km s⁻¹ (a factor of 5−10), indicating our assumption of spherical symmetry for expanding ejecta is not valid in such cases. This also implies that viewing angle effects are important for these collisions.

To cover the early evolution of the radiative properties of such collisions when they are the most luminous, we started our simulations at a post-collision time of 0.5 d. Unlike Dessart et al. (2024) in which the evolution was followed until the debris became optically thin (i.e., when the total radial electron-scattering optical depth drops below one), we evolved most of our simulations until the debris entered the H recombination phase, at a time when the luminosity from the debris has dropped by a factor of ten to a hundred. By that time, the associated transient would likely be undetectable at representative distances of 10−100 Mpc. The other reason for starting early is that we are interested in documenting the properties of our simulations when they are UV bright, which is also when they are the most luminous, and thus exhibit favorable properties for UV observations with future satellites like UVEX (Kulkarni et al. 2021) or ULTRASAT (Shvartzvald et al. 2024) – these space observatories will routinely detect and monitor transients merely hours after explosion, collision, or disruption, depending on the type of transient. With this aim, we could make another adjustment in setting up the initial models by truncating the inner, dense, and slow moving regions of the debris. This adjustment should not affect our main results for light curves and spectra over the duration considered (≲10 d). Indeed, at 0.5 d, these regions tend to be exceptionally dense and slow, where an optical depth of ≳10⁵, and containing merely a few percent of the total stored radiative energy at that time (see penultimate column in Table A.1.). This means that because of relatively long photon diffusion times, these regions would influence escaping radiation at later times – days and weeks – only marginally. A further complication is that at such later times, accretion of debris material into the supermassive black hole can lead to a power injection into the expanding debris (see Ryu et al. 2024b; Brutman et al. 2024) and alter its energy budget – this influence should be small or even negligible at earlier times when the debris are totally optically thick but would completely dominate at late times. The longer-time evolution on week-timescales of our models would require a more global approach to account for both the debris and the supermassive black hole, which we delay to a forthcoming study.

Collisions with small impact parameters and high stellar masses reached debris masses and kinetic energy comparable to what is found in Type II supernovae (i.e., about 10 M_⊙ and 10⁵¹ erg) but with a very small radiative energy content on the order of 10⁴⁸ erg more typical of (compact) blue-supergiant (Dessart & Hillier 2019) rather than (extended) red-supergiant star explosions (see, for example, the simulations of Dessart et al. 2013; Dessart & Hillier 2019). These quantities are stored in Table A.1.. In this work, we performed radiative-transfer calculations for 35 collisions but excluded the weakest collisions (some in the sample having small relative velocities and large impact parameters) with the lowest radiative-energy content since their transients would be too faint to be observed. In a few cases (e.g., model ms10/ms9/vrel2p5/b0p5), the model was not included because of numerical difficulties when evolving the ejecta with CMFGEN, in particular when the initial temperature structure was highly nonmonotonic (e.g., a strong temperature peak at intermediate velocities in the debris but low temperatures in the inner, dense debris regions). The suite of radiation transfer simulations is sufficiently large to establish correlations between collision parameters and observables.

For the CMFGEN calculations, we assume a solar-metallicity composition for all metals above O but take the composition profile from the AREPO simulations for species H, He, C, N, and O (the N/C abundance ratio is increased by a factor of a few in the H-depleted He-rich core of the original stellar models; see also Fig. 2). Specifically, we treat H, He, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca, Sc, Ti, Cr, Fe, Co, Ni, and Ba. Following Dessart et al. (2024), we include the following ions for the above species: H I, He I–II, C I–IV, N I–V, O I–V, Ne II–IV, Na I, Mg II–IV, Al II–III, Si I–IV, S II–IV, Ar I–III, Ca II–IV, Sc II–III, Ti II–III, Cr II–V, Fe I–V, Co II–V, Ni I–V, and Ba II. All associated isotopes being stable, there is no contribution from radioactive decay and no nonthermal effects associated with such decays. All spherical-averaged debris inputs are remapped onto a grid with 80 points uniformly spaced on a logarithmic optical-depth scale. Remapping the grid is also necessary at the start of every time step as well as during the calculation in order to track the H I and He I–He II recombination fronts that form as the debris cool through expansion and radiation escape. To cover the full evolution, a total of about 20−30 timesteps (depending on the timespan covered) were computed for most models. In some models, the calculation was stopped when the debris were becoming very faint (e.g., luminosity below 10⁶ L_⊙) during the recombination phase. The main CMFGEN output of interest for this study is the emergent flux, which is calculated at each time step in the observer's frame from the far-UV to the far-IR. This spectrum can be shown directly or used to compute the bolometric luminosity L_bol or the absolute magnitude in various filters. We selected the Swift UVW1-band and the optical V-band filters. The UVW1 filter is a close analog of the NUV filter that will be on board ULTRASAT (Shvartzvald et al. 2024).

3.2. Overview of photometric and spectroscopic evolution

The radiative properties of main-sequence star collisions are found to be qualitatively similar to those of giant-star collisions and discussed in detail in Dessart et al. (2024). We thus summarize the properties obtained for modelset # 1 as a representative case and compare the results obtained for different collision configurations. We summarize the results of our radiative transfer calculations in Table 1 and provide some additional results and illustrations for the full sample in Appendix B. The main novel aspect of this work is the extent of the model grid and the correlations drawn from it (see Sect. 4).

Figure 6 illustrates the photometric properties of models in the set # 1 characterized by M₁ = 1 M_⊙ and M₂ = 0.9 M_⊙. All models exhibit the same bolometric light curve morphology, with values in the range 10^41.5 to 10^42.5 erg s⁻¹ at 0.5 d and dropping rapidly (L_bol∝t⁻²−t⁻⁴ at 0.5−1 d) to a plateau at 10⁴⁰ to 10⁴¹ erg s⁻¹ at about 1.5 d. Given the fading rate at 0.5 d, the peak luminosities probably reach up to 10⁴²−10⁴³ erg s⁻¹. The small radiative energy content (roughly 10⁴⁶ to 10⁴⁷ erg) and mass (∼0.05 to ∼1 M_⊙) limit the luminous, optically thick phase to about a week (this was not calculated for all models in the set). All these models are less luminous than SN 1987A at the same epochs. The UV brightness (shown here with the Swift UVW1 filter) follows a similar evolution as the bolometric luminosity. The models are brightest in the UV at the earliest times. As a result of expansion and radiation cooling, they rapidly cool and redden and fade in the UV. For all models, the fading rate in the UV is large and roughly in the range 4−6 mag per day at one day. In contrast, the concomitant drop in luminosity and ejecta temperature leads to a quasi-plateau V-band light curve, with values in the range of −11 to −14 mag a few days after the onset of the collision. Here again, all these models are optically fainter than SN 1987A at the same epochs. Overall, the contrast between models follows directly from the initial conditions in the sense that the models with the greater stored radiative energy at 0.5 d are brighter (in any filter) or bolometrically more luminous (all models have the same color evolution). Lacking any long-lived power source, the present models become exceedingly faint at late times. Accretion onto the supermassive black hole would dramatically change this picture.

Fig. 6. Photometric properties for main-sequence star collisions from the modelset # 1. From top to bottom, we show the bolometric, UVW1-band, and V-band light curves. Filled dots indicate the corresponding quantities inferred for SN 1987A (after correction for reddening and distance). Corresponding results for all sets are shown in Appendix B.

Figure 7 illustrates the spectral evolution from UV to optical for model ms1/ms0p9/vrel10/b0p1 (also from set # 1) – the evolution for other models in any modelset is qualitatively similar and also analogous to results obtained previously by Dessart et al. (2024). It reveals in more details what is conveyed by the photometric evolution, with the rapid shift of the spectral energy distribution from the UV (brighter than the optical until about 2 d) to the optical. Initially strong and isolated lines from low-lying levels (e.g., He II λ 1640, C III λ 2297, C IV λλ 1548, 1551, or Mg II λλ 2795, 2802) coexist with extended spectral regions blanketed by a multitude of metal lines (in particular those due to Fe II–IV). By two days after the onset of the collision, most of the flux falls in the optical and exhibits the whole series of H I Balmer lines, with Hα the strongest and broadest of all lines. Strong Fe II lines at and below 5000 Å are present, together with Na I λλ 5896, 5890 or the Ca II NIR triplet near 8500 Å. At ∼6 d, the rapid drop in luminosity shown in Fig. 6 corresponds spectroscopically to the weakening of the continuum flux and the transition to the nebular phase. As the debris release the small residue of stored radiative energy, the spectrum exhibits a weakening continuum flux and strong, broad, and quasi-symmetric emission lines. This occurs because of the low optical depth of the debris. This apparent strength of lines is only relative since all fluxes are dropping (i.e., the bolometric flux drops but the flux in the lines drops slower than that in the continuum during that transition).

Fig. 7. Spectral evolution from the UV to the optical and from ∼0.5 until ∼7.2 d for model ms1/ms0p9/vrel10/b0p1 (modelset # 1). Fλ represents the spectral flux density or the flux per unit wavelength interval. Corresponding results for all models from the six sets are shown in Appendix B.

Figure 8 illustrates the full diversity of spectra from the UV to the optical at one day, including all models from all sets (a few models were not computed – see the list in Table A.1.) but grouped according to the value of the collision velocity V_rel and impact parameter b. To allow the simultaneous assessment of the UV and the optical, we show the quantity λ²F_λ. In essentially all models, the spectral energy distribution peaks around 2000 Å, although models with a greater V_rel and smaller b tend to be brighter in the UV. Most models exhibit broad lines at one day, but these Doppler-broadened lines tend to be narrower in configurations involving smaller collision velocities or higher stellar masses. The spectral properties obtained here are very similar to those of SN 1987A or Type II supernovae in general (see discussion in Dessart et al. 2024). Therefore, detailed studies of SNe II can be used as a source of information on what is displayed here (e.g., Dessart & Hillier 2011; Dessart et al. 2013).

Fig. 8. Spectral comparisons of models grouped by the same values of Vrel and b but differing in the stellar masses of the colliding pair. The time is one day after the onset of the collision. To better reveal the weaker flux in the optical, we show the quantity λ2Fλ, which is analogous to an Fν.

4. Correlations

In this section, we examine possible correlations between ejecta properties and collision parameters, which can be used for building semi-analytic models for destructive collisions, as well as those between observables and collision parameters, which can be used to infer collision parameters from observations. We used the ejecta properties at 0.5 d and stored in Table A.1. (see Sects. 4.1 and 4.2), and the radiative-transfer results at 1 d and stored in Table 1 (see Sect. 4.3). We also confronted these to their counterparts from simulations of collisions between giant stars (GS) (measured at 1 d) in Ryu et al. (2024b) and Dessart et al. (2024) – these works considered collisions between two identical 1 M_⊙ giants with different radii and impact parameters of b = 0.02−0.4⁵.

4.1. Ejecta mass

The ejecta mass M_ej (see Table A.1.) primarily depends on the impact parameter b, with larger M_ej resulting from more head-on collisions. It also exhibits a weak positive correlation with the initial kinetic energy E_kin. This trend is illustrated in Fig. 9. As shown in the top panel, more head-on collisions produce larger M_ej. M_ej falls into three distinct ranges depending on b: M_ej≳1 M_⊙ for b = 0.1, 0.1 M_⊙≲M_ej≲1 M_⊙ for b = 0.25, and M_ej≲0.1 M_⊙ for b = 0.5. However, the distribution of M_ej is rather horizontal, indicating its weak dependence on E_kin.

Fig. 9. Ejecta mass Mej (top) and its fraction of the total mass Mej/(M1+M2) (bottom), at 0.5 d after the collision as a function of initial kinetic energy Ekin,coll. We use different markers to distinguish collisions of different masses, while different colors represent different impact parameters. The horizontal dashed lines in the bottom panel correspond to estimates using the fitting formula in Eq. (1). The empty crosses indicate Mej for giant collisions considered in Ryu et al. (2024b).

These three mass ranges can also be expressed in terms of the fractional mass M_ej/(M₁+M₂), as shown in the bottom panel of Fig. 9. Generally, M_ej/(M₁+M₂)≳0.1 for b = 0.1, ≃0.1 for b = 0.25, and ≲0.05 for b = 0.5. The nearly constant mass fraction for each b suggests that M_ej is not directly correlated with M₁+M₂, although M₁+M₂ sets an upper limit on M_ej, i.e., M_ej<M₁+M₂. This plot highlights the effect of the central concentration of the colliding stars. Notably, for the collisions between 1 M_⊙ and 0.9 M_⊙ stars with b = 0.5 (the three blue circles at 0.03≲M_ej/(M₁+M₂)≲0.05), the mass fractions are particularly high compared to other b = 0.5 cases. We attribute this to how much a He-rich core develops within the colliding stars. As a He-rich core forms, the star becomes increasingly resistant to complete destruction. In 1 M_⊙ and 0.9 M_⊙ stars, the He-rich core constitutes a noticeably small fraction compared to more massive stars (see Fig. 2).

As a more extreme case than the relatively massive main-sequence stars with large He-rich cores, collisions between giants (Ryu et al. 2024b; Dessart et al. 2024) do not exhibit a correlation of M_ej with b or E_kin,coll, where M_ej/(M₁+M₂)≳80% for b = 0.02−0.4 (empty crosses in Fig. 9), roughly corresponding to the envelope mass of the giant stars. The different correlations are directly related to the structural differences between main-sequence and giant stars. Because almost entire main-sequence star is subject to being shocked and launched as ejecta upon collision, the ejecta mass generally increases as b decreases. While the central part of the star is more resistant to destruction than the envelope is, the density contrast is usually not large enough to make a significant difference. However, in giant stars, the He-core is highly resistance to destruction due to their compactness (≃10⁻² R_⊙ for 1 M_⊙ giants; see Fig. 2 in Ryu et al. 2024b) and high binding energy (≳10⁴⁹ erg for 1 M_⊙ giants). For example, collisions that can result in even partial destruction of the cores in two 1 M_⊙ giants occur only when b≲O(10⁻²) and V_rel≳5000 km s⁻¹. Consequently, most of the ejecta mass produced in giant collisions originally belonged to the envelope. Furthermore, because the envelope is relatively loosely bound, it can be easily stripped away by the collision, resulting in the weak dependence on b. The large size of the envelope in giants also leads to a more efficient conversion of kinetic energy to radiation energy, producing larger luminosities, which will be shown in the following section.

The dependence of M_ej/(M₁+M₂) on b for collisions between main-sequence stars can be fitted with the formula

$$\frac{M_{\mathrm{ej}}}{M_1+M_2}=exp\left[-0.63{\left(\frac{b}{0.1}\right)}^{1.3}\right],$$(1)

with a 1σ uncertainty of ≃40−50%⁶.

4.2. Ejecta energy

The ejecta kinetic energy E_kin,ej exhibits a similar trend with b as M_ej, yielding a larger E_kin,ej for more head-on collisions. However, its correlation with E_kin is stronger, following E_kin,ej∝E_kin, as depicted in the top panel of Fig. 10. As a consequence of the linearity between E_kin,ej and E_kin, their ratio is nearly independent of E_kin,ej, which is shown in the middle panel. In most of our models, E_kin,ej/E_kin≃0.6 for b = 0.1, ≃0.1 for b = 0.25, and ≃0.01 for b = 0.5. Accounting for the dependence on b, the correlation between E_kin,ej and E_kin for 0.1≲b≲0.5 can be expressed as,

Fig. 10. Same as Fig. 9, but for the kinetic energy Ekin,ej (top and middle) and radiation energy Erad,ej of the ejecta. The dashed lines indicate estimates by Eq. (2) for Ekin,ej and Eq. (3) for Erad,ej. Solid markers represent energy values directly obtained from 3D simulations, while open markers (positioned above or below the solid markers) indicate those from 1D profiles used for radiative-transfer simulations. The empty cross markers depict values for giant collisions considered in Ryu et al. (2024b).

$$E_{\mathrm{kin},\mathrm{ej}}=E_{\mathrm{kin}}exp\left[-0.5{\left(\frac{b}{0.1}\right)}^{1.4}\right],$$(2)

with a 1σ uncertainty of ≃25% (for b≲0.25) and ≃60% (for b = 0.5). Eq. (2) may not satisfactorily describe the correlation outside the considered range of b, in particular, for b<0.1. For b>0.5, we expect Eq. (2) to describe the continuous decay of E_kin,ej beyond b = 0.5 down to zero. However, it also predicts that E_kin,ej⟶E_kin as b⟶0. While this remains reasonable, our models cannot confirm it.

The radiation energy E_rad,ej trapped in the unbound ejecta shares similar correlations with E_kin as E_kin,ej – namely, larger E_rad,ej for larger E_kin and smaller b. However, E_rad,ej has a steeper correlation with E_kin, i.e., $E_{\mathrm{rad},\mathrm{ej}}\propto E_{\mathrm{kin}}^{1.3}$. The overall dependence of E_rad,ej on E_kin is expressed as,

$$E_{\mathrm{rad},\mathrm{ej}}=6\times {10}^{46}\mathrm{erg}{\left(\frac{E_{\mathrm{kin}}}{{10}^{50}\mathrm{erg}}\right)}^{1.3}exp\left[-0.6{\left(\frac{b}{0.1}\right)}^{1.1}\right].$$(3)

with a 1σ uncertainty of 30−60%.

The fraction of E_kin,coll converted into E_kin,ej for giant collisions (empty cross markers in the top and middle panels of Fig. 10) is consistently high, independent of b, i.e., E_kin,ej/E_kin,coll0.5−0.7, which corresponds to the highest fraction for main-sequence collisions with b = 0.1. Its dependence of E_kin,ej on E_kin,coll is reasonably described by Eq. (3) for b = 0.1. On the other hand, E_rad,ej for giant collisions is up to two orders of magnitude larger than that for main-sequence collisions at a given E_kin,coll. This indicates a greater conversion of E_kin,coll to E_rad,ej and brighter events in giant collisions (see Dessart et al. 2024).

4.3. Ejecta observables

Figure 11 shows three strong correlations with observables at 1 d. As illustrated in the top panels and bottom-left panel, the bolometric luminosity at 1d L_bol,1d ranges from 5×10³⁹ erg s⁻¹ to 3×10⁴² erg s⁻¹. More energetic collisions with smaller impact parameters (so more head-on) generate brighter events (L_bol,1d∝E_kinb⁻¹). This is expected because the collision kinetic energy sets the energy budget available for conversion into the radiation energy of the ejecta, while the impact parameter determines the amount of mass that can be shocked at the collision, so the energy conversion efficiency. However, the brightness of the events appears to show no correlations with either the stellar masses or the relative velocity when considered individually. The correlations with the original kinetic energy suggest that it may be possible to infer the kinetic energy from the observables. However, the lack of correlations with the masses and the relative velocity indicates that it would be challenging to infer them separately based solely on the luminosity or magnitudes.

Fig. 11. Bolometric luminosity Lbol,1d (top-left) and UVW1-band magnitudes MUVW1,1d (top-right), V-band magnitudes MV,1d (bottom-left) as a function of initial kinetic energy Ekin,coll and the Doppler velocity VHα,1d (bottom-right) at maximum absorption measured in the Hα line profile as a function of Vrel. All these values are measured at 1 d since collision. We distinguish collisions with different masses using different markers and those with different impact parameters using different colors. The black and blue diagonal lines show the relations given by the fitting formulae in Equations (4)–(7) for b = 0.1 and 0.5, respectively. In the top panels and bottom-left panel, the empty crosses represent values for giant collisions considered in Ryu et al. (2024b).

Note that the correlation between the event's brightness (L_bol,1d, M_UVW1,1d and M_V,1d) and E_kin for b = 0.25 appears noticeably weaker than that for the other two cases, b = 0.1 and 0.5. In other words, the scatter for b = 0.25 in the top panels and bottom-left panel is distributed more horizontally. We suspect that this weak correlation is largely due to the small sample size for b = 0.25 (only six cases, compared to 17 cases for b = 0.1 and 11 cases for b = 0.5).

In comparison, giant collisions, with their high values of E_kin,ej/E_kin,coll, produce significantly brighter transients. As shown in the top panels and bottom-left panels of Fig. 11 using empty cross markers, L_bol,1d for giant collisions is higher by up to two orders of magnitude and both M_UVW1,1d and M_V,1d are brighter by a factor of 2−3 mag.

However, the positive correlation between the Doppler velocity at maximum absorption measured in the Hα line profile and the relative collisional velocity, shown in the bottom panel of Fig. 11, may help break the degeneracy between the masses and relative velocity. This correlation, V_Hα,1d∝V_rel, can be understood by the fact that the debris expansion speed scales with the initial relative velocity, as has been also found in collisions involving giants (Ryu et al. 2024b). Consequently, accurate measurements of the Hα velocity, along with either the luminosity or magnitudes, may place valuable constraints on collision parameters, particularly the initial relative velocity and the masses (e.g., the reduced mass if the center-of-mass of the two colliding stars was initially negligible compared to the relative velocity). This implies that spectra are essential for parameter inference.

The fitting formulae that describe the correlations are

$$L_{\mathrm{bol},1\mathrm{d}}=8\times {10}^{40}\mathrm{erg}{\left(\frac{b}{0.1}\right)}^{-1}\left(\frac{E_{\mathrm{kin}}}{{10}^{50}\mathrm{erg}}\right),$$(4)

$$M_{\mathrm{UVW}1,1\mathrm{d}}=-2{log}_{10}\left(\frac{E_{\mathrm{kin}}}{{10}^{50}\mathrm{erg}}\right)-14.5{\left(\frac{b}{0.1}\right)}^{-0.1},$$(5)

$$M_{\mathrm{V},1\mathrm{d}}=-1.6{log}_{10}\left(\frac{E_{\mathrm{kin}}}{{10}^{50}\mathrm{erg}}\right)-12.5{\left(\frac{b}{0.1}\right)}^{-0.07},$$(6)

$$V_{\mathrm{H}\alpha ,1\mathrm{d}}=5500{\left(\frac{b}{0.1}\right)}^{-0.25}\left[1.3\left(\frac{V_{\mathrm{rel}}}{5000{\mathrm{km}\mathrm{s}}^{-1}}\right)+1\right].$$(7)

The 1σ uncertainty is 36%, 2−3%, and 10% for L_bol,1d, both of the magnitudes, and V_Hα,1d, respectively.

5. Summary and conclusion

In this paper, we investigated the properties of ejecta produced in high-velocity collisions between terminal-age main-sequence stars in galactic nuclei and their observables. We adopt the same methodology as in Dessart et al. (2024), combining 3D moving-mesh hydrodynamics and 1D radiation-transfer simulations. We consider a wide range of parameters for collisions, namely, stellar masses of 0.9−10 M_⊙, mass ratios of 0.1−0.9, impact parameters of (0.1−0.5)(R₁+R₂), and initial relative velocities of 2500−10 000 km s⁻¹.

When two main-sequence stars collide at high velocities, strong shocks form, driving quasi-spherically expanding ejecta. In most of our models, the collisions leave behind single He-rich remnants. If we use stars are at an earlier evolutionary stage (e.g., bona fide main-sequence stars), they are more likely to be completely destroyed because their core is less gravitationally bound.

Because of the large stored radiation energy, the debris from high-velocity stellar collisions release copious amounts of electromagnetic radiation. Their bolometric luminosities range from 5×10³⁹ erg s⁻¹ to 3×10⁴² erg s⁻¹ at 0.5 d after the collision, during which the flux falls in the UV. If the rate of fading at 0.5–1 d (L_bol∝t⁻²−t⁻⁴) holds at earlier times, the peak luminosity at the moment of collision might be significantly higher (potentially by a factor of 10 or 2 magnitude brighter), possibly reaching a brightness level comparable to tidal disruption events (Gezari 2021; Wevers & Ryu 2023). However, the peak duration is typically short on the order of O(0.1) d, so transients surveys with cadence less than a day would be able to observe these events at their luminosity peak. Even after a photometric discovery, prompt spectroscopic follow-up is essential to confirm the nature of the events. Within a week, the luminosity drops by more than an order of magnitude and levels off on a plateau, at which point the spectrum has shifted to the optical.

The observables of these collisions are primarily determined by the mass and energy content of the debris at the onset of their expansion. The ejecta mass depends mostly on the impact parameter b, with the fractional mass M_ej/(M₁+M₂) increasing as b decreases (∝exp[−b^1.3] in Eq. 1). On the other hand, the energy content of ejecta and the transient brightness (e.g., luminosity and magnitude) generally correlate positively with b and the initial kinetic energy E_kin. For example, at one day, the radiation energy of the expanding ejecta scales as $E_{\mathrm{rad},\mathrm{ej}}\propto E_{\mathrm{kin}}^{1.3}exp[-b^{1.1}]$ (Eq. 3) and the bolometric luminosity L_d,1d∝E_kinb⁻¹ (Eq. 4). The correlation with E_kin, rather than the masses of the colliding stars and initial relative velocity individually, suggests a degeneracy between the stellar masses and the relative velocity when inferred from observations. This degeneracy can be broken with velocity measurements done on H I Balmer lines like Hα, as those lines depend most strongly on the relative velocity and not on the stellar masses.

In conclusion, our detailed hydrodynamics and radiative transfer calculations suggest that transients generated by high-velocity collisions in galactic nuclei are observable nuclear transients. This is expected for collisions between main-sequence stars (this work) as well as between giant stars (Ryu et al. 2024b; Dessart et al. 2024). However, some of our assumptions require further investigation. For example, the quasi-spherical expansion found in our hydrodynamics simulations would be disrupted near the supermassive black hole, where tidal forces are strong enough to deform the expanding debris. In such cases, observable properties would depend on the viewing angle. Additionally, the supermassive black hole would influence the debris on a longer times scale, extending the optically thick and luminous phase, which must be thoroughly examined to develop a complete picture of transients associated with high-velocity collisions.

Data availability

All radiative-transfer simulations presented in this work are available at https://zenodo.org/communities/snrt

Acknowledgments

We thank the anonymous referee for kindly reviewing the manuscript and providing constructive comments. This research project was conducted using computational resources (and/or scientific computing services) at the Max-Planck Computing & Data Facility. The authors gratefully acknowledge the scientific support and HPC resources provided by the Erlangen National High Performance Computing Center (NHR@FAU) of the Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU) under the NHR projects b166ea10 and b222dd10. NHR funding is provided by federal and Bavarian state authorities. NHR@FAU hardware is partially funded by the German Research Foundation (DFG) – 440719683. In addition, some of the simulations were performed on the national supercomputer Hawk at the High Performance Computing Center Stuttgart (HLRS) under the grant number 44232. This work was supported in France by the “Programme National Hautes Energies” of CNRS/INSU co-funded by CEA and CNES. This work was granted access to the HPC resources of TGCC under the allocation 2023 – A0150410554 on Irene-Rome made by GENCI, France.

References

Appendix A: Model parameters

Tab. A.1. summarizes the model parameters for the collision configurations and corresponding ejecta characteristics.

Appendix B: Light curves and spectra

We provide the bolometric luminosity (Fig. B.1.), UVW1−band magnitude (Fig. B.2.), V−band magnitude (Fig. B.3.) for all our models as a function time. In addition, we also present their spectral evolution in Figs. B.4.–B.9..

Fig. B.1. Bolometric light curves for our six sets of models.

Fig. B.2. UVW1-band light curves for our six sets of models.

Fig. B.3. V-band light curves for our six sets of models.

Fig. B.4. Spectral montage for the modelset # 1.

Fig. B.5. Spectral montage for the modelset # 2.

Fig. B.6. Spectral montage for the modelset # 3.

Fig. B.7. Spectral montage for the modelset # 4.

Fig. B.8. Spectral montage for the modelset # 5.

Fig. B.9. Spectral montage for the modelset # 6.

All Tables

All Figures

Fig. 1. Ranges of the relative velocity between stars and stellar density for destructive (magenta-shaded) and non-destructive (gray-shaded) collisions, separated by the surface escape velocity of a 1 M⊙ main-sequence star. The velocity v and distance R from the black hole are related by $v/c=\sqrt{r_{\mathrm{g}}/R}$, with the black-shaded region (“Inside black hole”) indicating v/c>1. Here, rg is the gravitational radius. Black diagonal lines indicate nuclear clusters with the same mass as the black hole, assuming a Bahcall & Wolf (1976) profile extending inward (arrow) from the influence radius until 2 M⊙ is enclosed or the event horizon is reached. The red and blue crosses indicate the observationally constrained density at 0.25 pc from the Galactic center (Genzel et al. 2010) and the semimajor axis of the S2 star (Gillessen et al. 2009), respectively.

In the text

	Fig. 2. Internal structure of stellar models. Density (top) and H mass fraction (bottom) profiles of our main-sequence stellar models with a metallicity of Z = 0.02 and at the core H mass fraction of 0.01, as a function of enclosed mass (left) and radius (right).
In the text

	Fig. 3. Gas density distribution at successive moments in main-sequence star collisions with M1 = 1 M⊙, M2 = 0.9 M⊙, Vrel = 5000 km s−1, but different impact parameters. These AREPO simulations of stellar collisions are characterized by b = 0.1, 0.25, and 0.5 (from top to bottom) and shown here at t = 0.0, 0.01, 0.05, and 0.5 d (from left to right).
In the text

	Fig. 4. Same as Fig. 3, but for collisions with M1 = 10 M⊙ and M2 = 1 M⊙, Vrel = 5000 km s−1. In all cases, the collision unbinds only a fraction of the envelopes of both stars, leaving behind remnants.
In the text

Fig. 5. Evolution of the ejecta kinetic and radiative energies in a sample of stellar-collision simulations. We show radiation energy Erad,ej (top) and kinetic energy Ekin,ej (bottom) for modelsets #1 (left, 1 M⊙+0.9 M⊙) and #5 (right, 10 M⊙+1 M⊙). The legend indicates the initial relative velocity (e.g., vrel5 means Vrel = 5000 km s−1) and impact parameter (e.g., b0p5 means b = 0.5). The diagonal black lines in the top panels indicate the power law of t−1, expected for an adiabatic expansion of a radiation-dominated plasma.

In the text

	Fig. 6. Photometric properties for main-sequence star collisions from the modelset # 1. From top to bottom, we show the bolometric, UVW1-band, and V-band light curves. Filled dots indicate the corresponding quantities inferred for SN 1987A (after correction for reddening and distance). Corresponding results for all sets are shown in Appendix B.
In the text

	Fig. 7. Spectral evolution from the UV to the optical and from ∼0.5 until ∼7.2 d for model ms1/ms0p9/vrel10/b0p1 (modelset # 1). Fλ represents the spectral flux density or the flux per unit wavelength interval. Corresponding results for all models from the six sets are shown in Appendix B.
In the text

	Fig. 8. Spectral comparisons of models grouped by the same values of Vrel and b but differing in the stellar masses of the colliding pair. The time is one day after the onset of the collision. To better reveal the weaker flux in the optical, we show the quantity λ2Fλ, which is analogous to an Fν.
In the text

Fig. 9. Ejecta mass Mej (top) and its fraction of the total mass Mej/(M1+M2) (bottom), at 0.5 d after the collision as a function of initial kinetic energy Ekin,coll. We use different markers to distinguish collisions of different masses, while different colors represent different impact parameters. The horizontal dashed lines in the bottom panel correspond to estimates using the fitting formula in Eq. (1). The empty crosses indicate Mej for giant collisions considered in Ryu et al. (2024b).

In the text

Fig. 10. Same as Fig. 9, but for the kinetic energy Ekin,ej (top and middle) and radiation energy Erad,ej of the ejecta. The dashed lines indicate estimates by Eq. (2) for Ekin,ej and Eq. (3) for Erad,ej. Solid markers represent energy values directly obtained from 3D simulations, while open markers (positioned above or below the solid markers) indicate those from 1D profiles used for radiative-transfer simulations. The empty cross markers depict values for giant collisions considered in Ryu et al. (2024b).

In the text

Fig. 11. Bolometric luminosity Lbol,1d (top-left) and UVW1-band magnitudes MUVW1,1d (top-right), V-band magnitudes MV,1d (bottom-left) as a function of initial kinetic energy Ekin,coll and the Doppler velocity VHα,1d (bottom-right) at maximum absorption measured in the Hα line profile as a function of Vrel. All these values are measured at 1 d since collision. We distinguish collisions with different masses using different markers and those with different impact parameters using different colors. The black and blue diagonal lines show the relations given by the fitting formulae in Equations (4)–(7) for b = 0.1 and 0.5, respectively. In the top panels and bottom-left panel, the empty crosses represent values for giant collisions considered in Ryu et al. (2024b).

In the text

	Fig. B.1. Bolometric light curves for our six sets of models.
In the text

	Fig. B.2. UVW1-band light curves for our six sets of models.
In the text

	Fig. B.3. V-band light curves for our six sets of models.
In the text

	Fig. B.4. Spectral montage for the modelset # 1.
In the text

	Fig. B.5. Spectral montage for the modelset # 2.
In the text

	Fig. B.6. Spectral montage for the modelset # 3.
In the text

	Fig. B.7. Spectral montage for the modelset # 4.
In the text

	Fig. B.8. Spectral montage for the modelset # 5.
In the text

	Fig. B.9. Spectral montage for the modelset # 6.
In the text