Fates of Rotating Supergiants from Stellar Mergers and the Landscape of Transients upon Core-collapse

TL;DR

This study models rotating supergiants from binary mergers to understand their evolution and explosion outcomes, linking them to specific supernova types, gamma-ray bursts, and luminous transients.

Contribution

It provides a detailed grid of post-merger stellar models with rotation, exploring their core-collapse behavior and associated transient phenomena, a novel approach to merger product evolution.

Findings

Merger mass gain influences supergiant type (blue vs. red).

Rotating BSGs can produce SN 1987A-like supernovae.

Failed explosions lead to fallback accretion and potential jets.

Abstract

We present a grid of rotating supergiant models from post-main sequence binary merger products, constructed by the MESA stellar evolution code. We focus on the evolution of these stars until core-collapse, in addition to their rotation, which could influence their mass loss and explosion phenomenology. We find that (i) as in previous studies, larger mass gain by merger favors the production of blue supergiants (BSGs) over red supergiants, and (ii) merger products that end as BSGs at core collapse have rotating outer envelopes, with lower-mass BSGs having faster envelope rotation due to less wind mass loss. We model the expected transients from these BSGs upon core-collapse, considering cases where the neutrino-driven explosion is successful and unsuccessful. The successful explosions result in supernovae (SNe) with long-rising light curves of morphology similar to SN 1987A. Failed explosions of these BSGs result in envelope fallback of $\sim (0.1$- several) $~M_\odot$ over $10^3$-$10^5$ seconds that power strong ($10^{51}$-$10^{53}$ erg) accretion-driven outflows in winds and possibly jets, with relativistic jets (if formed) generally capable of breaking out of the BSG envelope. Our modeling points to these merger-origin BSGs as viable progenitors for SN 1987A-like SNe, ultra-long gamma-ray bursts, and some of the fast luminous transients found in high-cadence optical surveys.

Tables

Table 1. Table 1: Summary of our 32 models of single and merger-remnant stars computed by MESA and their post-main sequence properties.

$M_{1,\mathrm{ZAMS}}$ [$M_{\odot }$]	$M_2$ [$M_{\odot }$]	$M_{\mathrm{fin}}$ [$M_{\odot }$]	$M_{\mathrm{He},\mathrm{fin}}$ [$M_{\odot }$]	$R_{\mathrm{fin}}$ [$R_{\odot }$]	$T_{\mathrm{eff},Y_{\mathrm{He}}=0.5}$ [K]	$T_{\mathrm{eff},\mathrm{fin}}$ [K]	$⟨j_{\mathrm{rot}}⟩$ [${10}^{17}{\mathrm{cm}}^2{s}^{-1}$]	$M_{j>j_{\mathrm{ISCO},\mathrm{Sch}}}$ [$M_{\odot }$]
10	0	9.6	3.0	658	3720	3337	1.90	5.28
10	2	11.1	3.0	644	3768	3387	2.97	6.92
10	4	12.9	3.3	678	3808	3410	3.64	8.33
10	6	14.9	3.5	698	19116	3457	3.55	9.56
10	8	16.6	3.6	34	21207	15466	3.51	8.42
15	0	13.6	4.8	992	3572	3318	1.92	6.38
15	2	15.0	4.9	997	3578	3339	1.64	6.27
15	4	17.8	5.4	1050	14570	3363	1.64	6.70
15	6	19.6	5.2	971	19776	3436	2.11	8.52
15	8	21.4	5.1	50	21805	15128	2.13	5.47
15	10	23.2	4.9	36	23568	18761	2.03	5.41
15	12	25.2	4.7	29	25156	21367	2.09	5.69
20	0	17.7	7.2	1346	8775	3318	0.05	0.00
20	2	19.9	7.3	1338	12047	3341	0.11	0.00
20	4	22.2	7.4	1324	15154	3370	0.36	0.00
20	6	24.4	7.2	1254	18800	3424	0.99	1.05
20	8	26.1	6.8	139	20811	10356	1.01	1.80
20	10	27.9	6.7	81	22445	13960	1.06	1.93
20	12	29.7	6.4	59	23966	16809	1.11	1.89
20	14	31.7	6.5	51	25283	18434	1.23	2.07
20	16	33.7	6.2	45	26458	20124	1.28	1.96
25	0	20.6	9.5	1602	8354	3338	0.06	0.00
25	2	23.4	9.5	1593	11224	3355	0.07	0.00
25	4	26.4	9.5	1563	14188	3385	0.08	0.00
25	6	28.8	9.3	1502	16911	3427	0.19	0.00
25	8	30.6	9.0	1460	19031	3454	0.45	0.00
25	10	32.4	8.8	445	20681	6436	0.41	0.35
25	12	34.1	8.6	219	21995	9400	0.48	0.23
25	14	35.9	8.5	158	23326	11256	0.53	0.22
25	16	38.0	8.2	131	24380	12559	0.65	0.33
25	18	40.0	8.3	108	25389	14103	0.73	0.40
25	20	42.0	8.3	89	26275	15898	0.80	0.42

Equations

$$\dot{M}=\max \left[{\dot{M}}_{\mathrm{init}},\frac{{\dot{M}}_{\mathrm{add}}}{1+(\Delta M_{\mathrm{crit}}/\Delta M)}\right]$$

$$j_{\mathrm{acc}}$$

$$j_{\mathrm{crit}}$$

$$j_{\mathrm{surf}}\sim (0.3-2)\times 10^{19}\mathrm{c}{\mathrm{m}}^2{\mathrm{s}}^{-1}{\left(\frac{⟨\sin i⟩}{\pi /4}\right)}^{-1}\left(\frac{R_{\ast }}{50R_{\odot }}\right),$$

$$⟨j_{\mathrm{rot}}⟩=\frac{1}{M_{\ast }}{\int }_0^{M_{\ast }}j(M_r)d{M}_r,$$

$$\frac{d(⟨j_{\mathrm{rot}}⟩M_{\ast })}{dt}=-\frac{d{M}_{\ast }}{dt}\left(\frac{2}{3}\Omega R_{\ast }^2\right).$$

$$\frac{d⟨j_{\mathrm{rot}}⟩/dt}{⟨j_{\mathrm{rot}}⟩}=\left(\frac{1}{k^2}-\frac{2}{3}\right)\frac{d{M}_{\ast }/dt}{M_{\ast }},$$

$$⟨j_{\mathrm{rot}}⟩=j_0{\left(\frac{M_{\ast }}{M_{\ast ,0}}\right)}^{(1/k^2-2/3)}.$$

$$j_{\mathrm{ISCO}}(M)$$

$$r_{\mathrm{ISCO}}^{\prime }$$

$$z_1$$

$$z_2$$

$$t_{\mathrm{ff}}\approx \sqrt{\frac{{\pi }^2{R}_{\ast }^3}{8G{M}_{\ast }}}\sim 1\mathrm{day}{\left(\frac{R_{\ast }}{50R_{\odot }}\right)}^{3/2}{\left(\frac{M_{\ast }}{30M_{\odot }}\right)}^{-1/2}.$$

$$r_{\mathrm{disk}}\approx \frac{j^2}{G{M}_{\mathrm{BH}}}\sim 4\times 10^8\mathrm{cm}{\left(\frac{j}{10^{18}\mathrm{c}{\mathrm{m}}^2{\mathrm{s}}^{-1}}\right)}^2{\left(\frac{M_{\mathrm{BH}}}{20M_{\odot }}\right)}^{-1},$$

$$t_{\mathrm{visc}}$$

$${\dot{M}}_{\mathrm{wind}}(t)$$

$$L_{\mathrm{wind}}(t)$$

$${\dot{M}}_{\mathrm{BH}}(t)={\dot{M}}_{\mathrm{acc}}-{\dot{M}}_{\mathrm{wind}}=\frac{M_d}{t_{\mathrm{visc}}}{\left(\frac{r_{\mathrm{ISCO}}}{r_{\mathrm{disk}}}\right)}^s.$$

$$E_{\mathrm{rot}}=\left(1-\sqrt{\frac{1+\sqrt{1-a_{\mathrm{BH}}^2}}{2}}\right)M_{\mathrm{BH}}{c}^2,$$

$$t_{\mathrm{oc}}\approx \sqrt{\frac{\kappa M_c}{4\pi v_{c}c}}\sim 4\mathrm{day}{\left(\frac{v_{\mathrm{c}}}{0.1c}\right)}^{-1/2}{\left(\frac{M_{\mathrm{c}}}{0.2M_{\odot }}\right)}^{1/2},$$

$$L_{\mathrm{oc}}$$

$$T_{\mathrm{eff},\mathrm{oc}}$$

$$t_{\mathrm{diff}}$$

$$r_{\mathrm{inj}}\approx \max (R_{\ast },v_{\mathrm{ej}}{t}_{\mathrm{acc}}).$$

$$L$$

$$L_{\mathrm{jet}}(t)=\frac{1}{2}{\eta }_{\mathrm{jet}}{\dot{M}}_{\mathrm{BH}}(t)c^2$$

$$\frac{d{r}_{\mathrm{h}}}{dt}=\frac{N_{s}c}{1+{\tilde{L}}^{-1/2}}$$

$$\tilde{L}(t)=\frac{L_{\mathrm{jet}}(t^{\prime }=t-r_{\mathrm{h}}/c)}{{\Sigma }_{\mathrm{jet}}{\rho }_{\ast }(r=r_{\mathrm{h}})c^3},$$

$$E_c(t)\approx {\int }_0^{t-r_h/c}d{t}^{\prime }{L}_{\mathrm{jet}}(t^{\prime })$$

$$\frac{d{r}_{\mathrm{c}}}{dt}=\sqrt{\frac{P_c}{{\bar{\rho }}_{\ast }}},$$

$${\bar{\rho }}_{\ast }=\frac{M_c}{V_c}\approx \frac{{\int }_{r_{h,0}/c}^{t}d{t}^{\prime }[\pi r_c^2{\rho }_{\ast }{v}_h]}{{\int }_{r_{h,0}/c}^{t}d{t}^{\prime }[\pi r_c^2{v}_h]},$$

$$P_c(t)=\frac{E_c}{3V_c}.$$

Full text

††thanks: ∗E-mail: tsuna@caltech.edu

1 Introduction

Massive stars are observed in the local universe in large numbers, occupying diverse regions of the Hertzsprung-Russell (HR) diagram. The diversity is shaped by various processes in the star, such as nuclear burning, convection, radiation transport and mass loss. These stars also often exist in binaries or higher-order multiples (Duchêne & Kraus, 2013; Sana et al., 2014), and binary interactions can drastically affect the evolutionary paths throughtout their lives (Sana et al., 2012).

In the local universe, there has been a long-standing puzzle on the abundant population of blue supergiants (BSGs; e.g., Humphreys & McElroy 1984; Fitzpatrick & Garmany 1990; Castro et al. 2014; Patrick et al. 2025). They reside in a region of the HR diagram where single star evolution theory predicts a gap (Hertzsprung gap), due to the rapid evolution of these stars to red supergiants (RSGs) once they leave the main sequence. Among the channels proposed for long-lived BSGs, a promising one is mergers/accretion in binaries (Braun & Langer, 1995; Glebbeek et al., 2013; Vanbeveren et al., 2013; Justham et al., 2014; Farrell et al., 2019; Menon et al., 2024; Schneider et al., 2024), which generally leads to stars with undermassive helium cores and potentially higher surface helium/nitrogen abundance due to dredge-up of the burned layer to the surface.

While there have been successful attempts in reproducing the BSG population via binary merger products, the final fates of these stars are less explored. The best-studied case is SN 1987A, an explosion of a BSG in the Large Magellanic Cloud (LMC; West et al., 1987). The equatorial ring seen in the supernova (SN) remnant phase supports the BSG progenitor being a product of a binary merger that occurred $\sim 10^{4}$ yrs before core-collapse (e.g., Hillebrandt & Meyer, 1989; Podsiadlowski et al., 1990; Podsiadlowski, 1992; Morris & Podsiadlowski, 2009; Menon & Heger, 2017; Urushibata et al., 2018). On the other hand, failed SNe of rotating BSGs result in fallback accretion onto the compact remnant, and have also been suggested as progenitors of various high-energy transients (Woosley & Heger, 2012; Dexter & Kasen, 2013; Kashiyama et al., 2013; Nakauchi et al., 2013; Kashiyama & Quataert, 2015; Perna et al., 2018; Margutti et al., 2019; Tsuna et al., 2021). While stellar mergers provide a natural pathway for forming rotating supergiants, detailed modelling of such transients from these merger products has not been systematically explored.

In this work we model rotating evolved supergiants from post-main sequence stellar mergers, aiming to obtain a comprehensive understanding of their final fates. We couple a 1D framework of constructing and evolving merger products with prescriptions for angular momentum (AM) transport inside the star. Using these stellar models, we consider diverse transients that can be expected from the deaths of these stars. We find that BSGs are produced by gain of a significant fraction of its mass due to accretion. Furthermore, the merger products tend to maintain rapid envelope rotation until core-collapse, which can produce not only 1987A-like SNe but also accretion-driven high-energy transients in case of a failed SN.

In Section 2 we first describe our modeling of the stellar merger products and their evolution towards core-collapse, and in Section 3 discuss the diverse evolution of the structure and rotation of these stars until core-collapse. We then discuss in Section 4 the expected transients upon core-collapse using our stellar models, showing that these post-main sequence merger products are viable progenitors for various classes of transients, from 1987A-like SNe to engine-powered high-energy transients such as ultra-long gamma-ray bursts (GRBs) and fast luminous transients like AT2018cow. We conclude in Section 5, with directions for future studies.

2 Stellar Evolution Models

We evolve massive stars with the open-source stellar evolution code MESA r23.05.1 (Paxton et al., 2011, 2013, 2015, 2018, 2019; Jermyn et al., 2023)111A Zenodo link for the MESA inlists and the stellar models at core-collapse will be released upon acceptance of this manuscript.. Throughout this work we fix the stars to a metallicity of $Z=0.007$, comparable to the environment of the LMC. This is also close to the mean/median metallicities inferred from observations of the explosion sites of local core-collapse SNe (12+$\log_{10}({\rm O/H})\approx$ 8.35–8.5; Galbany et al., 2016; Pessi et al., 2023), and could be regarded as representative of the local core-collapse SN population.

We consider rigidly rotating stars at zero age main sequence (ZAMS), with surface rotation velocity of 20% of the critical velocity. The resulting surface rotation speeds of $\sim 100$ km s*-1* are similar to those observed for OB stars in LMC (Hunter et al., 2008; Ramachandran et al., 2018, 2019). We evolve the rotation profile of the star based on recent prescriptions of AM transport (Fuller et al., 2019; Ma & Fuller, 2019), with parameters adopted in Fuller & Lu (2022). Mass loss is included under the Dutch wind prescription (de Jager et al., 1988; Nugis & Lamers, 2000; Vink et al., 2001), with an efficiency (Dutch_scaling_factor) of 0.5. The mass loss rate also depends on the surface rotation velocity, as described in Section 6.4 of Paxton et al. (2013).

For convection we adopt a mixing length parameter of $\alpha_{\rm MLT}=1.5$, and a Ledoux criterion with semi-convection efficiency parameter $\alpha_{\rm sc}=10$ motivated from the properties of post-main sequence massive stars in the Small Magellanic Cloud (Schootemeijer et al., 2019)222This high $\alpha_{\rm sc}$, close to adopting a Schwarzschild criterion, is inferred from the observed fraction of BSGs in the SMC (Schootemeijer et al., 2019). A high $\alpha_{\rm sc}$ as in the Schwarzschild criterion is also argued from independent stellar modeling, as a consequence of entrainment of the semi-convective zone by the convective zone (e.g., Anders et al., 2022) and/or double-diffusive instabilities occurring in semi-convective regions (e.g., Moore & Garaud, 2016)].. Convective core overshooting is taken into account as an exponential scheme in MESA with parameters $(f,f_{0})=(0.02,0.005)$. We turn off overshooting after core helium depletion, by setting the parameter min_overshoot_q in MESA to be a value above unity.

We mimic the merger process in MESA by rapid mass accretion to the primary, as done in previous works (e.g., Justham et al., 2014; Rui & Fuller, 2021; Schneider et al., 2024; Menon et al., 2024). We consider the most likely case of (early) Case B merger products, where merger occurs during when the star rapidly expands across the Hertzsprung gap (see e.g. Figure 1 of Justham et al., 2014). When the star expands to a radius $R_{*}=50\ R_{\odot}$, we add mass with a rate of $\dot{M}_{\rm add}\approx 10^{-2}\ M_{\odot}\ {\rm yr}^{-1}$ until it accretes the mass of the companion (using max_star_mass_for_gain=$M_{2}$). This choice of $R_{*}$ is motivated by the primary to retain most of its envelope upon merger, before it expands too much and becomes loosely bound such that it may be ejected by common envelope evolution (e.g., Klencki et al., 2021). In Appendix B we verify by running additional MESA models that the final outcome is not strongly sensitive to the detailed choice of the radius at the onset of mass accretion.

To make the accretion process numerically stable in MESA, we initially add mass at a lower rate $\dot{M}_{\rm init}=10^{-5}\ M_{\odot}\ {\rm yr}^{-1}$, and then evolve $\dot{M}$ to asymptote to the desired $\dot{M}_{\rm add}$ when the star roughly accretes a mass of $\Delta M_{\rm crit}$. We prescribe $\dot{M}$ as a function of cumulative accreted mass $\Delta M$ by

[TABLE]

where we take $\Delta M_{\rm crit}=10^{-2}~M_{\odot}$. In MESA, mass being added is set to be the same specific entropy as that at the surface (Paxton et al., 2015) due to the short thermal timescale expected at the photosphere. We also assume that the composition of the added mass is same as that at the photosphere of the primary, which is almost identical to the composition at ZAMS (see Justham et al. 2014; Schneider et al. 2024 and Section 5 for discussions on this assumption).

For completeness we consider mergers with mass ratios $q=M_{2}/M_{1}$ up to $0.8$, which may occur when the secondary expands and also fills its own Roche Lobe upon mass transfer from the primary (Pols, 1994; Justham et al., 2014). However, we note that a relatively high $q$ may result in stable mass transfer that avoid merger depending on the donor/accretor properties and uncertain stability boundary (e.g., Temmink et al., 2023; Henneco et al., 2024; Ercolino et al., 2024). We assume that the entire secondary’s mass $M_{2}$ is accreted, and no mass is lost from the system upon merger. If one assumes a fraction $f_{\rm acc}<1$ of the secondary mass ends up accreted onto the primary, our models are essentially equivalent to merging a secondary with mass $M_{2}/f_{\rm acc}$ (a similar argument also applies if a fraction of the primary’s mass is lost during the merger process). We prescribe the specific AM of the accreted material $j_{\rm acc}$ by the orbital AM being fully inherited to the merger product,

[TABLE]

where $G$ is the gravitational constant and here $R_{*}=50~R_{\odot}$.

After accretion we run the models until core carbon depletion, defined as when the central carbon mass fraction drops to $10^{-3}$. Afterwards we expect the envelope structure in the MESA models to be mostly unchanged, as the time until core-collapse ($\lesssim 10$ yrs) is much less than the Kelvin-Helmholtz (KH) timescale of supergiants (order $10^{2}$–$10^{3}$ yrs). We adopt the approx21_cr60_plus_co56 nuclear reaction network, appropriate for modeling the star up to this evolutionary stage.

We consider four primaries with initial masses of $M_{1}=[10,15,20,25]~M_{\odot}$, and accreted masses $M_{2}$ in steps of $2~M_{\odot}$ from $0~M_{\odot}$ (i.e. single star model) up to $0.8M_{1}$. This results in a total of 32 models, whose properties are summarized in Table 1.

3 Results

We first show evolutions of representative stellar models starting from a $15~M_{\odot}$ and $25~M_{\odot}$ primary to illustrate the underlying physics of the evolution of the surface properties and the AM profile. We then present the results for our grid of models to discuss the dependence of the merger product’s fate on both $M_{1}$ and $M_{2}$.

3.1 Evolution of Representative Models

3.1.1 Evolution on the HR Diagram

Figure 1 shows the evolutionary tracks of the representative stellar models in the HR diagram. The solid lines in the two panels show the single star models, the dashed and dotted lines show merger models with different $M_{2}$, that ends its life either as a RSG or a BSG.

The single stars expand once they leave the main sequence, with slight differences in the detailed radial evolution for the $15~M_{\odot}$ and $25~M_{\odot}$ models. The former immediately expands and becomes a RSG as it crosses the Hertzsprung gap, while the latter expands more gradually and becomes a RSG near the end of core helium burning. Such variations in radial evolution with mass have been seen in previous post-main sequence stellar models for low metallicities (e.g., Figures 1, C1 of Klencki et al., 2020, also Georgy et al. 2013; Schootemeijer et al. 2019). Despite these subtle differences with mass, all of our single star models eventually become RSGs well before core-collapse, with effective temperatures of $\approx 3300$ K and radii of $\approx 600$–$1600~R_{\odot}$.

The merger models have a more complicated evolution in the HR diagram, whose final fate depends on the mass $M_{2}$ being added to the star. The rapid accretion onto the primary, with accretion timescales much shorter than its KH timescale of the primary, initially inflates the envelope as the star cannot thermally adjust its structure (e.g., Kippenhahn & Meyer-Hofmeister, 1977; Neo et al., 1977; Schürmann & Langer, 2024; Lau et al., 2024). The detailed surface properties at this phase would be sensitive to the prescription for the mass accretion rate.

We expect the envelope structure to be more robust from a few KH times after the end of mass accretion, when the merger product thermally relaxes. The merger models move to the blue side of the HR diagram over a timescale of $10^{3}$–$10^{4}$ years, and initiate core helium burning as a BSG. Models with low $M_{2}$, such as $15M_{\odot}+4M_{\odot}$ and $25M_{\odot}+8M_{\odot}$ in Figure 1, significantly expand soon after core helium depletion, spending the rest of their lives as a RSG. On the other hand, the models with more massive $M_{2}$, such as the $15M_{\odot}+8M_{\odot}$ and $25M_{\odot}+16M_{\odot}$ models, expand less after core helium depletion and remain as a BSG until core-collapse.

Such a behavior of stars turning blue with increasing mass accretion has been explained in the literature in multiple ways. The expansion of massive stars into RSGs have been attributed to the large pressure gradient established around the H-burning shell, when the core contracts soon after main sequence. A large instantaneous mass accretion increases the envelope mass (and gravity) while keeping the He core mass fixed, which makes the star harder to expand while maintaining hydrostatic equilibrium. Another picture was suggested in Schneider et al. (2024), where higher mass accretion increases the temperature at the base of the hydrogen envelope, turning the hydrogen shell-burning layer convective. We also observe the hydrogen shell-burning layer being convective for our BSG models of large $M_{2}$, with reduced pressure gradient around the shell. For both pictures, the resulting structure of such merger products is analogous to a less-evolved massive star in its main-sequence phase, explaining the compact, radiative envelope required for BSGs.

3.1.2 Rotation Profile

Figure 2 shows the specific AM profile for two representative merger product models, of a $15~M_{\odot}$ star accreting $4~M_{\odot}$ and $8~M_{\odot}$. The star before mass accretion has a profile sharply rising at the outermost mass coordinate. The sharp rise is expected for a radiative envelope with density sharply dropping with radius near the surface (i.e. $r$ rapidly increasing with increased $M_{r}$), and obeying rigid rotation (i.e. uniform $\Omega\approx j_{\rm rot}/(2/3)r^{2}$ – note the weak dependence of the factor 2/3 on $\Omega$; Paxton et al. 2019). AM transport arising from magnetic field amplification by Taylor instability (Spruit, 2002; Fuller et al., 2019) is efficient in radiative zones with largest shear in composition interfaces. This realizes a nearly rigid rotation in the compositionally uniform H envelope, and discontinuities in the AM profile near the edge of the He core at a mass coordinate of around $4~M_{\odot}$.

After mass accretion the merger product begins to contract, and sheds mass due to the large AM gained from the merger. The critical specific AM at the surface where centrifugal force balances with gravity is

[TABLE]

where $\Gamma_{\rm Edd}$ is the surface Eddington ratio which is well below unity ($0.1$–$0.5$) for our stellar models. As the star contracts to a BSG $j_{\rm crit}$ decreases to $\lesssim 10^{20}\ {\rm cm\ s^{-1}}$, resulting in the outer layers with $j>j_{\rm crit}$ lost as a (rotationally-enhanced) stellar wind. This results in a nearly instantaneous mass loss of typically $1$–$2~M_{\odot}$ in our models, and further mass loss by stellar winds brings the merger product to a sub-critical rotation. This kind of mass shedding accompanying envelope contraction has been suggested to be responsible for the equatorial circumstellar ring seen in SN 1987A (Morris & Podsiadlowski, 2007), although for our stellar models the entire process occurs (by construction) much earlier than suggested for the progenitor of SN 1987A, at $\sim 10^{6}$ years before core-collapse. A qualitatively similar phenomena is also seen in simulations of main-sequence merger products (Schneider et al., 2019).

For the $15~M_{\odot}+4~M_{\odot}$ model, the subsequent expansion to a RSG after core helium depletion changes the AM profile to a flatter one, reflecting the flattening in the density structure of the envelope as it becomes convective. The AM profile of the $15~M_{\odot}+8~M_{\odot}$ model is nearly unchanged till core-collapse as the envelope remains radiative.

While our prescriptions for AM evolution and transport do not capture the complicated three-dimensional processes that operate during the merger, our stellar models are found to be consistent with measured rotation rates of BSGs. The projected surface rotation velocities of observed BSGs in the LMC have a range of $v_{\rm surf}\sin i=(10$–$70)$ km s*-1* (Menon et al., 2024). Taking an angle-average as a representative value, $\langle\sin i\rangle=\frac{1}{4\pi}\int_{0}^{2\pi}d\phi\int_{0}^{\pi}\sin i(\sin i)di=\frac{\pi}{4},$ the resulting observationally estimated surface AM is

[TABLE]

consistent with our stellar models during core helium burning at an order of magnitude level.

3.2 Model Grid

In Figure 3 we plot the effective temperature $T_{\rm eff}$ of the full stellar models at two stages of evolution, in the middle of core helium burning (with central helium mass fraction of $0.5$), and in the end of our simulations at core carbon depletion. The former is more relevant for the supergiants we observe locally (e.g. in the LMC), while the latter is relevant for the transients upon core-collapse.

We find the trend on $M_{2}$ for BSGs described in Section 3.1.1 to be generally true across the adopted range of $M_{1}$. When looked at the core He burning phase, the threshold of $M_{2}$ for forming a BSG is lower for larger primary mass. This reflects the delayed expansion for more massive stars at LMC metallicities described in Section 3.1.1, with massive primaries only requiring a small amount of (or even zero) accretion to keep them blue during the helium burning phase. However, this trend is then reversed at core carbon depletion, where a larger $M_{2}$ is needed for forming BSGs for larger primary mass. With the exception of the $M_{1}=10~M_{\odot}$ models, the threshold for $M_{2}$ to form long-lived BSGs is roughly found as $M_{2}/M_{1}\gtrsim 0.5$, shown as a dotted line in the right panel of Figure 3.

We compare our models with recent parameter studies that follow similar approaches to evolve merger-produced BSGs (Schneider et al., 2024; Menon et al., 2024). Our model qualitatively agrees with these works that larger accretion ($M_{2}$) more likely leads to BSGs, but our result lies somewhere in the middle of Schneider et al. (2024) and Menon et al. (2024) in terms of how much $M_{2}$ is required to create long-lived BSGs. In the range of $M_{1}$ we consider, we produce more BSGs than Schneider et al. (2024) and less BSGs than Menon et al. (2024). A potential reason is that we consider a lower metallicity and higher $\alpha_{\rm sc}$ than Schneider et al. (2024), that both favor production of BSGs. On the other hand we neglect the dredge-up of He core and the hydrogen-burning shell during merger (as parameterized in e.g., Menon et al. 2024), which is important for reproducing the surface abundances of helium and CNO elements in observed BSGs. While the detailed level of dredge-up driven by merger is highly uncertain for these systems, including this effect will make the merger product bluer due to the increased mean molecular weight in the envelope.

Figure 4 shows the AM profiles of the full stellar models at core carbon depletion, with solid and dotted lines indicating BSG of $T_{\rm eff}>10^{3.9}$ K and red/yellow supergiants of $<10^{3.9}$ K respectively. The cooler supergiants have a flatter rotation profile across the envelope, while the BSGs have a sharp rise near the surface as expected for a radiative envelope with nearly rigid rotation (Section 3.1.2).

Comparing the BSG progenitors across different $M_{1}$, we find a lower specific AM for progenitors with higher $M_{1}$. This is shown in Figure 5 with the stellar models color-coded by the mass-averaged specific AM,

[TABLE]

that varies by a factor of several across our adopted range of $M_{1}$, from $3.5\times 10^{17}\ {\rm cm^{2}\ s^{-1}}$ for the BSG with $M_{1}=10~M_{\odot}$ to ($0.5$–$0.8)\times 10^{17}\ {\rm cm^{2}\ s^{-1}}$ for the models with $M_{1}=25~M_{\odot}$ (Table 1). We believe this arises mainly due to mass loss after these progenitors contract to BSGs. For a rigidly rotating envelope with gyration radius $kR_{*}$ and angular velocity $\Omega$, mass loss affects the rotation as

[TABLE]

Using $\langle j_{\rm rot}\rangle=\Omega(kR_{*})^{2}$, we obtain

[TABLE]

which can be solved, for initial values $\langle j_{\rm rot}\rangle=j_{0},M_{*}=M_{*,0}$, as

[TABLE]

Inspecting our BSG models (solid lines in Figure 4), the fractional mass lost during the BSG phase is similar among the progenitor models with a range $5$–$10~\%$. The key difference is the gyration factor $k$ measured in our models from $k^{2}=\langle j_{\rm rot}\rangle/(\Omega R_{*}^{2})$, ranging from $k^{2}\approx 0.03$ for the BSG model of $M_{1}=10M_{\odot}$ to $k^{2}\approx 0.01$ for BSGs of $M_{1}=25M_{\odot}$. As $k^{2}\ll 1$, this results in a huge difference in the power-law index in equation 8. While the above framework is approximate, this qualitatively explains the much slower rotation obtained for the stellar models with higher $M_{1}$. This variation with $M_{1}$ is also seen in models that end as RSGs (dotted lines in Figure 4), but with a stronger dependence aided by the higher mass loss rates of RSG winds compared to BSG winds.

The grey shaded region on Figure 4 is the threshold AM for a shell to circularize around a BH of mass $M$ and dimensionless spin $a_{\rm BH}=0$–$1$ (Fuller & Lu, 2022),

[TABLE]

where $c$ is the speed of light, and $r^{\prime}_{\rm ISCO}\equiv r_{\rm ISCO}/(GM/c^{2})$ ($=1$–$6$) is the innermost stable circular orbit of a BH scaled by $GM/c^{2}$,

[TABLE]

As we discuss in Section 4, comparison with $j_{\rm ISCO}$ is most relevant for cases where the progenitor fails to explode upon core-collapse and forms a BH. This $j_{\rm ISCO}$ serves as a benchmark separating the inner part of the star that falls into the newborn BH almost radially, and the outer part that can circularize around the BH to form an accretion disk. We find that the variations in the AM profile for different $M_{1}$ could result in huge changes in the amount of mass that can circularize around the newborn BH. For example, for the $M_{1}=10~M_{\odot}$ model around half of the star’s mass can have specific AM exceeding $j_{\rm ISCO}$, while for the $M_{1}=25~M_{\odot}$ models only the outermost $0.1$–$1~M_{\odot}$ of the star can have such large specific AM (see also Table 1). Such difference in mass have important consequences in the resulting high-energy transients, as we discuss in Section 4.2.

Since the mass loss affects the final AM, the degree of rotation depends on the (unfortunately uncertain) wind prescription. We check this by running additional calculations for the fiducial $15M_{\odot}+8M_{\odot}$ BSG model, with wind efficiency factor varied by factor 2 on both sides (i.e. Dutch_scaling_factor = 0.25, 1). The average AM $\langle j_{\rm rot}\rangle$ is found to mildly vary by $\approx 30\%$ on both sides (Appendix B), and other surface properties vary only by a few %. We thus expect our qualitative picture of the model grid to still hold.

This additionally indicates a metallicity dependence of our finding, where lower metallicities leading to weaker winds would be more favorable to produce rapidly rotating BSGs. However, it is also noteworthy from our models that the stellar merger channel does not require extremely low metallicities to produce fast-rotating BSGs. Previous stellar modeling, mainly in the context of ultra-long GRBs, have adopted much lower metallicities than ours when producing rotating BSGs ($Z=0$–$2\times 10^{-4}$; Kashiyama et al. 2013; Perna et al. 2018; Ror et al. 2024). Such extremely low metallicities, which are rather unrealistic in the redshifts these events are observed (e.g., Figure 6 of Neijssel et al., 2019), are not a requirement for the stellar merger channel.

4 Transients upon Core-collapse

Here we consider the expected transients from the core-collapse of these RSGs and BSGs from merger products. Besides the diversity due to the progenitors, we further expect a more crucial dichotomy in the transient signatures from whether the neutrino-driven explosion at the inner part of the star is successful or not. We respectively discuss the two cases where the neutrino heating drives a successful explosion (Section 4.1), and where the neutrino-driven explosion fails and a BH forms in the center (Section 4.2).

As a whole, we expect most of our stellar models to successfully explode by neutrino heating, while some of them fail to explode by neutrinos. Which stars do one or the other is more difficult to predict, as it requires core-collapse SN simulations with accurate initial conditions for the innermost core structure. Absent of first-principle simulations for neutrino-driven explosions of merger products, we here take a model-agnostic approach and include all of our progenitors as candidates for either a successful or failed explosion. Future numerical works with the aforementioned capabilities would be able to narrow down the transient landscape further (see also discussion in Section 5).

4.1 The Case of Successful Neutrino-driven Explosions: From SN II-P to 1987A-like SNe

4.1.1 Supernova Light Curves

It is well known that the light curves from explosions of RSGs and BSGs greatly differ, due to differences in their radii. The former, typically a Type II-P SNe, shows a characteristic plateau in the light curve powered by diffusion of internal energy deposited by the SN shock (Litvinova & Nadezhin, 1983; Popov, 1993; Kasen & Woosley, 2009). In the case of more compact BSG progenitors, the internal energy from the SN shock is quickly lost by adiabatic expansion, and the main powering source of the SN is the radioactive decay of 56Ni and 56Co. The interplay of radioactive heating and the gradual recombination of the ejecta creates a light curve rising slowly over months, as seen in SN 1987A and reproduced in a number of light curve models (e.g., Arnett & Fu, 1989; Shigeyama & Nomoto, 1990; Utrobin, 1993; Blinnikov et al., 2000; Moriya & Menon, 2024; Pumo & Cosentino, 2025).

To demonstrate this for our stellar models, we calculate bolometric light curves from explosions of our pre-SN models using a 1D radiation hydrodynamical code SNEC (Morozova et al., 2015). Assuming spherical symmetry as initial conditions is a good approximation even for our rotating models, as the stars are rotating at speeds much lower than their mass-shedding limits, and centrifugal forces that deform the star are weak. While multi-D models have demonstrated that the early emission around shock breakout is affected by asymmetries in the near-surface stellar structure or in the explosions themselves (Goldberg et al., 2022a; Irwin et al., 2021; Vartanyan et al., 2025), we expect the later emission at the main peak or plateau to be more robust to the effects of asymmetry.

For all the simulated explosions, we excise the innermost $1.4~M_{\odot}$ of the pre-SN models assuming it becomes a neutron star. We then adopt SNEC’s thermal bomb formulation to inject internal energy over 1 sec in the innermost $0.1~M_{\odot}$ of the leftover star, with energy set such that the final energy of the ejecta becomes a user-specified value $E_{\rm ej}$. The computational region is sampled by 3000 Lagrangian cells, with finer sampling near the center and the surface.

We fix the mass of 56Ni in the ejecta as $0.07~M_{\odot}$, a value inferred in SN 1987A from the 56Co-powered tail. The degree of mixing of 56Ni in the ejecta affects the light curve morphology and duration, particularly for BSG explosions. While mixing itself is a multi-D phenomenon, some 1D codes like SNEC allow users to parameterize mixing of 56Ni by hand. We set 56Ni to be efficiently mixed in the ejecta out to $90\%$ in mass, motivated by 3D simulations on core-collapse SNe from both BSGs and RSGs (e.g., Hammer et al., 2010; Utrobin et al., 2015; Stockinger et al., 2020; Vartanyan et al., 2025).

The left panel of Figure 6 compares the explosions of four merger models with $M_{1}=15~M_{\odot}$ and $E_{\rm ej}=10^{51}$ erg, with two stars ($M_{2}=4,6~M_{\odot}$) being RSGs and the other two stars ($M_{2}=8,12~M_{\odot}$) being BSGs. There is a clear difference in the light curves of these two kinds of explosions over the first $\sim 100$ days, with a IIP-like morphology in the two RSG models and a slowly rising SN 1987A-like morphology for the two BSG models. For a fixed $E_{\rm ej}$, the more massive merger products produce a slightly dimmer light curve with longer duration, due to the larger ejecta mass resulting in a larger photon diffusion time.

The right panel of Figure 6 shows three BSG explosion models with varied progenitor and explosion energy, and compared to the bolometric light curves of SN 1987A (Hamuy et al., 1988), SN 2005ci (Taddia et al., 2016), and SN 2009mw (Takáts et al., 2016), which all have comparable luminosities at tail such that the required mass of 56Ni are similar. While variations in the aforementioned parameters slightly change the onset and morphology of the main peak, the overall light curve evolution of (i) a fast initial drop from shock cooling emission of the outermost ejecta, (ii) a slow rise to main peak, and (iii) fast fall to the 56Co-powered tail due to hydrogen recombination, is in good qualitative agreement to the observed events. An additional parameter not captured in the SNEC calculations is the mixing of hydrogen in the ejecta, which can affect the duration of the main peak governed by hydrogen recombination. Modeling of light curve and nebular spectra for SN 1987A imply some mixing of hydrogen in the inner core (Shigeyama & Nomoto, 1990; Utrobin, 1993; Kozma & Fransson, 1998), but detailed tuning to match the light curve of SN 1987A is beyond the scope of this work (see e.g. Utrobin et al. 2015; Menon et al. 2019 for such attempts).

4.1.2 Event rates of 1987A-like SNe

The rate of 1987A-like SNe from BSG progenitors is suggested to be 1–3 % of core-collapse SNe (Smartt et al., 2009; Kleiser et al., 2011; Pastorello et al., 2012; Taddia et al., 2016; Sit et al., 2023), with a preference for sub-solar metallicities comparable to LMC (Taddia et al., 2013; Sit et al., 2023). As we cover the mass range of progenitors where we expect a large fraction of them to explode, the observed rates is a good cross check of our binary models. To estimate the event rate of these SNe from our binary channel, we conduct a similar exercise to Justham et al. (2014).

From the results in Section 3.2, we infer that binaries with $q\gtrsim 0.5$–$0.6$ end up as BSGs upon core-collapse. Our MESA models do not constrain the initial separation of the binaries required to produce BSGs, but this can be roughly predicted from binary evolution. At large separations of $a\gtrsim$ (a few) $10^{3}~R_{\odot}$, the two stars do not interact, and evolve as single stars that become RSGs at core-collapse (or Wolf-Rayet stars for ZAMS mass $\gtrsim 35~M_{\odot}$; Sukhbold et al. 2016). At intermediate separations of $a\approx$ (a few) $100$– (a few) $1000~R_{\odot}$ the binary can interact as the primary crosses the Hertzsprung gap, but a clear core-envelope boundary would exist in the primary. The primary’s envelope is loosely bound such that it is likely to be ejected in the common envelope phase (e.g., Klencki et al., 2021) at such high values of $q$, making the primary a hydrogen-poor star. At the tightest separations of $a\lesssim 10~R_{\odot}$ the binaries merge earlier in the main sequence (e.g., de Mink et al., 2013; Kinugawa et al., 2024), and rejuvenates as a single main-sequence star that most likely ends its life as a RSG or a Wolf-Rayet star (Schneider et al., 2024). We therefore expect (early) Case B mergers leading to BSGs to occur in binaries with initial separations of $a\approx 10$ – (a few) $100~R_{\odot}$. For a log-uniform distribution of $a$ spanning 5–6 dex commonly adopted in population synthesis modeling, binaries in this separation range occupy $f_{a}\sim 20$–$30\%$ of the entire binary population.

Further adopting a binary fraction of $f_{\rm bin}\sim 50$% for massive stars in LMC-like environments (Sana et al., 2013; Dunstall et al., 2015) and a flat mass ratio distribution (Shenar et al., 2022) with binaries of $0.5(0.6)\leq q\leq 0.8$ occupying $f_{q}\approx 30\%(20\%)$, we obtain the core-collapses of BSGs from our merger channel to be $f_{a}f_{\rm bin}f_{q}\approx$ $2$–$4\%$ of all core-collapses. Our estimate reasonably agrees with the observed fraction of 1987A-like SNe, given the expectation that a majority of our BSGs explode. This implies that binaries in sub-solar metallicity environments play a key role in producing 1987A-like SNe.

As the threshold for $q$ to undergo a contact phase is somewhat uncertain, the rate would be reduced if one considers a lower threshold of $q$ for a contact phase to occur. Meanwhile we expect an additional contribution from wider, low-$q$ binaries undergoing Case C mergers (as proposed for SN 1987A), although it requires the separation to be further tuned to a narrow range for unstable Case C mass transfer to occur while avoiding unstable Case B mass transfer.

4.2 The Case of Failed Neutrino-driven Explosions: Accretion-driven Transients

If the neutrino-driven explosion fails and a BH forms, the fate of the leftover star is completely different from a successful SNe. The loss of pressure support in the center would lead to infall of the inner part of the star onto the newborn BH, while mass loss due to neutrino emission preceding BH formation could eject the outermost part of the star (Nadezhin, 1980; Lovegrove & Woosley, 2013; Fernández et al., 2018; Ivanov & Fernández, 2021). In our progenitors with rotating envelopes, we may expect an energetic transient powered by accretion of stellar material onto the BH (e.g., Dexter & Kasen, 2013; Kashiyama & Quataert, 2015; Moriya et al., 2018).

For RSGs however, we expect that a large fraction (or all) of the outer envelope is ejected due to neutrino mass loss (Lovegrove & Woosley, 2013; Fernández et al., 2018; Ivanov & Fernández, 2021). Even if the inner envelope could fall onto the BH, the AM from rotation ($\lesssim 10^{17}$–$10^{18}$ cm2 s*-1*; Figure 4) is subdominant with respect to that expected from random convective motion ($10^{18}$–$10^{19}$ cm2 s*-1*; Goldberg et al. 2022b; Antoni & Quataert 2022). Therefore coherent circularization of the envelope in the form of an accretion disk is likely not realized for RSG progenitors, and it has instead been suggested that the randomly-oriented infall will unbind the envelope via a weak explosion of energies $10^{48}$–$10^{49}$ ergs (Antoni & Quataert, 2023). The resulting transient has been investigated in past works, and is expected to have a shock breakout emission of luminosity $10^{7}$–$10^{8}L_{\odot}$ that lasts for days to weeks (Piro, 2013; Lovegrove et al., 2017; Tsuna et al., 2025), followed by a months to year-long plateau emission with luminosity $10^{6}$–$10^{7}L_{\odot}$ (Lovegrove & Woosley, 2013; Antoni & Quataert, 2023).

On the other hand, BSGs have compact envelopes that results in weak mass ejection ($<0.1~M_{\odot}$) by the neutrino mass loss (Fernández et al., 2018; Ivanov & Fernández, 2021). As the envelope is radiative and rotation dominates its AM content, coherent circularization of the envelope would be realized. The envelope falls back over the free-fall time at the stellar surface

[TABLE]

We find in Figure 4 that the material able to circularize outside the BH has masses of $M_{\rm circ}\sim 0.1$–$10~M_{\odot}$, local free-fall times of $10^{3}$–$10^{5}$ sec, and bulk specific AM of $j\sim 10^{17}$–$10^{18}$ cm2 s*-1*. The circularization of fallback matter forms an accretion disk at a characteristic radius

[TABLE]

and accretes onto the BH at a viscous timescale of

[TABLE]

where $\alpha$ is the viscosity parameter (Shakura & Sunyaev, 1973), and $H/R$ is the height-to-radius ratio. As $t_{\rm visc}\ll t_{\rm ff}$, the accretion rate is regulated not by viscosity but by the infall of the envelope, with accretion rates of $\dot{M}_{\rm acc}\sim M_{\rm circ}/t_{\rm ff}\sim 10^{-5}$–$10^{-3}\ M_{\odot}\ {\rm s^{-1}}$. This is orders of magnitude higher than the Eddington-limited accretion rate where photons can cool the accretion flow, yet lower than the rate where neutrinos can cool the flow ($\gtrsim 10^{-2}\ M_{\odot}\ {\rm s}^{-1}$; Chen & Beloborodov 2007; Kawanaka et al. 2013). For such radiatively inefficient accretion flows, we expect outflows from the disk (Shakura & Sunyaev, 1973; Blandford & Begelman, 1999), with high velocities of the order of the escape velocity within the disk $v_{\rm wind}\gtrsim\sqrt{GM_{\rm BH}/r_{\rm disk}}\sim 0.1c$ (see also Strubbe & Quataert, 2009; Kashiyama & Quataert, 2015).

A subset of the BSG models with the largest accretion rates may further launch relativistic jets, analogous to the collapsar scenario for long GRBs (MacFadyen & Woosley, 1999). Such “BSG collapsars” have previously been proposed as channels for ultra-long GRBs (e.g., Woosley & Heger, 2012; Gendre et al., 2013; Kashiyama et al., 2013; Nakauchi et al., 2013; Perna et al., 2018), as they can naturally realize engine durations ($\sim t_{\rm ff}$) consistent with the observed durations of these GRBs. While there are only a handful of observed ultra-long GRBs, one event (GRB 111209A) has a superluminous SN-like optical counterpart (SN 2011kl), potentially consistent with this scenario (Levan et al., 2014; Greiner et al., 2015).

4.2.1 BSG Collapsars: Envelope Accretion

To investigate the expected transients in detail, we model the time-dependent accretion of BSG collapsars using the model of Fuller & Lu (2022) (their Appendix A). The disk is evolved by a one-zone approach that includes infall of stellar material, viscous accretion, and the wind outflow. The disk forms when the specific AM $j$ of the infalling shell of material first exceeds the $j_{\rm ISCO}$ of the BH, for a BH mass and spin given by assuming all the material interior to this shell collapses into a BH. The disk material is then continuously sourced by infall of stellar material with $j>j_{\rm ISCO}$, and is taken away by both accretion onto the BH and the disk wind.

At any given time the disk is characterized by two parameters, its mass $M_{d}$ and total AM $J_{d}$. The viscous evolution of the disk is calculated for $\alpha(H/R)^{2}=0.03$. In our BSG models where the accretion rate is regulated by the free-fall time of the stellar material rather than the viscous time, the accretion is insensitive to this choice.

Radiatively inefficient accretion flows are expected to be geometrically thick, with aspect ratios of $H/R\sim 0.3$–$1$, and are prone to powerful outflows. The model assumes that once the disk forms, only the equatorial material within a given solid angle of the disk manage to fall back onto the disk, and polar material outside this solid angle is removed by mechanical feedback from the disk wind/jet. We adopt a constant disk opening angle of $\theta_{\rm disk}=45$ degrees, corresponding to $H/R=1$. This results in a fraction $\sin(\theta_{\rm disk})\approx 0.7$ of stellar material with $j>j_{\rm ISCO}$ contributing to the disk.

The outermost bit of the envelope is expected to be ejected by the aforementioned neutrino mass loss, and will not contribute to the infalling mass. The ejected mass most sensitively depends on the compactness of the envelope $\xi_{\rm env}=(M_{*}/M_{\odot})/(R_{*}/R_{\odot})$ as defined in Fernández et al. (2018), which is expected since the binding energy of an outermost shell of a given mass $\Delta M(\ll M_{*})$ scales with $M_{*}/R_{*}$. For each stellar model we adopt an ejected mass inversely scaling with $\xi_{\rm env}$ as $M_{\rm ej}\approx 0.01\xi_{\rm env}^{-1}\ M_{\odot}$ based on the simulations of Ivanov & Fernández (2021), which results in $M_{\rm ej}\approx 0.01$–$0.1~M_{\odot}$ for our BSG models.

With these prescriptions, we solve for the time-evolution of the disk parameters using eq (A1)–(A4) of Fuller & Lu (2022). The radial mass-accretion rate in the disk down to $r_{\rm ISCO}$ is expected to follow a power-law $\dot{M}(r)=(M_{d}/t_{\rm visc})(r/r_{d})^{s}$ (Blandford & Begelman, 1999), where $r$ is the radius from the BH and we set $s=0.5$ motivated from recent simulations of radiatively inefficient accretion flows (Cho et al., 2024; Guo et al., 2024). The time-dependent mass-loss rate and kinetic luminosity of the wind is then given by integration over $r$ as

[TABLE]

Finally the accretion onto the BH is given by

[TABLE]

Figure 7 shows the accretion history onto the BH $\dot{M}_{\rm BH}(t)$ from disk formation, when the first (innermost) infalling material circularizes around the BH’s ISCO. We find $\dot{M}_{\rm BH}$ to be most sensitive to the primary mass, as expected from the large differences of the specific AM among these masses. The four $M_{1}=10,15\ M_{\odot}$ models (dash-dotted/dotted lines), the five $M_{1}=20\ M_{\odot}$ models (dashed lines), and the five $M_{1}=25\ M_{\odot}$ models (solid lines) respectively predict accretion rates of $10^{-4}$–$10^{-3}M_{\odot}\ {\rm s}^{-1}$, $10^{-5}$–$10^{-4}M_{\odot}\ {\rm s}^{-1}$, and $10^{-6}$–$10^{-5}M_{\odot}\ {\rm s}^{-1}$. The durations of the high accretion phase are in the range of $1$–$10$ hours, longer for stars with higher $M_{1}$ due to their larger radii (Table 1; equation 11).

4.2.2 BSG Collapsars: Implications for Ultra-long GRBs

Ultra-long GRBs are a class of GRBs with very long durations of hours and isotropic-equivalent luminosities of $\sim 10^{49}$–$10^{50}$ erg s*-1* (Levan et al., 2014). The $\gamma/$X-ray emission of the burst has non-thermal spectra and rapid variability similar to classical GRBs, indicating beamed emission from a relativistic jet.

In classical GRBs, the large energies in the relativistic jets have led to suggestions of them being generated electromagnetically (e.g., Lyutikov & Blandford, 2003; Kawanaka et al., 2013; Liu et al., 2015; Kumar & Zhang, 2015), where the jet is powered by the extraction of rotational energy from the BH assisted by large-scale magnetic fields (Blandford & Znajek, 1977). Figure 8 shows the BH spin of our BSGs models at the start of disk formation, calculated by the model in Section 4.2.1. The BH has moderate spin for all cases with $a_{\rm BH}\approx 0.2$–$0.7$, inherited from the finite AM of stellar material plunging into the BH. The rotational energies that can in principle be extracted from the BHs (e.g., equation 18 of Rees, 1984),

[TABLE]

are of the order of $10^{54}$ erg, sufficient to power a strong jet with energies observed in ultra-long GRBs.

Recent studies (Jacquemin-Ide et al., 2024; Lowell et al., 2024) have considered the BH spin evolution in collapsar disks, under the scenario that the magnetic field is strong enough such that the BH-disk system could reach the magnetically-arrested disk (MAD) state (Narayan et al., 2003; Tchekhovskoy et al., 2011). In this MAD regime, two key findings are: (i) The jet luminosity scales with the accretion rate with a spin-dependent efficiency $L_{\rm jet}=\eta_{\rm jet}(a_{\rm BH})\dot{M}c^{2}$, and (ii) The spin of the BH is quickly reduced to $a_{\rm BH}\approx 0.1$–$0.2$ even if only a small fraction of the original mass is accreted, which drops the jet efficiency to $\eta_{\rm jet}\sim 0.4$–$2\%$. More recent work (Wu et al., 2025) extended this framework to non-MAD disks, finding a maximal jet efficiency of $\eta_{\rm jet}\sim 2\%$ and lesser degree of spindown when the system deviates from MAD.

If we assume that the efficiencies for long GRBs similarly apply in our case of BSG collapsars (although with a few orders of magnitude lower accretion rate), $L_{\rm jet}$ could be calculated as shown in Figure 7 for $\eta_{\rm jet}=1\%$. The beaming-corrected average X-ray luminosity of $\gtrsim 10^{47}{\rm erg\ s^{-1}}$ inferred in ultra-long GRBs (Levan et al., 2014), or a jet power of $\gtrsim 5\times 10^{47}{\rm erg\ s^{-1}}$ for a radiative efficiency of $20\%$ (Beniamini et al., 2016), can be achieved for BSG models with $M_{1}=10,15,20~M_{\odot}$. These also meet the required duration of hours.

To observe the ultra-long GRB, the relativistic jet formed near the BH has to break out of the BSG envelope before the accretion ceases. To check this constraint, we calculate (see Appendix A) the propagation of a jet inside the star with a model based on Bromberg et al. (2011). Figure 9 shows the comparison between the time it takes for a jet to break out the BSG envelope $t_{\rm bo}$, and the lifetime of the jet $t_{\rm jet}$. The latter is set by a threshold on the jet power of $L_{\rm jet}>3\times 10^{47}$ erg s*-1*, but the dependence of $t_{\rm jet}$ on this is weak because of the sharp drop in accretion at late times. For all models we find $t_{\rm jet}>t_{\rm bo}$, i.e. the jet breaks out the star before it shuts off, and emission from the jet is observable for a duration $\approx t_{\rm jet}-t_{\rm bo}\sim 10^{3}$–$10^{4}$ sec.

While the jet propagates through the star, a fraction $\sim t_{\rm bo}/t_{\rm jet}\sim 10\%$ of the jet energy is deposited into the cocoon. We expect thermal emission from the nearly isotropic (outer) cocoon, especially visible for off-axis observers (e.g., Zheng et al., 2025). The cocoon’s mass from our jet model is found to range as $M_{c}\approx 0.1$–$0.6~M_{\odot}$, with asymptotic velocities $v_{c}\approx\sqrt{2E_{c}/M_{c}}$ approaching $\sim 0.1c$. The timescale of the outer cocoon emission is set by its diffusion time

[TABLE]

where we adopted a uniform opacity of $\kappa\approx 0.3\ {\rm cm^{2}\ g^{-1}}$. The cocoon adiabatically expands from an initial volume at breakout of $V_{c,0}=f_{c}R_{*}^{3}$, where $f_{c}\approx 0.08$–$0.13$ in our models, and radiates its remaining internal energy over a timescale of $t_{\rm oc}$ when it expands to a volume of $\sim 4\pi(t_{\rm oc}v_{c})^{3}/3\gg V_{c,0}$. This leads to a luminosity

[TABLE]

and effective temperature

[TABLE]

where $\sigma_{\rm SB}$ is the Stefan-Boltzmann constant. We thus predict an emission peaking in near-UV with AB magnitude of $\approx$ -17 mag at a wavelength of 250 nm, detectable out to $z\approx 0.2$ with ULTRASAT for a single-visit depth of 22.5 mag (Shvartzvald et al., 2024).

A final important constraint comes from the event rates of these ultra-long GRBs. While the number of detected ultra-long GRBs are much smaller than the classical long GRBs, ultra-long GRBs have lower luminosities and is visible on a smaller horizon. Taking this into account, Levan et al. (2014) inferred the event rate of ultra-long GRBs to be comparable or down to a factor several lower than long GRBs. The local long GRB rate is estimated to be $\sim 1\ {\rm Gpc}^{-3}\ {\rm yr}^{-1}$ (Wanderman & Piran, 2010; Lien et al., 2014), i.e. $\sim 0.1\%$ of all core-collapse SNe for a beaming fraction of $10^{-2}$, while the core-collapse of BSGs are estimated in Section 4.1.2 to be (a few) % of core-collapse SNe. We thus expect at most a few % of our BSG core-collapses will actually produce the ultra-long GRBs modeled here, regardless of whether the BSG collapsars are the progenitors of ultra-long GRBs (see e.g. MacLeod et al., 2014; Ioka et al., 2016; Marchant & Moriya, 2020; Hutchinson-Smith et al., 2024, for alternative scenarios). This upper limit is rather constraining, given that the rate of failed SNe could be as high as $\sim 20\%$ of massive star core collapse (Neustadt et al., 2021). This may indicate that a failed explosion and BH formation have a ZAMS mass distribution biased towards more massive BSGs, that would be less likely to have a massive rotating envelope to realize the extreme accretion.

4.2.3 BSG Collapsars: Wind-driven Luminous Transients

Regardless of whether a relativistic jet is launched, the mildly relativistic wind of $v_{w}\sim 0.1$–$0.3c$ from super-Eddington accretion could also power a bright thermal transient on its own. The wind will eject most of the equatorial disk because $r_{\rm disk}\gg r_{\rm ISCO}$, and the polar part of the star will also be unbound.

Figure 10 shows the kinetic luminosity of the disk wind $L_{\rm wind}$. The time-integrated energy of the wind varies with $4\times 10^{51}$ – $2\times 10^{53}$ erg, with the largest energies achieved by the $M_{1}=10,15~M_{\odot}$ models ejecting $M_{\rm ej}\sim 10~M_{\odot}$ of the star, and the least energetic being the $M_{1}=25~M_{\odot}$ models that eject $M_{\rm ej}\sim 0.2$–$0.7~M_{\odot}$.

We estimate the radiative output of this explosion below. As the wind breaks out of the star and the resulting ejecta expands, thermal emission is radiated over a diffusion timescale of

[TABLE]

where we have again used $\kappa\approx 0.3\ {\rm cm^{2}\ g^{-1}}$. As $t_{\rm diff}$ is longer than the duration of the disk wind $t_{\rm acc}\sim 10^{3}$–$10^{4}$ s, the ejecta radiates after a phase of adiabatic expansion, the energy supplied by the wind that occurs at a characteristic radius

[TABLE]

The former case of $r_{\rm inj}\approx R_{*}$ is valid for instantaneous energy injection (e.g. core-collapse SNe), while the latter case of $r_{\rm inj}\approx v_{\rm ej}t_{\rm acc}$ is valid for long-lasting energy injection (e.g. magnetar engine; Kasen & Bildsten 2010). Our case is marginally the latter, where $v_{\rm ej}t_{\rm acc}/R_{*}$ ranges from unity to a few. Since internal energy drops with radius as $\propto r^{-1}$ by adiabatic expansion, the resulting luminosity is

[TABLE]

The timescale, luminosity, and the photospheric velocity ($\sim v_{\rm ej}$) predicted from the low-$M_{1}$ models align with the “super-luminous SN” component seen in an ultra-long GRB 111209A (SN 2011kl; Greiner et al. 2015). This matches well with the picture raised in the previous section, of these low-mass BSGs also being more favorable progenitors for ultra-long GRBs.

On the other hand, the $M_{1}=25~M_{\odot}$ models predict a $\lesssim$ week-long emission with luminosity reaching $\sim 10^{44}$ erg s*-1*, overlapping with the parameter space of AT2018cow-like fast blue optical transients (FBOTs; Ho et al., 2023). As initially proposed for AT2018cow (e.g., Margutti et al., 2019), the BSG collapsar model could explain many of the key observational features of FBOTs, such as the presence of hydrogen and helium, fast and asymmetric ejecta, low 56Ni yield, and the star-forming host environment.

Still, the current model faces two potential challenges when interpreting the multi-wavelength phenomenology of FBOTs, which merit further investigations. First is the X-ray emission lasting for up to years (e.g., Coppejans et al., 2020; Yao et al., 2022; Migliori et al., 2024), which suggests long-lived accretion onto the remnant compact object. While our model does not predict strong accretion over such timescales (Figure 7), a key component not included in our model is fallback material of the neutrino-driven weak explosion (Fernández et al., 2018). This can create a power-law tail in the accretion rate, that can keep the accretion super-Eddington for years after core-collapse. This material is guaranteed to have a large AM of $10^{18}$–$10^{19}\ {\rm cm^{2}\ s^{-1}}$, enough to sustain an accretion disk.

Second is the existence of fairly dense circumstellar matter inferred from radio observations (e.g., Ho et al., 2019; Margutti et al., 2019; Coppejans et al., 2020; Ho et al., 2020, 2022; Yao et al., 2022). The inferred mass-loss rates at $\sim 10^{17}$ cm are $\sim 10^{-5}$–$10^{-4}~M_{\odot}\ {\rm yr}^{-1}$, when scaled for a BSG-like wind velocity of $100$ km s*-1*. The mass-loss rates predicted by the prescriptions in our MESA models at the last centuries of evolution are typically $\sim 10^{-6}M_{\odot}$ yr*-1*, at least an order of magnitude lower. One possibility to resolve this discrepancy is that the mass loss is underestimated, as the default Dutch wind prescription does not include LBV-like mass loss. Inspecting the stellar tracks of our $M_{1}=25M_{\odot}$ BSG models, we find that they typically enter the empirically-known LBV instability region (Smith et al., 2004, see Figure 5) when they slightly expand at helium shell burning. The measured mass loss rates for LBV winds is $\sim 10^{-5}$–$10^{-4}~M_{\odot}\ {\rm yr}^{-1}$ (Smith, 2014), which matches with the values inferred from radio. However, such extreme mass-loss over the final $\sim 10^{3}$–$10^{4}$ yrs of the star from helium shell burning to core-collapse can be problematic, as it could entirely remove the rotating outermost envelope of the star, required for triggering accretion that powers the FBOTs themselves.

5 Conclusion

We have presented rotating supergiant models from post-main sequence stellar mergers, using the stellar evolution code MESA with updated prescriptions for AM transport inside the star. In our model, the key parameter that set the dichotomy of a RSG or BSG at death was the mass accreted onto the primary, with larger binary mass ratios (typically $\gtrsim 0.5$–$0.6$) more likely leading to BSGs. Though the merger product loses some mass and AM due to centrifugal mass loss, our BSG models tend to maintain fast rotation in the outer envelope until core-collapse, and models with smaller ZAMS masses lead to faster rotation due to the effects of stellar winds during the helium burning phase and beyond.

Using the models evolved until close to core-collapse we made predictions for the transients of these stars upon core-collapse, considering both cases where the canonical neutrino-driven explosion succeeds and fails. For the former case of canonical SNe, light curve modeling successfully reproduced the distinct SN light curves from RSGs and BSGs, i.e. the Type II-P and the Type II-pec (or 1987A-like) SNe.

The failed SN case is found to have a richer phenomenology, in particular for BSG progenitors due to it being capable of forming an accretion disk around the newborn BH. Using a model for time-dependent accretion of the rotating envelope onto the BH, we find extreme accretion rates of $\sim 10^{-6}$–$10^{-3}M_{\odot}$ s*-1* over hours, with higher rates achieved for BSGs with lower ZAMS mass that maintain faster rotation. We considered the observational implications of the relativistic jet and the mildly relativistic disk wind, expected for such extreme accretion rates in the context of the collapsar scenario in GRBs. Our modeling points to these merger-origin BSGs as viable progenitors of ultra-long GRBs, as well as fast luminous transients like AT2018cow. A schematic flowchart for the transient landscape derived in this study is depicted in Figure 11.

We conclude by raising suggestions for future improvements and applications to our modeling. First, our prescription for merger of rapid accretion onto the primary is admittedly a simplification of the various physics involved in the merger process. This most importantly neglects the details of the secondary’s structure. When the mass ratio of the binary is closer to unity, the secondary is evolved in the main sequence. The secondary can then have a helium-rich core with low entropy that can preferentially sink towards the core, affecting both the core and the envelope helium abundance of the merger product. An enhanced helium abundance in the envelope leads to hotter envelopes, potentially increasing the parameter space that lead to BSGs. Recent work comparing various merger prescriptions of a single binary merger model (Patton et al., 2025) find differences in the core structures of the “merger” products, while their surface properties after thermal relaxation are qualitatively similar. However the outcome should depend on the parameters governing the merger such as the mass ratio and the binary separation, and a future work of constructing a model grid with these parameters would be important.

In this work we refrained from detailed assumptions of which of the progenitors end their lives as NSs and BHs, which remains an open topic of research. Recent simulations (e.g., Burrows et al., 2023, Figure 15) predict an “island of BH formation” of progenitors with He core masses of $3$–$4~M_{\odot}$ (corresponding to ZAMS masses of $12$–$15~M_{\odot}$) in their stellar models, and possibly at higher ZAMS masses around $20~M_{\odot}$ (e.g., Boccioli & Fragione, 2024). For our model that includes convective overshoot the corresponding ZAMS mass may be slightly lower, while additional variations on explodability is expected in the merger models that experience accretion (e.g., Justham et al., 2014; Schneider et al., 2024). It would be interesting in future work to run our models further up to core-collapse to obtain detailed predictions of the innermost core, which could give insights on the explodability of these stars through e.g. their profiles at the silicon-oxygen interface (e.g., Wang et al., 2022; Boccioli et al., 2023). This requires accurately modeling (and resolving) the innermost core structures, by running with a more sophisticated nuclear reaction network with a much larger number of isotopes than that adopted in this study (Farmer et al., 2016; Renzo et al., 2024)333The choice of nuclear networks also affected disk formation for the collapsar models considered in Renzo et al. (2024). For our BSG models at central carbon depletion, the innermost regions do not have enough AM for disk formation, and the outer edge of the carbon-oxygen core has an AM transport timescale of $r^{2}/\nu_{\rm AM}\sim 10^{4}$ yrs, much longer than the time left till core-collapse. We thus expect that evolving the BSG models further would not alter our predictions for the mass accretion outside the newborn BH..

For the BH spins obtained by the one-zone accretion disk model, we generally see an anti-correlation with the final BH mass (Figure 8). While the merger products from field binaries would lead to an isolated single BH, those in dynamical environments can further interact and merge with other compact objects. Massive star mergers can occur frequently in the centers of clusters, if they gather in the center upon formation by (primordial) mass segregation. They are potentially important for producing massive mass-gap BHs (e.g., Portegies Zwart et al., 2004; Kremer et al., 2020; González et al., 2021) and spinning BHs observed in some gravitational-wave events, which are generally difficult to produce by multiple generations of BH mergers due to the large recoil kicks expected during the hierarchical merger process. Investigations of the final masses and spins of stellar merger products, using an expanded grid of mass and metallicity, would be an important piece to better understand such channels for very massive and/or spinning BHs.

Acknowledgements

D.T. is supported by the Sherman Fairchild Postdoctoral Fellowship at Caltech. We thank Samantha Wu, Ryan Chornock, Kaustav Das, Kazumi Kashiyama, Kyle Kremer, Anthony Piro, Toshikazu Shigeyama, Mathieu Renzo, and David Vartanyan for valuable discussions. Finally, we deeply thank Bill Paxton for his continuous developments of the MESA code, without which this work and many other of our works have not been possible – the MESA code will live on and continue to shine.

Appendix A Jet propagation in BSG Envelope

A semi-analytical model for the interaction of a relativistic jet with the surrounding stellar envelope was provided by Bromberg et al. (2011). Here we adopt this model for the propagation of the accretion-powered jet through the BSG envelope, to estimate when/if the jet can break out from the star. The jet power is set as

[TABLE]

where $\eta_{\rm jet}=0.01$ (see Section 4.2.2) is the jet efficiency, and here a factor $1/2$ is included as we consider only one side of the jet (towards the observer). A top-hat jet with opening angle of $\theta_{0}$ at launch, drills through the stellar material with a density profile $\rho_{*}(r)$. We set $\theta_{0}=0.3$ rad following Perna et al. (2018), which was motivated by values constrained from ultra-long GRBs (e.g., Levan et al., 2014; Kann et al., 2018).

We assume the jet launched at $t=0$ travels unimpeded (with velocity $c$) out to the innermost radius $r_{\rm in,circ}$ where the fallback material can circularize, and starts interacting with the BSG material at time $t=r_{\rm in,circ}/c$ from launch. The initial radius of the jet head and cocoon is set as $r_{h,0}=r_{\rm in,circ}$ and $r_{c,0}=(\tan\theta_{0})r_{\rm in,circ}$ respectively. For simplicity, we assume the BSG envelope as static and neglect the potential infall of the BSG envelope exterior to $r_{\rm in,circ}$.

At a given time, the velocity of the jet’s head is set by pressure balance as (Matzner, 2003)

[TABLE]

where $\tilde{L}$ is a dimensionless parameter

[TABLE]

$N_{s}$ is an order-unity factor derived from calibration to simulations (eq 13 of Harrison et al., 2018), and $\Sigma_{\rm jet}$ is the jet’s cross section set by whether the jet is collimated by the cocoon’s pressure $P_{c}$ as derived below. The jet continuously dissipates its energy into the cocoon, and neglecting adiabatic losses (see Appendix of Bromberg et al., 2011) the cocoon’s energy (focusing only on the observer’s side) is

[TABLE]

Due to its pressure the cocoon expands sideways with velocity

[TABLE]

where $\bar{\rho}_{*}$ is the mean density in the volume occupied by the cocoon

[TABLE]

where $V_{c}$ and $M_{c}$ are respectively the volume and mass of the cocoon. The pressure of the cocoon is dominated by radiation

[TABLE]

The collimation shock by the cocoon splits the jet into an unshocked region and shocked region. The unshocked jet’s height is given as a function of $L_{\rm jet}$ and $P_{c}$ as (Bromberg et al., 2011, their Figure 2)

[TABLE]

with the jet being uncollimated below $\approx\hat{z}/2$ and collimated above. Comparison of $\hat{z}/2$ and $r_{h}$ gives the criterion of collimation at the jet head, and the jet cross section at the head is then

[TABLE]

The jet breaks out the star when the jet head $r_{h}$ reaches the stellar radius $R_{*}$, and we define the time this occurs as the breakout time $t_{\rm bo}$ in Section 4.2.2. We note that the model is not strictly applicable close to the stellar surface, where the density profile $\rho_{*}(r)$ sharply drops (faster than $r^{-3}$). The forward shock accelerates and loses causal contact with the jet, and the assumption of pressure balance can break down (Bromberg et al., 2011). Comparison of the analytical model with hydrodynamical simulations, for a GRB-like jet propagating through a Wolf-Rayet star, nevertheless shows reasonable agreement for the estimation of $t_{\rm bo}$ (Harrison et al., 2018).

Appendix B Model variations

To check how the results are sensitive to various prescriptions, we re-run our baseline BSG model of $M_{1}=15~M_{\odot}$, $M_{2}=8~M_{\odot}$ with variations in the input physics as follows:

we change the onset of mass accretion from a primary radius from our fiducial value of $50~R_{\odot}$ to [$30~R_{\odot},100~R_{\odot}$]. 2. 2.

we change the wind efficiency (Dutch_scaling_factor) that controls the mass loss rate, by a factor 2 on both sides from our fiducial value of $0.5$. 3. 3.

we update the AM transport scheme from the Fuller et al. (2019) formalism to the default prescription based on the Taylor-Spruit (TS) dynamo (Spruit, 2002). 4. 4.

we increase the mass resolution by a factor of 2, by lowering two parameters in MESA mesh_delta_coeff and mesh_delta_coeff_for_highT from our default 0.8 to 0.4.

Figure 12 shows the differences in four parameters of the pre-SN star most relevant for our work – mass, stellar radius, effective temperature, and mass-averaged specific AM ($M_{\rm fin},R_{\rm fin},T_{\rm eff,fin}$ and $\langle j_{\rm rot}\rangle$ in Table 1). We find the stellar parameters ($M_{\rm fin},R_{\rm fin},T_{\rm eff,fin}$) to weakly vary within $10\%$, suggesting our results and qualitative picture of the RSG/BSG dichotomy are not subject to these uncertainties. The largest variations are seen in the mass-averaged AM $\langle j_{\rm rot}\rangle$, with larger radii at merger increasing the AM and larger mass loss rate decreasing the AM, both expected from our framework of AM evolution (equations 2, 6). The effects of AM prescription are less clear, but the larger AM content for the TS dynamo case may be explained by a smaller AM loss by winds, due to less efficient AM transport to the surface. While they are important model uncertainties, we expect them to not change our qualitative picture of the transient landscape expected from these stars.

Finally, our parameters are unchanged within $5\%$ when the mass resolution is increased by a factor of 2, and we consider the models to be sufficiently converged in resolution for our purpose.

Appendix C Summary of Model Grid

In Table 1 we show our model grid and the resulting stellar properties at core-collapse.