© The Authors 2025

1. Introduction

Gamma-ray bursts (GRBs) are some of the most energetic events in the Universe, lasting between a fraction of a second and several hundreds of seconds. Powered by a central engine, jets form and propagate through surrounding material, emitting intense gamma rays. Their duration has become the defining feature in their classification into short and long GRBs, which is often linked to their formation pathway. Several lines of evidence (Berger 2014) consistently point to short GRBs (≲2 seconds) being produced after mergers of binary neutron stars, in line with the early suggestion by Eichler et al. (1989). The link has been verified in one case by the joint observations of the gravitational wave event GW170817 (Abbott et al. 2017c), the kilonova AT2017gfo/SSS17a (Abbott et al. 2017a), and the gamma-ray burst GRB170817A (Abbott et al. 2017b). Several long GRBs (LGRBs), on the other hand, have been associated with Type Ic broad-line supernovae (SN Ic-BL; Modjaz et al. 2016), indicating a hydrogen- and helium-poor progenitor massive star. The connection is thought to encompass most LGRBs, even though observations of some LGRBs have shown photometric and spectroscopic optical and near-infrared signatures compatible with kilonova emission (Rastinejad et al. 2022; Levan et al. 2024).

Due to their highly beamed and energetic nature, LGRBs have been observed up to z ∼ 9.4 (Cucchiara et al. 2011) and have been used to estimate star formation rates at high redshifts. However, the LGRB rate rises more rapidly over z ≤ 4 than other star formation estimators (Yüksel et al. 2008; Kistler et al. 2009; Ghirlanda & Salvaterra 2022). Furthermore, the local population of LGRB hosts (z ≤ 2) tends towards lower mass hosts compared to the general population of star-forming galaxies (Krühler et al. 2015; Vergani et al. 2015; Japelj et al. 2016; Perley et al. 2016). With more complete and unbiased samples of LGRBs and their host galaxies, Greiner et al. (2015) have shown that between z = 3 and z = 5, the LGRB rate traces the UV metrics of the cosmic star formation. Together with average metallicity measurements of LGRB host galaxies, this has been interpreted as a low-metallicity preference for the formation of LGRBs (Levesque et al. 2010a,b; Graham & Fruchter 2013; Perley et al. 2016; Vergani et al. 2017; Palmerio et al. 2019). A sharp cut-off of LGRB formation at a sub-solar metallicity can explain these observed phenomena over redshift, although its exact value ranges from 0.3 Z_⊙ to 0.7 Z_⊙, depending on the observational survey (Modjaz et al. 2008; Salvaterra et al. 2012; Graham & Fruchter 2013; Perley et al. 2016; Vergani et al. 2017; Palmerio et al. 2019) and metallicity calibration used (Kewley & Ellison 2008).

Despite the preference of LGRBs for sub-solar metallicity host galaxies, GRB 020819 (Levesque et al. 2010c), GRB 111005A (MichałowskI et al. 2018), and GRB 160804A (Heintz et al. 2018) have been observed to have solar or super-solar metallicity explosion sites. Furthermore, for z < 2.5, the fraction of near-solar metallicity LGRB host galaxies remains nearly constant at 10% to 20% (Vergani et al. 2015; Palmerio et al. 2019; Graham et al. 2023). A similar fraction of 20% of super-solar metallicity LGRB host galaxies has been observed by Krühler et al. (2015) with VLT X-SHOOTER emission-line spectroscopy. However, these fractions might be significantly lower due to observational biases. A metallicity measurement of the host galaxy requires a spectrum, which is easier to obtain for more massive, brighter hosts, which are also generally more metal-rich. Although the LGRB sample might be complete, the LGRB host sample with measured metallicities is not (for example, see Graham et al. 2023). Secondly, if a metallicity measurement is done, in most cases it only probes the average line-of-sight metallicity of the LGRB host due to their distance, but their metallicity distribution can be highly non-uniform (Metha & Trenti 2020). While the fraction of high-metallicity hosts is, thus, difficult to constrain, the direct detection of the three high-metallicity explosion sites provides evidence that at least some LGRBs originate from metal-rich environments.

The exact criterion for an LGRB to occur from the collapse of a massive star depends on the assumed central engine powering the GRB. The two best-studied scenarios involve either a highly magnetized, rapidly rotating proto-neutron star (a proto-magnetar) or a spinning black hole (BH) with an accretion disk (the standard ‘collapsar’ scenario). In the former scenario, the proto-magnetar powers the LGRB through its fast rotation and high magnetic field strength (Usov 1992; Thompson 1994; Thompson et al. 2004; Metzger et al. 2011; Mazzali et al. 2014). In the collapsar scenario a BH is instead formed, and sufficient angular momentum is available in the stellar envelope to form an accretion disk around it, allowing jets to form either by the Blandford & Znajek (1977) mechanism, or because of energy deposition by annihilations of neutrino-antineutrino pairs formed in the inner accretion disk (Eichler et al. 1989; Woosley 1993; MacFadyen & Woosley 1999). Independent of either central engine scenario, the progenitor star has to (I) reach iron core collapse stripped of hydrogen and helium, (II) contain sufficient angular momentum to either spin up the proto-magnetar or to form an accretion disk, and (III) have a successful and collimated jet breakout through the stellar structure and produce gamma-ray emission (Matzner 2003; Modjaz et al. 2016; Salafia et al. 2020). Combined, these three criteria require the stellar progenitor to be a rapidly rotating, highly stripped star, where the outer hydrogen and most of the helium envelopes have been removed. This outcome gives rise to two main formation pathway questions and conditions: (a) How does the stellar core retain or gain sufficient angular momentum at core collapse? (b) How does the star lose its envelope?

From a single-star perspective, an initially rapidly rotating zero-age main-sequence (ZAMS) star will expand to a red supergiant after finishing core hydrogen burning. The subsequent evolution of its rotation depends heavily on the assumed coupling between the core and outer layers of the star (Qin et al. 2019). Astroseismology observations of slightly less massive giant stars (Fuller et al. 2014; Cantiello et al. 2014; den Hartogh et al. 2020) indicate a strong internal coupling (Eggenberger et al. 2022; Moyano et al. 2023). As the star expands post-main sequence, this results in an efficient redistribution of angular momentum from the core to the outer stellar layers and leads to the formation of slowly spinning BHs at the collapse of the core (Qin et al. 2018; Fuller & Ma 2019). Strong stellar winds will further enhance this by removing mass and angular momentum from the outer layers. For a single star to form a rapidly rotating, highly stripped star, the metallicity has to be low to minimize mass and angular momentum loss (Z < 0.3Z_⊙), and the initial rotation needs to be sufficient for rotational mixing to keep the star close to chemical homogeneity (quasi-chemically homogeneous; Maeder 1987; Meynet & Maeder 2007). As a result of the strong mixing, the star burns nearly all its hydrogen and avoids the giant phase. Without an extended envelope, angular momentum loss from the core is minimized (Yoon & Langer 2005; Yoon et al. 2006; Woosley & Heger 2006). In a sub-solar metallicity environment, the Wolf-Rayet winds are sufficiently weak (Vink & De Koter 2005) for these quasi-chemically homogeneous stars to retain sufficient angular momentum in their core to power a LGRB.

Binary systems provide additional formation pathways for stars to gain angular momentum and be stripped during their evolution instead of requiring the star to be rotating at birth. On the other hand, processes such as common envelope evolution can remove angular momentum from the star and restrict the formation of LGRBs. In general, binary formation channels invoke chemically homogenous evolution or late-time tidal spin-up, although several other mechanisms have been proposed throughout the years (see Fryer et al. 2025, for an extensive overview). First, close massive binary systems can undergo mass transfer and spin up the accreting star sufficiently to undergo quasi-chemically homogeneous evolution at metallicities below Z = 1/3 Z_⊙ (Cantiello et al. 2007; Yoon et al. 2008, 2010; Eldridge et al. 2011; Ghodla et al. 2023), as long as the accretor has not yet formed a well-defined stellar core. Secondly, in very tight near-equal-mass-ratio contact binaries, tides can spin up both stars, where the chemically homogenous evolution keeps both stars from expanding while forming the LGRB progenitors (Marchant et al. 2016; Song et al. 2016; Mandel & de Mink 2016). Additionally, if the mass ratio is further from unity, just the initially more massive star in a binary can undergo this tidal spin up (de Mink et al. 2009; Marchant et al. 2017). Finally, tidal spin-up can also occur after a helium star has formed in tight binaries with a compact object companion (van den Heuvel & Yoon 2007; Detmers et al. 2008; Qin et al. 2018). These four LGRB formation channels are more common at sub-solar metallicities.

The tidal spin-up and the chemically homogeneous evolution formation channels have been investigated using several binary population synthesis codes. Izzard et al. (2004) and Chrimes et al. (2020) looked at the effect of spin-orbit coupling in LGRB formation using binary population synthesis, although the exact nature of systems experiencing the tidal effects was not explored. The latter includes the treatment of accretion-induced chemically homogeneous evolution. However, this was done using an approximate prescription (described in Eldridge et al. 2011), which was calibrated on detailed stellar models from Yoon et al. (2006). In the context of binary black hole mergers, Bavera et al. (2022) investigated LGRB formation through chemically homogenous evolution of close binaries and tidal spin-up in BH-helium star binaries. Similarly, Ghodla et al. (2023) showed the effect that tidally and accretion-induced chemically homogenous evolution have on the expected LGRB population.

Several additional scenarios for LGRB formation have been proposed over the years that generate the LGRB through other mechanisms, such as the merger between a helium star and a BH (Fryer & Woosley 1998) or an explosive common envelope ejection (Podsiadlowski et al. 2010). The contribution of the former to the cosmic LGRB rate has not yet been quantified, while the latter could produce ≃10⁻⁶ LGRBs yr⁻¹ at solar metallicity, but is expected to be more efficient at lower metallicities (Podsiadlowski et al. 2010). Finally, a merger between two helium cores has been suggested to generate a rapidly rotating core (Ivanova et al. 2003; Fryer & Heger 2005), but its population contribution is poorly studied.

The fact that most theoretical formation channels predict an inverse relation between LGRB progenitors and metallicity aligns well with the observational preference for sub-solar metallicity host galaxies. However, the three (super-)solar metallicity explosion sites and the apparent lack of evolution in the high-metallicity fraction of LGRB host galaxies across different redshifts (Graham et al. 2023) suggest additional formation mechanisms might be at play. In this work, we propose a previously unexplored formation channel through binary interaction that exhibits a positive relationship between metallicity and LGRB progenitors that could explain these effects. For binaries with an initial mass ratio close to unity, after a phase of stable mass transfer from the primary onto the companion, a second phase of stable mass transfer from the companion onto the now stripped primary star spins it up sufficiently to power an LGRB. We refer to this mass transfer from the initially less massive star (secondary) onto the initially more massive star (primary) as “reverse mass transfer”. We consistently refer to the initially more massive star as the primary and to the initially less massive star as the secondary or companion. We determine the occurrence of this LGRB formation channel using the next-generation binary population synthesis code, POSYDON, which uses grids of detailed single- and binary-star models, as described in Section 2. In Section 3, we explain the formation pathway in detail using an example model and show why the efficiency of this LGRB channel has a positive dependence on metallicity. Furthermore, we explore the available energy to power a jet in potential LGRB progenitors, and how different criteria for defining a successful LGRB affect the predicted rate density from the stable reverse-mass-transfer channel. We discuss uncertainties in this stable reverse-mass-transfer formation pathway and the rate density evolution in Section 4, while we summarize our conclusions in Section 5.

2. Method

We used the public state-of-the-art binary population synthesis code POSYDON¹ (Fragos et al. 2023; Andrews et al. 2025), which employs grids of detailed single- and binary-star models computed with the 1D stellar evolution code MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023) to create synthetic populations of evolved binary systems. POSYDON v2 release² contains grids at eight metallicities: 10⁻⁴, 10⁻³, 10⁻², 0.2, 0.45, 0.1, 1, 2 Z_⊙ with Z_⊙ = 0.0142, which include self-consistent modeling of internal rotation, as well as angular momentum transport in the stellar interior and between the stars and the orbit of the binary.

2.1. Stellar and binary evolution physics

We summarize the essential stellar and binary physics model assumptions in this section – see Fragos et al. (2023) and Andrews et al. (2025) for more detailed descriptions. Rotational mixing and angular momentum transport are implemented using the diffusion approximation and a Spruit–Tayler dynamo (Spruit 2002) is included, which results in strong coupling between the stellar core and its envelope during the post-main sequence, and an overall efficient angular momentum transport (Heger et al. 2005).

The POSYDON H-rich+H-rich binary star grids, which start their evolution as two ZAMS stars, combine different formulations for mass transfer depending on the evolutionary state of the stars when Roche lobe overflow (RLOF) occurs (for details, see Fragos et al. 2023). This grid self-consistently treats rotation, tides, and mass transfer within the binary, as the orbital evolution and stellar structure equations of both stars are solved with MESA. For stars on the main sequence, the contact scheme is used, which keeps the stars inside their Roche lobe (Marchant et al. 2016). The kolb scheme (Kolb & Ritter 1990) is used for stars undergoing RLOF post-main-sequence to allow for radial expansion past its Roche lobe since the envelope has a low density. The grids use the correction described in Xing et al. (2024) to allow mass transfer from the companion onto the primary star with the kolb scheme in MESA. Either star can initiate the mass transfer, and multiple mass-transfer phases can occur in the lifetime of a binary if the interaction remains stable. Stellar rotation enhances the stellar winds to keep rotation below the critical threshold ω_s/ω_s, crit ≤ 1. The latter is also a mechanism that can limit the accretion efficiency of a stable RLOF mass-transfer phase.

POSYDON considers mass transfer in the H-rich+H-rich grids as unstable when (i) a donor mass transfer rate of 0.1 M_⊙ yr⁻¹ is reached, (ii) when mass is lost from the L₂ point, or (iii) when both stars fill their Roche lobe while either has evolved off the main sequence³. If a contact phase happens while both stars are on the main sequence, this is modeled following Marchant et al. (2016) and can result in a stable overcontact phase. Binary MESA models reaching one of the unstable mass transfer criteria are stopped and evolved using POSYDON through a common envelope phase and subsequent detached or interaction phases. If none of the instability conditions are met, the binary evolution is tracked within MESA until the end of core carbon burning⁴, where the final fate of each star as a neutron star or BH is determined. At this point, we determine whether it is an LGRB progenitor.

2.2. Core collapse and long gamma-ray burst formation

In the collapsar scenario, the relativistic jet is powered by the accretion of material from the collapsing star onto the newly-formed BH through an accretion disk. We use the remnant mass prescription from Patton & Sukhbold (2020) to determine whether the core collapse produces a BH or a NS, assuming always a prompt collapse for BH formation. We adopt the methodology from Batta & Ramirez-Ruiz (2019) to track the formation and evolution of the accretion disk during the infall without wind feedback. In their approach, the accretion disk mass is determined by material in the stellar interior whose specific angular momentum exceeds that of the innermost stable circular orbit (ISCO) of the BH. Following Bavera et al. (2020), we split the stellar structure into shells of mass m_shell. We assume the innermost 2.5 M_⊙ of the star to collapse directly to form a BH. For each consecutive i-th shell, the fraction of the shell with sufficient angular momentum to form a disk is given by m_disk, i = m_shellcos(θ_disk, i), where

$$\begin{array}{c}\hfill {\theta }_{\mathrm{disk},i}=arcsin\left({\left[\frac{j_{\mathrm{ISCO},\mathrm{i}}}{j_{\max ,i}}\right]}^{1/2}\right),\end{array}$$(1)

with j_ISCO, i being the specific angular momentum of the inner-most stable orbit (ISCO), and j_max, i = Ω(r_i)r_i² that of material at the equator of the i-th mass shell, located at a stellar radius r_i. We use equation (5) from Batta & Ramirez-Ruiz (2019) for j_ISCO. The amount of angular momentum gained by the BH after accreting the i-th shell is ΔJ = m_disk, ij_ISCO. We perform the calculation per mass shell to account for the change in BH spin and the ISCO during the collapse.

We calculate the jet power assuming the Blandford-Znajek mechanism, which powers the GRB by extracting rotational energy from the BH (Blandford & Znajek 1977). We decompose the total jet luminosity using two efficiency factors, one based on the magnetic field and one on the BH spin, following Gottlieb et al. (2023, 2024):

$$\begin{array}{c}\hfill L_{\mathrm{BZ}}={\eta }_{\varphi }{\eta }_a\dot{M}{c}^2\mathrm{erg}{\mathrm{s}}^{-1},\end{array}$$(2)

where Ṁ is the accretion rate through the disk, η_a is an efficiency factor that depends on the BH spin, and η_ϕ is an additional efficiency factor that depends on the magnetic field flux at the BH horizon. The efficiency of extraction depends highly on the state of the disk and requires magnetohydrodynamic simulations to follow the full disk and jet formation, which is beyond the scope of this work, but we discuss the choices in more detail in Section 4. We followed the original spin dependence from Blandford & Znajek (1977) of L_BZ ∝ a² and set η_a = a², although see Tchekhovskoy et al. (2010, 2012) or Lowell et al. (2024) for a more complex dependence on the BH spin that also depends on the aspect ratio of the accretion disk. We track how this efficiency changes as the star collapses and the BH spin increases, but we do not model the extraction of angular momentum from the powering of the jet. We sum over the mass of each shell going to the accretion disk to get the total power of the jet available in the collapse:

$$\begin{array}{c}\hfill E_{\mathrm{BZ}}={\displaystyle \sum _i}{m}_{\mathrm{disk},i}{a}_{i-1}^2{c}^2,\end{array}$$(3)

where a_i is the BH spin parameter after accreting the i-th shell. While this does not track the change of efficiency from the evolving magnetic field flux at the BH horizon, different η_ϕ values should provide a reasonable approximation for the energy budget. We discuss the choices for the η_ϕ efficiency in Section 3.2.

2.3. Binary population properties

We create stellar populations at each metallicity in POSYDON following a Kroupa (2001) initial mass function between 7 M_⊙ and 120 M_⊙, with a flat mass ratio, and log-period distributions. We initiate 1 000 000 binaries per metallicity and evolve them using the POSYDON nearest neighbor interpolation scheme.

POSYDON implements a nearest neighbor and an initial-final interpolation method across its detailed model grids, which are evolved to either the end of core collapse or the onset of unstable mass transfer. The initial-final interpolator maps the initial properties of a system to its final state at these endpoints. Additionally, it includes interpolation across compact object formation to include specific properties, such as BH spin (for more details, see Fragos et al. 2023) However, here, we require a stellar profile to perform the collapse, as described in Section 2.2. To obtain this, we use the nearest neighbor scheme, where each system is matched to its closest MESA model in the mass-q-period space of the relevant grid. This provides a stellar profile at carbon exhaustion for binaries that have undergone stable mass transfer or no interaction⁵. We only select systems with stable reverse mass transfer and sufficient angular momentum in the pre-collapse structure to form an accretion disk during the collapse.

2.4. Star formation and metallicity history

We combine the binary populations with metallicity-dependent star formation history prescriptions to determine the occurrence rate of stable reverse-mass-transfer LGRB formation over cosmic time, where we assume a binary fraction of 0.7. The star formation rate and metallicity evolution are extracted from the TNG-100 IllustrisTNG simulation (Springel et al. 2018; Nelson et al. 2018; Pillepich et al. 2018; Naiman et al. 2018; Marinacci et al. 2018). This cosmological simulation contains baryonic and dark matter particles of masses 1.4 × 10⁶ M_⊙ and 7.5 × 10⁶ M_⊙, respectively, in a comoving volume of (75 Mpc/h)³.

The IllustrisTNG simulation is based on the moving-mesh code AREPO (Springel 2010), which combines Lagrangian particle-based methods and Eulerian mesh-based methods. Because the former methods tend to underestimate the mixing between cells (Wiersma et al. 2009), while the latter tend to overestimate it (Sarmento et al. 2017), the net effect on the metallicity of galaxies in the IllustrisTNG simulation is an active area of research. For these reasons, in Section 4 we implement two additional empirical star formation histories based on Madau & Fragos (2017) and Neijssel et al. (2019), for comparison. The total star formation rate from Neijssel et al. (2019) was calibrated to produce the observed binary black hole distribution⁶ with a log-normal metallicity distribution with σ = 0.35, $\mu ={log}_{10}(\langle Z\rangle )-{\sigma }^2\frac{ln(10)}{2}$ with $\langle Z\rangle =0.035\times {10}^{-0.23z}$, and Madau & Fragos (2017) has a higher star formation rate at high-redshift with $\mu =ln(\langle Z\rangle )-{\sigma }^2/2$ and $\langle Z\rangle ={10}^{0.153-0.074z^{1.34}}{\mathrm{Z}}_{\odot }$, which spreads out the metallicity-specific star formation rate compared to the Neijssel et al. (2019) relationship. We follow the parameters from Appendix A in Bavera et al. (2020), where they set σ = 0.5. With both these prescriptions, the cosmic metallicity evolution proceeds slower than in the Illustris-TNG simulation and contains significantly more star formation at very high redshift (z > 15). Adopting the cosmological parameter values from Planck Collaboration VI (2020), we convolve these star formation histories with the POSYDON populations to calculate a total stable reverse-mass-transfer LGRB rate.

3. Results

Since LGRB formation through stable reverse mass transfer is constrained by multiple stellar and binary evolutionary processes, we first examine a representative evolutionary model to outline the key factors at play. We then shift our attention to the LGRB population and how each process limits their formation in the proposed channel.

3.1. Stable reverse mass transfer towards LGRB formation

3.1.1. Example model

From the binaries in our population that undergo stable reverse mass transfer, here we take an example binary with sufficient angular momentum at the collapse of the primary star to produce an accretion disk around the BH. It consists of a primary and secondary star of initial masses M₁ ≈ 25.1 M_⊙ and M₂ ≈ 23.8 M_⊙, respectively, at solar metallicity, orbiting each other at an initial period of P ≈ 26.8 days. Figure 1 shows the time evolution of several properties of the binary, during a time window that begins close to the end of the main sequence phase. The binary’s evolution can be separated into five main phases:

Fig. 1. Typical example of the evolution of a potential LGRB progenitor through the stable reverse-mass-transfer channel. The initial properties of the binary are Z = Z⊙, M1, ZAMS ≈ 25.1 M⊙, q = 0.95, and P ≈ 26.8 days. In all panels, red (blue) solid lines refer to the primary (secondary) star. The columns are split into five different evolutionary phases, focusing on the mass transfer phases. The second and fourth columns show the mass transfer and ‘reverse’ mass transfer phases, respectively. The first, third and last columns show the long, nuclear timescale, detached evolution of the system. Top two rows: Total mass (solid lines), He core mass (dashed lines) and CO core mass (dotted lines). The He (CO) core boundaries are set where the hydrogen (helium) fraction falls below 0.1. Third row: Mass transfer rate from Roche lobe overflow. The periods of Roche lobe overflow are marked by a red (blue) shading when the donor is the primary (secondary). Fourth row: Evolution of stellar radius (red: primary; blue: secondary) and separation of the system (black). Fifth row: Surface rotational angular frequency in units of the critical (i.e., mass shed) value. Sixth row: the dimensionless spin parameter cJ/GM2 of the star.

1. Initially, the binary evolves detached until the primary reaches core hydrogen depletion at 7.6 Myr. At this stage, the primary has lost ≃2.97 M_⊙ due to the stellar winds. The primary evolves off the main sequence and expands past its Roche lobe radius, thus initiating the first mass transfer phase.

2. During the mass transfer phase of ∼16 kyr, the companion is spun up rapidly, while accreting only ≃0.27 M_⊙. Any additional transferred material (≃9.49 M_⊙) is quickly lost by the rotation-enhanced winds. This initial mass transfer phase removes most of the hydrogen envelope from the primary star before detaching, leaving a hydrogen-rich envelope of ≃2.43 M_⊙ and a surface hydrogen fraction of ≃0.33.

3. The low surface hydrogen fraction signals that Wolf-Rayet winds have kicked in and the remaining hydrogen envelope is quickly removed in the subsequent detached phase until a fully stripped helium star remains.

4. At ∼8.0 Myr the secondary, now 21.4 M_⊙, expands due to core hydrogen exhaustion and initiates a second mass transfer phase. It loses ≃9.69 M_⊙, while the fully stripped primary star only accretes ≃0.04 M_⊙ before it reaches critical rotation. When the secondary reaches core helium ignition, it exits the mass transfer phase at a stage similar to the primary with a limited hydrogen envelope of ≃2.22 M_⊙ and a similar surface hydrogen fraction of ≃0.32. The properties between the primary and secondary are similar because they start with a mass ratio of 0.95 but the primary is now rapidly rotating.

5. The primary continues to evolve for another ∼100 kyr and remains compact until core carbon exhaustion. It loses only 0.06 M_⊙ during this phase and the primary maintains a sufficiently large angular momentum reservoir in the internal structure to form a 1.46 M_⊙ accretion disk during the collapse into the BH, as shown in Figure 2. Stellar mass shells with angular momentum above j_isco form an accretion disk around the BH, as described in Section 2.2. This accretion disk results in a total energy budget of E_BZ/η_ϕ ≈ 2.60 × 10⁵⁴ erg. With an η_ϕ of 10⁻³, this system has a typical LGRB energy of ∼10⁵¹ erg and, with the binding energy of the primary only being ∼10⁵⁰ erg, the jet is likely to breakout successfully.

   Fig. 2. Specific angular momentum distribution of the mass shells of the primary star in the example binary shown in Figure 1, at core carbon depletion in red. The kink at ≃7.8 M⊙ is the transition between the carbon-oxygen core and helium shell. It forms an accretion disk of 1.46 M⊙ during the collapse. For comparison, we plot the specific angular momentum of the ISCO, jISCO, for a non-rotating Schwarzshild BH (dashed black) and a maximally rotating Kerr BH (dotted black), as a function of BH mass. Additionally, we track the evolution of jISCO for the example primary during the collapse as its mass and spin change (dashed green). It tracks jISCO at the collapse of each mass shell and starts at 3 M⊙ with a BH of 2.5 M⊙ and 0.5 M⊙ neutrino mass loss.

The essential components of this LGRB formation channel are the occurrence of stable reverse mass transfer and the angular momentum reservoir at collapse. Since these have complex dependencies on the initial properties of the binary and vary across metallicity, in what follows we first discuss the conditions for the secondary to fill its Roche lobe before the primary star has collapsed at solar metallicity. Then we discuss for what initial condition this interaction is stable and leads to disk formation at collapse. Finally, we investigate why this LGRB formation channel only occurs at near-solar metallicities.

3.1.2. Dependence on mass ratio, period, and donor mass

In our example model, both stars fill their Roche lobe on the Hertzsprung gap after completing hydrogen burning in their cores. This moment of interaction is common in the proposed LGRB formation channel and requires the main-sequence lifetimes of the binary components to be similar. Otherwise, the primary star will have reached carbon exhaustion before the secondary has left the main sequence. This evolutionary timescale difference is influenced by the mass of the primary and the binary mass ratio. Specifically, the more equal the masses of the stars in the binary, the smaller the lifetime difference. As a result, many systems with q close to unity have the secondary fill its Roche lobe before the primary completes its evolution. This effect can be seen in Figure 3, where example grid slices of solar metallicity POSYDON detailed binary models with initial mass ratios of q = 0.85, q = 0.90, q = 0.95, and q = 0.99 are shown. The green markers indicate systems with the secondary expanding past its Roche lobe, while the orange markers indicate systems with only the primary star filling its Roche lobe. Systems with both stars filling their Roche lobe at the same time are marked with purple, as contact systems. As the mass ratio approaches 1, the green area becomes larger, indicating that more systems have reverse mass transfer. Additionally, the minimum mass ratio for the secondary to fill its Roche lobe depends on the absolute mass of the donor star. Massive stars are dominated by radiation pressure, which causes their main-sequence lifetimes to become similar towards higher masses. This means that binaries with more massive primaries have a broader range of mass ratios where the secondary can initiate a reverse mass transfer phase. As a result, the highest mass primaries in Figure 3, still have reverse mass transfer at q = 0.85, while the lower mass systems do not.

Fig. 3. Four grid slices at Z⊙ with q = 0.85, 0.90, 0.95, and 0.99. Orange symbols represent systems where only the primary fills its Roche lobe throughout the evolution of the binary model. Green squares indicate systems where also the secondary fills its Roche lobe, but not at the same time as the primary. Purple symbols indicate systems where both stars fill their Roche lobe at the same time. The type of symbol shows the stability of the interaction: squares indicate stable mass transfer, while binary models reaching instability criteria are indicated with a diamond. Gray squares indicate no mass transfer, while black dots indicate a non-converged model or Roche lobe overflow when the model is initiated. For systems where the companion fills its Roche lobe and the interaction remains stable, we indicate if an accretion disk is formed during the collapse of the model. This is indicated by a coloured circle inside the square symbol, where the colour indicates the mass of the accretion disk.

In Figure 3, it is important to address that the secondary does fill its Roche lobe in contact systems, but often reaches the L₂ overflow criteria for unstable mass transfer or, at longer periods, one of the stars leaves the main sequence while in contact. In general, only when the initial period is long enough can the secondary fill its Roche lobe without the binary reaching an unstable contact state.

A stable mass transfer phase is a requirement to provide sufficient angular momentum to the primary star to produce an accretion disk at collapse. In Figure 3 at q = 0.95, two islands of unstable reverse mass transfer are present in the regime for Roche lobe overflow initiated by the secondary. The higher period island around 10³ days is a result of reaching the maximum mass transfer rate, while the lower period island near 10² days is a result of L₂ overflow. While reverse mass transfer occurs more often nearing q = 1, the grid at q = 0.99 in Figure 3 at solar metallicity shows that this is often unstable. Many binaries at q = 0.99 reach a contact phase due to both stars evolving on a similar timescale, which causes many massive donors to reach contact off the main sequence even for large separations. The island of L₂ overflow remains and the diagonal unstable mass transfer is caused by the maximum mass transfer rate. There is still a large range of stable reverse mass transfer present at q = 0.99, but not as extensive as at q = 0.95 and many of the primary masses are too small to form a BH at collapse.

The final limitation in the formation of LGRB through the stable reverse-mass-transfer channel is the angular momentum content of the primary at collapse. The circles inside the squares in grid slices in Figure 3 show the stable reverse-mass-transfer binaries that form an accretion disk during their collapse. Their colour indicates the mass of the accretion disk that is formed. The q = 0.90 and q = 0.95 slices show a gradient in the disk mass with donor mass and period. In general, wider systems produce a less massive accretion disk than shorter-period binaries, and the widest periods with stable reverse mass transfer do not produce a rapidly rotating star at collapse. The boundary of accretion disk formation moves to shorter periods as the donor mass increases and is a result of hydrogen being present when the reverse mass transfer occurs, which inhibits the transfer of angular momentum into the core of the primary.

We explore the gradient in the accretion disk mass as a function of donor mass and period in Appendix A. In short, at the same mass and mass ratio, the accreted angular momentum and mass during the stable reverse-mass-transfer phase decrease with increasing period due to the mass transfer being shorter and, thus, there being less time to diffuse the angular momentum into deeper layers of the accretor. After the stable reverse-mass-transfer phase, the fraction of angular momentum loss is similar across different initial periods, around 80%. Thus, short-period systems, which gained more angular momentum through mass transfer, reach core carbon depletion with a larger angular momentum reservoir, sufficient to produce an accretion disk during collapse.

The gradient as a function of donor mass originates from a slightly different physical phenomenon. At the same period and mass ratio with an increasing initial donor mass, more angular momentum is gained by the primary during the reverse-mass-transfer phase. However, the initially more massive primary loses more angular momentum and mass due to stellar winds after the stable reverse mass transfer. A lower total angular momentum amount and the increased total mass at core carbon depletion limit the mass of the formed accretion disk. This effect is visible in Figure 3, where at the same period, an initially more massive star in the regime of stable reverse mass transfer produces a smaller or no accretion disk compared to a lower mass primary.

3.1.3. Metallicity dependence

The criteria for stable reverse mass transfer with sufficient angular momentum at collapse provide a limited parameter space at solar metallicity to potentially produce LGRBs, resulting in a rate of ≃10⁻⁵ events yr⁻¹ Gpc⁻³. The efficiency of this formation channel decreases with decreasing metallicity, with the rate going to zero below 0.2 Z_⊙. We attribute this to a combination of effects based on metallicity, which we highlight in Figure 4, where we show grid slices of detailed binary models with an initial binary mass ratio of q = 0.95, for a range of metallicities.

Fig. 4. Grids slices at q = 0.95 at Z = 10−3 Z⊙, 0.2 Z⊙, and 0.45 Z⊙. The same symbol and colours are used as in Figure 3. Less stable reverse mass transfer (green squares) occurs as the metallicity decreases and those that do, do not produce accretion disks upon their collapse.

At Z = 0.45 Z_⊙, the region with sufficient angular momentum at collapse in the stable reverse-mass-transfer channel shrinks, as shown in Figure 4. Weaker main-sequence winds and partial stripping of the primary star during the first mass transfer phase (Götberg et al. 2017) lead to the formation of a more massive helium core. Additionally, if a hydrogen envelope is still present during the reverse mass transfer, the angular momentum transferred to the primary does not propagate into the core due to the abundance gradient between the core and the hydrogen shell. This results in less angular momentum gain during the reverse mass transfer, and as the primary loses its hydrogen envelope before collapse, it loses the majority of its angular momentum, preventing the formation of an accretion disk. At the same time, the islands of unstable mass transfer from L₂ overflow and the maximum Ṁ are expanding.

Moving to Z = 0.2 Z_⊙, two additional effects are starting to show: short-period systems no longer experience reverse mass transfer, while long-period systems are now experiencing a contact phase, which becomes dominant towards even lower metallicities, as the Z = 10⁻³ Z_⊙ panel in Figure 4 shows. In wider binaries, the mass transfer no longer sufficiently strips the primary star for it to reach the Wolf-Rayet wind regime, instead, it finds a thermal equilibrium close to or at its Roche lobe (Klencki et al. 2022; Ercolino et al. 2024). As the companion evolves off its main sequence and also fills its Roche lobe, the system goes into contact either because the primary is already filling its Roche lobe or because its equilibrium state is near the Roche lobe so that the radial expansion from mass accretion causes the primary to fill its Roche lobe. On the other hand, short-period binaries no longer experience reverse mass transfer, because the companion remains compact and hydrogen-burning until the primary reaches carbon exhaustion. Due to lower stellar wind mass loss, these accretors retain the gained angular momentum from the initial mass transfer phase. At lower metallicities, the radial expansion of the stars is already less compared to solar metallicity due to the lower opacity, and the increased rotational mixing further limits the radial expansion, which restricts a reverse-mass-transfer phase. The combination of these effects creates a LGRB formation channel that becomes more efficient as metallicity increases, while completely disappearing at metallicities below 0.2 Z_⊙.

3.2. Stable reverse-mass transfer LGRB rate

The combined requirements of stable reverse mass transfer and the angular momentum reservoir at collapse provide strong constraints on the formation of LGRB progenitors: the reverse-mass-transfer channel only contributes at Z ≥ 0.2 Z_⊙. As described in Section 2, we use POSYDON to generate burst stellar populations at eight different metallicities and select systems undergoing stable reverse mass transfer that will form an accretion disk on collapse. We convolve these populations with the detailed star formation rate and metallicity distribution across cosmic time from the IllustrisTNG simulation in order to compute the cosmic LGRB rate from the stable reverse-mass-transfer channel. The black line in Figure 5 shows the resulting total cosmic rate of systems that produce an accretion disk on collapse. This is the maximum rate that can be achieved by reverse-mass-transfer progenitors according to our model, in the presence of a highly efficient jet-launching and powering mechanism, where every system with an accretion disk successfully produces a LGRB. This would correspond to a rate of ≃100 events per cubic gigaparsec per year at z = 0. For comparison, the most recent LGRB rate density estimate based on detailed modeling of the observed population, carefully accounting for selection effects, put the total intrinsic cosmic rate at z = 0 at ${79}_{-33}^{+57}$ Gpc⁻³ yr⁻¹ (Ghirlanda & Salvaterra 2022).

Fig. 5. Redshift evolution of the stable reverse-mass-transfer LGRB event rate density for the IllustrisTNG metallicity and star formation history. The total potential population is indicated in black, while the dashed red line shows the intrinsic LGRB rate from Ghirlanda & Salvaterra (2022) with the best-fit parameters based on observations. The red shaded area is 10% to 20% of this rate. The five green shaded lines indicate the rate density for an energy limit of 1051 erg and different ηϕ.

Given that most observed LGRBs occur in low-metallicity environments, this indicates that some effect must be limiting a large fraction of our stable reverse mass-transfer model progenitors from producing a LGRB. Theoretically, one would expect that such an effect follows from the required energy for the jet to propagate through the stellar envelope and successfully break out of it (e.g., Matzner 2003; Bromberg et al. 2011; Salafia et al. 2020). The binding energy of the stellar structure at core carbon exhaustion is less than the values of E_BZ/η_ϕ, when η_ϕ = 1. While the actual energy breakout in general depends on the properties of the stellar envelope and also on those of the jet at launch. Here, for simplicity, we model both as a universal threshold on the jet energy. Because we left the efficiency factor η_ϕ unspecified (see Section 2), we set an energy limit for a successful LGRB of 10⁵¹ erg based on the typical breakout energy from jet propagation simulation (for example, see Harrison et al. 2018). In what follows, we explore five values for the efficiency, η_ϕ, namely 1, 0.1, 10⁻², 10⁻³, and 5 × 10⁻⁴.

The total rate of disk formation through the stable reverse-mass-transfer channel is comparable to the local rate density of observed LGRBs (Ghirlanda & Salvaterra 2022). Given that the inferred fraction of high-metallicity host galaxies is 10%–20%, this requires an efficiency that removes 80%–90% of events. The most restrictive efficiency with η_ϕ = 5 × 10⁻⁴ results in a volumetric rate of ≃2.01 Gpc⁻³ yr⁻¹ at z = 0; with a total of ≃15.0 Gpc⁻³ yr⁻¹ events below z < 2.5. With η_ϕ = 10⁻³, the reverse-mass-transfer channel rate is around 18% of the latest estimate of the observed cosmic LGRB rate. This is remarkably close to the observational fraction of solar and super-solar host galaxies of 10%–20%. However, because the normalization of the observed rate is uncertain due to selection effects, such as the intrinsic energy of the jet, one should be careful to interpret a comparison of our theoretical rates against these observations. These cannot be constrained without considering all formation channels concurrently, which is beyond the scope of this work. The intrinsic rate density evolution with redshift for different choices in the efficiency is shown in Figure 5. We report their rate densities at z = 0 and averaged within z = 2.5 in Table 1.

In Figure 6 we show the cumulative distribution of E_BZ with η_ϕ = 10⁻³ for all our simulated progenitors within z < 2.5 (black dashed line) and for progenitors with specific metallicities within the same redshift range (coloured lines). We limit the redshift range to z < 2.5 to cover most observable LGRBs. Including higher redshift systems reduces the average volumetric rate since fewer events from the stable reverse-mass-transfer channel occur at low metallicity.

Fig. 6. Cumulative rate density (R[E > Ex]) within z < 2.5 from the IllustrisTNG simulation split up per metallicity with ηϕ = 10−3. The vertical line indicates a typical breakout energy of 1051 erg.

Although the progenitor masses at collapse only cover a range from 10 to 40 M_⊙, E_BZ covers several orders of magnitude from 10⁴⁴ to 10⁵² erg. The most energetic events come from the solar metallicity population, although their contribution to the rate changes depending on the star formation history (SFH; see Section 4). The median value of E_BZ for these models is around ≃10⁴⁹ erg, but a tail of much higher-energy jets is present, with clustering around E_BZ ∼ 10⁵¹ erg, which causes the steep decrease in the cumulative rate from ≃100 Gpc⁻³ yr⁻¹ to ≃10 Gpc⁻³ yr⁻¹ around this energy.

4. Discussion

While LGRBs have been found in high-metallicity environments, no collapsar formation channel with a positive metallicity relationship has previously been proposed. We have introduced a new formation channel capable of forming sufficient LGRBs to match observational or theoretical fractional rates in high-metallicity environments. However, much uncertainty remains in this fraction, as well as in the absolute rate of the total LGRBs and this channel. Since the main contributors to this uncertainty are the observational biases and hidden LGRB fraction, which would require a population study that includes all possible LGRB formation pathways in order to be explored, we instead focus on the physics which might impact the presence and occurrence of the reverse-mass-transfer LGRB formation channel.

4.1. Implications of energy limit

We have implemented a universal energy limit as a proxy for successful jet breakout with different efficiencies, which is equivalent to different energy limits with a constant efficiency, i.e. a 10⁵¹ erg limit with a η_ϕ = 10⁻³ selects the same progenitors as a 10⁵⁴ erg limit with a η_ϕ = 1. For the stable reverse-mass-transfer channel to produce 10–20% of the local LGRB rate, a low disk-to-jet-energy conversion or a high jet-breakout-energy threshold is required. This implies that successful jets are a minority, while choked jets occur an order of magnitude more frequently in the proposed channel. This fraction aligns with Bromberg et al. (2012) and from observed low-luminosity GRBs the inferred intrinsic rate can be up to an order of magnitude higher than their high-luminosity counterpart (Liang et al. 2007; Guetta & Valle 2007; Virgili et al. 2009). It suggests the existence of a large population of relativistic supernovae with potentially low-luminosity gamma-ray emission caused by these choked jets (see, for example, Soderberg et al. 2006; Wang et al. 2007; Nakar & Sari 2012).

Despite the broad range of jet energies, a high jet-breakout-energy limit or low disk-to-jet-energy conversion leads to a narrow energy range for successful jets, which declines steeply above the minimum energy. Moreover, successful jets possess only slightly more energy than needed for breakout implying the presence of a strong cocoon with comparable energy to the jet. For a structured jet, this results in a shallow angular profile (Salafia et al. 2020). Thus, the apparent diversity of LGRBs from this channel would primarily be due to viewing angles rather than an actual broad intrinsic energy distribution.

4.2. Type Ic broad-line and potential LGRB companions

Long gamma-ray bursts and SN Ic-BL have been associated with one another, but this relationship is not universal; not all Type Ic-BL are accompanied by a LGRB and vice versa (Modjaz et al. 2016; Japelj et al. 2018). Type Ic-BL without a LGRBs have distinct spectral features from those with a LGRB (Modjaz et al. 2016), and could indicate a failed jet breakout (see Corsi & Lazzati 2021, and references therein). In general, Type Ic-BL are associated more with super-solar host galaxies compared to LGRBs (Japelj et al. 2018), which aligns with the positive metallicity dependence of the stable reverse-mass-transfer LGRB channel. Additionally, if this channel produces the observed high-metallicity population of LGRBs, it produces an order of magnitude more choked jets compared to successful jets (see Section 4.1). Although the observed relative rate of SN Ic-BL to core collapse SNe is poorly constrained, based on a single SN Ic-BL observation in the volume-limited Lick Observatory Supernova Search, SN Ic-BL make up 4% of stripped-envelope supernovae (Shivvers et al. 2017). Combined with the stripped-envelope supernova rate from Frohmaier et al. (2021), this results in a rough SN Ic-BL rate estimate of ∼300 Gpc⁻³ yr⁻¹. Assuming that all disk-forming systems produce a jet and η_ϕ = 10⁻³, ≃86 choked jets occur per Gpc³ per year at z = 0 from the stable reverse-mass-transfer channel. This is comparable to the observationally inferred rate and could imply that either more dominant LGRB formation channels must have a smaller ratio of choked-to-successful jets to not overproduce the SN Ic-BL or not all failed LGRB accretion disk systems produce a SN Ic-BL.

More puzzling are the three local LGRB events without an associated Type Ic-BL (e.g., Fynbo et al. 2006; Gal-Yam et al. 2006; Xu et al. 2009). GRB 111005A, one of the three LGRB in a high-metallicity environment, is one such LGRB (Tanga et al. 2018). Their origin is poorly understood in the framework of the jet moving outwards through the stellar structure, producing a Type Ic-BL supernova, and could indicate a new LGRB formation mechanism (e.g., Gal-Yam et al. 2006). Either way, the near-solar metallicity environment does not have to be an indicator for a unique mechanism, as the stable reverse-mass-transfer channel can produce LGRBs at those metallicities through the standard collapsar scenario.

A special feature of the stable reverse-mass-transfer channel is the presence of a (partially) stripped companion at collapse. Its mass ranges from 10 M_⊙ to 50 M_⊙ with most potential LGRB progenitors having a companion around 15 M_⊙ at collapse. While no clear LGRB in a binary system has been observed (although see Ivy Wang et al. 2022), the presence of a companion will affect the GRB emission (Zou et al. 2021). Furthermore, several local Type Ic-BL supernovae have been investigated for companions with limited success. For SN 2002ap, an upper limit ∼10 M_⊙ for a companion has been set (Zapartas et al. 2017), while a potential companion in SN 2016coj is under debate (Terreran et al. 2019). The stripped nature of the companion could complicate its observability. However, a 15 M_⊙ helium star would appear as a Wolf-Rayet star at Z_⊙ (Shenar et al. 2020). The detection of or upper limits on a companion in nearby LGRB and SN Ic-BL can provide an independent observable on their formation mechanisms.

4.3. Jet physics

The relationship between the progenitor star and jet energy remains highly uncertain and general relativistic, radiation, magnetohydrodynamic simulations have shown that the efficiencies depend on the state of the disk (Tchekhovskoy et al. 2010). A popular approach is the magnetically arrested disk (MAD) state because it produces powerful jets even for slowly rotating BHs and causes rapidly spinning BHs to spin down (Jacquemin-Ide et al. 2024). In this state, the mass accretion rate and spin-down timescale determine the jet power. In such disks η_ϕ = 1, while η_a can be fitted following Lowell et al. (2024). This can result in efficiencies above 1 due to the extraction of angular momentum from the BH and results in high-energy jets. These even occur for slowly spinning BHs as long as the progenitor has sufficient angular momentum to form an accretion disk (Fryer et al. 2025). However, recent simulations of collapsars show that the accretion rate remains sufficiently high that neutrino cooling keeps the disk thin and the MAD state is not reached when assuming a proto-neutron star is the origin of the BH magnetic field (Gottlieb et al. 2024).

In a non-MAD state, the magnetic field is thought to determine the jet power, which Kawanaka et al. (2013) estimated using the ram pressure at the inner-most stable orbit around the BH following Beckwith et al. (2009). This provides a $L_{\mathrm{BZ}}\propto \dot{M}$ relationship, similar to that implemented here and as found by Blandford & Znajek (1977) and Heger et al. (2003), although the latter uses a 10^{0.1/(1 − a_BH)−0.1} spin dependence. Although dependent on the progenitor structure, many of the energetic LGRBs have a high accretion rate through their collapse, making the chosen relationship for L_BZ reasonable. On the other hand, E ∝ a² might require improvement (see, for example: Tchekhovskoy et al. 2010; Russell et al. 2013). Furthermore, we do not model the extraction of angular momentum from the BH while it collapses, which will reduce the total extracted energy. On the other hand, we partially account for this with different energy limits. As a consensus arises on the jet power mechanism in the future, this relationship can be improved.

Besides the potential energy budget to form a jet, an essential component in LGRB formation in the collapsar scenario is the jet breakout through the stellar structure. The jet breakout depends on the injected energy and active duration of the central engine, and the density profile of the stellar structure. While the internal stellar density profile at core-carbon exhaustion is available, we do not evolve the final evolutionary phases due to the required computational time, and, thus, do not have a representative profile at collapse available. We verify that the free-fall timescales of the core-carbon exhausted models are similar to observed T₉₀ values (See Appendix C). Furthermore, the injected energy depends on several poorly constrained parameters, such as the injection opening angle, the disk mass to jet energy conversion, and the total jet energy emitted in the gamma-ray. A full end-to-end pipeline with computationally intensive hydrodynamical simulations is required to accurately model the link between the progenitor model and the jet properties from the collapsar scenario (Gottlieb et al. 2024).

In general, failed jet breakouts are more likely to occur in the low-potential energy models, as these models provide insufficient mass to sustain a successful jet propagation through the stellar envelope. In contrast, a high-potential energy model is more likely to overcome the energy requirement for jet breakout and cause an observable LGRB. The proposed limits show that only the most energetic systems from the stable reverse-mass-transfer channel must be able to lead to LGRBs in order to match the observed fractions of high-metallicity ones at z < 2, which can be achieved with an η_ϕ ≲ 10⁻³.

4.4. Star formation history

Although the total star formation rate density in the local Universe is relatively well constrained, the metallicity evolution and spread are less understood. It is the high-metallicity-specific star formation rate that gives rise to uncertainty in the shape of the stable reverse-mass-transfer LGRBs and in the weighting of each model over redshift. The former can be seen in Figure 7, where we compare the LGRB rates over redshift from the IllustrisTNG simulation with the two empirical star formation history prescriptions described in Section 2. The metallicity distribution redistributes the weight between populations of different metallicities but does not alter the shape of the potential LGRB energy distribution per metallicity (see Appendix B). The ln-normal metallicity distribution from Neijssel et al. (2019) produces a very peaked metallicity-specific star formation rate, reducing the contribution from the reverse-mass-transfer LGRB channel at high redshift to the LGRB rate, while the higher metallicity and wider distribution from Bavera et al. (2020) allow for more near-solar metallicity star formation at higher redshifts.

Fig. 7. Total potential stable reverse-mass-transfer LGRB rate density over cosmic time from the IllustrisTNG simulation, Madau & Fragos (2017) with a log-normal metallicity distribution with σ = 0.5, and Neijssel et al. (2019) with a ln-normal distribution with σ = 0.39. No energy requirements are implemented.

The IllustrisTNG simulation, on the other hand, is rapidly enriched compared to the Neijssel et al. (2019) prescription, which boosts the LGRBs from the stable reverse-mass-transfer channel at high redshift and is similar to the metallicity distribution from Madau & Fragos (2017). Due to its metallicity distribution, the Madau & Fragos (2017) rate extends high-metallicity star formation to higher redshifts, resulting in even higher rates than the fast-enriching TNG simulation. Despite these differences at high redshift, the differences in the local Universe (z < 2.5) are of a factor of 2.

A stronger effect in the local Universe comes from the energy cuts. Due to the reweighing of each LGRB progenitor, the cumulative potential energy distribution is significantly different per star formation history. Thus, the energy cut can suddenly remove more events. For example, the most restrictive cut in Figure 6 reduces the stable reverse-mass-transfer LGRB rate to 2.01 LGRB Gpc⁻³ yr⁻¹, while the Neijssel et al. (2019) prescription in Figure B.4 reduces it to 0.17 LGRB Gpc⁻³ yr⁻¹. This is due to the energy distribution at each metallicity, which each SFH rebalances, as can be seen in Appendix B. As a result, the choice in mean metallicity and metallicity spread has a major impact on the final rate and the limit of the jet energy.

4.5. Additional caveats

An essential step in the metallicity-dependence of the stable nature of the reverse mass transfer is the further stripping of the partially stripped primary star by stellar winds. In most cases, before the first mass transfer, the donor star lies in the metallicity-dependent mass loss regime for OB stars from Vink et al. (2001). At high metallicities, this already removes significant material from the primary star, such that it reaches the strong Wolf-Rayet winds after the first mass transfer phase ends. These winds quickly remove any remaining hydrogen envelope. Their strength determines if the star becomes fully stripped and the amount of angular momentum remaining after a reverse mass transfer event. The Wolf-Rayet and the main-sequence massive star mass loss decreases with decreasing metallicity, reducing the possibility for the primary to become fully stripped. At the same time, the mass transfer is less efficient at removing the complete hydrogen envelope at lower metallicities than solar. This is due to the primary star finding an equilibrium state close to its Roche lobe with a significant hydrogen envelope left Klencki et al. (2022), which is an internal property of the star and not dependent on the stellar wind prescriptions.

If the metallicity dependence for O stars follows a shallower relation of Ṁ ∼ Z^0.42 as suggested by Vink & Sander (2021), the stable reverse-mass-transfer scenario is extended to lower metallicities. On the other hand, empirically derived mass loss rates of the Large and Small Magellanic Clouds have indicated stronger metallicity-dependencies for massive stars of Z² and Z^0.83, respectively (Mokiem et al. 2007; Ramachandran et al. 2019). Either way, the number of stable reverse-mass-transfer systems will increase or decrease depending on the metallicity-dependence of the massive star winds, which will impact the preferred potential energy limit to match the observed high-metallicity LGRB fraction. Irrespective of the metallicity-dependence of the mass-loss rate, the mass-loss rates for OB stars at solar metallicity remain similar (Vink & Sander 2021), thus minimally affecting the contribution of the near-solar metallicity populations.

After the stable reverse-mass-transfer phase, the surface hydrogen fraction of the primary can return to above 0.4, which is the boundary for the Wolf-Rayet winds. As such, the stellar winds for OB stars are applied to this phase of the evolution, which will transition back into the Wolf-Rayet winds as the small hydrogen layer is stripped and the surface hydrogen fraction reaches below 0.4. The transition to the Vink et al. (2001) winds means a lower mass loss in comparison to the Nugis & Lamers (2000) Wolf-Rayet winds are applied. If the mass loss during this phase were treated as a Wolf-Rayet wind, more mass and angular momentum would be lost in the final phase, which would shrink the number of primaries with sufficient angular momentum to produce an accretion disk. On the other hand, the total mass and angular momentum loss strongly depend on how far the core helium burning has progressed when the stable reverse mass transfer takes place. Generally, for potential LGRB progenitors in the stable reverse-mass-transfer channel, the interaction occurs late in the core helium-burning phase of the primary. Otherwise, the hydrogen envelope is not yet fully stripped, which inhibits the transport of angular momentum into the helium core. With no clear observations of hydrogen-accreted helium stars, the exact treatment of their mass loss is unclear.

The first mass transfer in the binary can result in rejuvenation of the companion, where the efficiency of semi-convection affects the parameter range for which this occurs (e.g., Braun & Langer 1995). The POSYDON v2 grids use α_sc = 0.1 (for more details, see Andrews et al. (2025)), a relatively inefficient semi-convection, which limits the parameter space for rejuvenation. Furthermore, the stable reverse-mass-transfer channel predominantly occurs around q ≈ 0.95 at Z_⊙. This implies that the companion is far along the main-sequence when the first mass transfer takes place. This limits the effect of the accreted material on the further evolution and the timing of the reverse mass transfer phase. It is unlikely that more efficient semi-convection can lead to more rejuvenation and less stable reverse mass transfer in this regime. At low metallicities, the first mass transfer does affect the companions, and they remain compact until the primary reaches carbon-exhaustion. More efficient semi-convection could result in more rejuvenation, while less efficient semi-convection could limit the rejuvenation and allow for a reverse mass transfer before the primary collapses. However, these interactions are likely to be unstable, since the primary will not be fully stripped at the moment of the interaction. As such, the choice of semi-convection parameter should not affect the overall results of this work.

Although we self-consistently calculate the collapse of the progenitor into the BH, we first determine the formation of a BH using the Patton & Sukhbold (2020) prescription. Individual models in this explodability landscape are sensitive to the numerical setup, such as the nuclear reaction network and spatial resolution (Farmer et al. 2016). These variations could be indicative of inherent chaos in the final stages of massive star evolution (Sukhbold et al. 2018), but introduce uncertainty into the mapping between the core-carbon depleted structure and compact object formation. Additionally, this prescription assumes that implosions lead to BH formation and explosions lead to neutron star formation.

While the majority of LGRB progenitors from the stable reverse-mass-transfer channel are well into the BH formation regime, their parameter space at low masses is limited by neutron star formation, as shown in the right grid slice in Figure 3. Allowing for more BH formation would increase the contribution of this formation channel. However, the additional parameter space is limited, because the timescale difference between the binary components is too large at low masses for reverse mass transfer to occur. Additionally, the Patton & Sukhbold (2020) prescription uses the average central carbon abundance and carbon-oxygen core mass at helium depletion to determine the explodability. Late-stage binary interactions, such as mass transfer after core helium ignition, can alter the stellar structure sufficiently to boost explodability (Laplace et al. 2021). Despite that, the POSYDON MESA models are evolved until carbon depletion, due to the Patton & Sukhbold (2020) prescription, late-stage mass transfer effects are not considered.

Switching the remnant mass prescription to the rapid or delayed prescriptions from Fryer et al. (2012) should give an indicator of how different prescriptions for estimating the remnant masses affect our results. These prescriptions use the core properties at carbon depletion to determine BH formation through fallback onto a proto-compact object. Figure 8 shows minimal differences between their metallicity bias functions, which propagate to nearly the same cosmic rates. Slightly less stable reverse-mass-transfer LGRBs occur with the Fryer et al. (2012) prescriptions. The star formation history, metallicity distributions, and energy cut have a significantly stronger impact on this rate than what stars collapse as BHs.

Fig. 8. Number of stable reverse-mass-transfer LGRBs per solar mass star formation (efficiency) as a function of metallicity. These metallicity bias functions are shown for the two Fryer et al. (2012) (rapid in purple; delayed in green) and the Patton & Sukhbold (2020) (orange) prescriptions.

5. Conclusions

Although LGRBs are more frequently observed in low metallicity environments, 10 to 20% of LGRB host galaxies at z < 2 exhibit solar and super-solar metallicities (Vergani et al. 2015; Palmerio et al. 2019; Graham et al. 2023). Moreover, several LGRB events have been associated with high-metallicity explosion sites (Levesque et al. 2010a; MichałowskI et al. 2018; Heintz et al. 2018). Traditional binary LGRB formation pathways predominantly have a low-metallicity preference and are rare to non-existent at solar metallicity. In this work, we have proposed a new collapsar LGRB formation pathway through stable reverse mass transfer that is limited to near solar-metallicity stellar binaries, creating a unique formation mechanism which has the opposite dependence on metallicity compared to traditional scenarios.

In this formation scenario, a partially stripped star is formed by mass transfer after hydrogen depletion (case B). At solar metallicity, it becomes completely stripped by stellar winds from before (OB-type stellar winds) and after the mass transfer (Wolf-Rayet star winds). These wind regimes introduce the metallicity dependence to this LGRB formation channel. The stripped star becomes rapidly rotating when case B mass transfer is initiated by the initially less massive star in the binary system. Mass accretion onto the stripped star carries sufficient angular momentum, which the star is able to maintain till the end of its evolution and form a substantial accretion disk around the proto-BH during collapse.

The existence of this stable reverse-mass-transfer pathway for forming rapidly rotating stripped stars at core collapse is a robust prediction of population synthesis with POSYDON. However, linking these progenitors to jet breakout and observable LGRBs on a population scale is significantly affected by theoretical uncertainties and observational biases. As such, the absolute rate is highly dependent on the star formation history, metallicity evolution and successful LGRB energy limits. Nevertheless, sufficient angular momentum, and thus energy, is available in the structure to power jets with LGRB energies around 10⁵² erg. Moreover, the rate is sufficiently high to produce 10 to 20% of the intrinsic observed LGRB rate from Ghirlanda & Salvaterra (2022). This contribution is consistent with the fraction of LGRBs observed in high-metallicity environments.

However, the absolute rate of LGRBs from the stable reverse-mass-transfer channel can vary between 1 and 100 LGRB Gpc⁻³ yr⁻¹ depending on the star formation rate, metallicity evolution, and energy threshold for a successful LGRB. To refine the contribution of this high-metallicity LGRB formation pathway, a comprehensive, simultaneous population analysis of all LGRB formation channels is required, which will benefit from a more extensive exploration of the link between rapidly rotating stellar models and successful LGRB observations.

Acknowledgments

We thank the anonymous reviewer for their insightful comments on the manuscript. The POSYDON project is supported primarily by two sources: the Swiss National Science Foundation (PI Fragos, project numbers PP00P2_211006 and CRSII5_213497) and the Gordon and Betty Moore Foundation (PI Kalogera, grant award GBMF8477). MMB was supported by the Boninchi Foundation, the project number CRSII5_21349, and the Swiss Government Excellence Scholarship. OSS acknowledges funding by the INAF Finanziamento per la Ricerca Fondamentale grant no. 1.05.23.04.04. GG and OSS acknowledge financial support by the INAF Grant “GRB PrOmpt Emis- sion Modular Simulator” (POEMS – 1.05.23.06.04) and European Union-Next Generation EU, PRIN 2022 RFF M4C21.1 (PEACE – 202298J7KT). JJA acknowledges support for Program number (JWST-AR-04369.001-A) provided through a grant from the STScI under NASA contract NAS5-03127. PMS, SG, KAR, and MS were supported by the project numbers GBMF8477 and GBMF12341. K.A.R. is also supported by the Riedel Family Fellowship and thanks the LSSTC Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining Grant No. 1829740, the Brinson Foundation, and the Moore Foundation; their participation in the program has benefited this work. KK is supported by a fellowship program at the Institute of Space Sciences (ICE-CSIC) funded by the program Unidad de Excelencia María de Maeztu CEX2020-001058-M. Z.X. acknowledges support from the China Scholarship Council (CSC). EZ acknowledges support from the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “3rd Call for H.F.R.I. Research Projects to support Post-Doctoral Researchers” (Project No: 7933).

References

Appendix A: Accretion disk dependence on initial primary mass and period

Within the stable reverse-mass-transfer region, the mass of the accretion disk at collapse depends on the initial primary mass and period. As the initial period or donor mass increase, the accretion disk mass decreases. As such, the lowest mass and shortest period systems within the stable reverse-mass-transfer region, in general, have the most massive accretion disk when collapsing. In this section, we will first cover the initial period dependence and then the initial primary mass dependence. We pick representative slices of the stable reverse-mass-transfer region to show the evolution of different stable reverse-mass-transfer binaries.

A.1. Initial period dependence

At Z_⊙, we select systems with M_ZAMS, 1 ≈ 39.87 and q = 0.95, which are undergoing stable reverse mass transfer. These systems are illustrated in Figure 3 and exhibit a correlation between decreasing disk mass and increasing initial period. Only the widest two systems that form an accretion disk at collapse are reversed in disk mass.

In Figure A.1, we show the evolution of these models focusing on the mass transfer phases. As expected, the shortest period systems (dark blue) fill their Roche lobe already on the main-sequence, as shown in the third row. They shrink back inside their Roche lobe until core hydrogen depletion is reached. At this stage, all primaries expand, and for the wider systems (green and yellow), this is their first mass transfer phase. The hydrogen envelopes of the primaries in long-period binaries experience less stripping compared to those in short-period systems; consequently, they exhibit a greater radial expansion after the mass transfer. The P_ZAMS = 37.28 days binary exits the first mass transfer with 4.39 M_⊙ of hydrogen envelope left and a radius of ∼84R_⊙, while the binary with P_ZAMS = 517.95 days has 8.27 M_⊙ of hydrogen envelope left and a radius of ∼546R_⊙. This remains until the remaining hydrogen is lost through stellar winds, which also removes all the angular momentum gained during the first mass transfer phase.

Fig. A.1. Evolution of stable reverse-mass-transfer systems with MZAMS, 1 ≈ 39.87 and q = 0.95 and different initial periods, including the first mass transfer and reverse phase. The rows from top to bottom show: the total angular momentum in the primary star, the radius of the primary star, and finally the relative Roche lobe filling factor of the primary and secondary. The colour map indicates the initial period at the zero-age main-sequence, with blue indicating the tightest and yellow the widest systems. The horizontal dashed black lines indicate when the star fills its Roche lobe.

With a mass ratio of q = 0.95, the secondary evolves off the main-sequence soon after. In the shortest period system, P_ZAMS ≈ 26.83 days, the secondary fills its Roche lobe while still on the main sequence, leading to a relatively prolonged and stable mass transfer phase. This results in the transfer of more mass and angular momentum compared to wider-orbit systems, which typically experience their initial reverse interaction upon core hydrogen depletion. In Figure A.2, we zoom in on the stable reverse mass transfer and exclude this short-period binary for clarity.

Fig. A.2. Stable-reverse-mass-transfer phase for the same models as in Figure A.1 with the same colour map, except for P = 26.83 days, which is not included. The models are aligned when starting the stable reverse-mass-transfer phase. The rows from top to bottom show the total angular momentum in the primary, the relative Roche lobe overflow factor of the secondary, the total hydrogen in the primary, and the fraction of critical rotation at the surface. In the top figure, the shaded areas indicate the reverse mass transfer.

In the top row of Figure A.2, we show the total angular momentum in the primary star and the duration of the stable reverse-mass-transfer phase as the shaded region. The shortest period systems (dark blue; ∼12000 yr) have a 6 times longer mass transfer than the widest periods (yellow; ∼5000 yr). Over time, angular momentum at the surface is able to diffuse inwards to spin up deeper layers in the star. For a longer duration mass transfer, the ongoing mass transfer will replenish the angular momentum of the outer layers of the accretor. As a result, the short-period systems are able to accrete more mass and angular momentum, as the total hydrogen in the third row in Figure A.2 shows, despite the surface being near critical rotation for the majority of the mass transfer. The amount of AM per accreted unit mass is constant around ∼7 g cm² s⁻¹ M_⊙. In Figure A.3, we show the total amount of angular momentum gained during the stable reverse-mass-transfer phase in blue. The total angular momentum gain in short-period binaries is higher than in long-period systems, but the angular momentum loss in the later evolutionary phases is higher for short-period binaries. Independent of the initial period, approximately 80% of the accreted angular momentum is lost between the stable reverse mass transfer and core carbon depletion.

With the evolutionary state of the primary being similar for different initial periods, the main difference is the duration of the stable reverse mass transfer leading to more mass accretion and angular momentum gain. A long-period binary has a shorter stable reverse-mass-transfer phase, and gains less angular momentum and mass. With a similar fraction of angular momentum lost across different initial periods, this puts most of the long-period primaries outside the regime for disk formation during collapse, while short-period primaries retain sufficient angular momentum to produce an accretion disk, despite the higher absolute angular momentum loss.

Fig. A.3. Amount of angular momentum (AM) gained during the stable reverse-mass-transfer phase (blue; dots) and the amount of AM lost between the end of the mass transfer and core carbon depletion (red, cross) for the models in Figure A.2. The shortest-period systems gain the most angular momentum, but also lose more angular momentum due to their large radial extension compared to long-period binaries. The fractional loss is approximately 80% of the AM gained at each initial period.

A.2. Initial primary mass dependence

The dependence on the initial primary mass is more straightforward compared to the period dependence. At the same period and mass ratio, a more massive primary gains more mass and angular momentum during the stable reverse-mass-transfer phase compared to a lower mass companion, as shown as the blue line in Figure A.4 for P_ZAMS = 19.31 days and q = 0.95 at Z = Z_⊙. At the same time, the interaction occurs earlier in the evolution of the primary due to its larger radius. As such, the duration between the reverse mass transfer and core carbon depletion is shorter for lower mass primaries compared to a higher mass primary, despite the longer evolutionary timescale. Together with a smaller radius and lower mass loss through stellar winds, this leads to less angular momentum loss in lower-mass primaries, as the red line in Figure A.4 shows. In conclusion, as the initial mass of the primary star increases in a binary, the angular momentum loss before core carbon depletion increases, which decreases the mass of an accretion disk at collapse.

Fig. A.4. Angular momentum gain (blue) and loss (red) during the reverse mass transfer and post-reverse mass transfer, respectively, as a function of initial donor mass for PZAMS = 19.31 days and q = 0.95 at Z = Z⊙. Only the binaries undergoing stable reverse mass transfer at this period are shown.

Appendix B: Star formation histories

Fig. B.1. Similar to Figure 6 but for the cumulative rate density from the Neijssel et al. (2019) star formation rate.

Since the star formation rate density and metallicity evolution impact the cosmic rate density of stable reverse-mass-transfer collapsar LGRB significantly, we perform the same potential energy budget calculations as discussed in Section 2 on the Madau & Fragos (2017) and Neijssel et al. (2019) star formation rate prescriptions. Their cumulative energy distributions are shown in Figures B.1 and B.2. Because they originate from the same underlying POSYDON stellar populations, the shape of each metallicity distribution is similar, but the star formation rate of each metallicity is redistributed. For example, Z_⊙ events have a total cumulative rate of ≃290 Gpc⁻³ yr⁻¹ in the Madau & Fragos (2017) population, but ≃600 Gpc⁻³ yr⁻¹ in the Neijssel et al. (2019) population. This is an effect of the latter prescription having a narrower metallicity distribution compared to Madau & Fragos (2017), which results in more high-metallicity star formation in the local Universe. Figures B.4 and B.3 show the decrease of the rate density from the energy limits. This effect is similar to that from the Illustris TNG simulation as shown in Figure 5.

Fig. B.2. Similar to Figure 6 but for the cumulative rate density from the Madau & Fragos (2017) star formation rate.

Fig. B.3. Stable reverse-mass-transfer LGRB rate density from the Madau & Fragos (2017) star formation history.

Fig. B.4. Stable reverse-mass-transfer LGRB rate density from the Neijssel et al. (2019) star formation history.

Appendix C: Free-fall collapse timescale

For each of the primary stars undergoing collapse through the stable reverse-mass-transfer channel, we calculate their free-fall timescale using

$$\begin{array}{c}\hfill t_{\mathrm{ff}}=\sqrt{\frac{3\pi }{32G\rho }},\end{array}$$(C.1)

where G is the gravitational constant and ρ is the average density of the star at collapse. We use the stellar radius and mass at carbon depletion to calculate ρ. In Figure C.1, the duration of each model is weighted by its contribution to the event rate at z = 0. The free-fall collapse timescale of the stable reverse-mass-transfer LGRB progenitors lies between ≃108 and ≃2204 seconds with the majority of durations around ≃120 seconds.

Observations use the T₉₀ duration to quantify the duration of the event, which is defined as the duration between the emission of 5% and 95% of the total measured counts of the LGRB. For LGRBs these peak around ∼50 seconds (see, for example Ghirlanda & Salvaterra 2022), but extend up to ∼1000 seconds.

The free-fall collapse durations of our models are at the higher end of the observed T₉₀ distribution, since we do not account for energies below the noise level nor do we consider only 90% of the counts. Additionally, the time for the accretion disk to form and the breakout time of the jet are not included in our calculation.

These effects would shorten the perceived duration of the collapse compared to the free-fall duration, although the infall via the disk can extend the duration for the fastest rotating models (see, for example Fryer et al. 2025). In general, the free-fall collapse durations of our LGRBs progenitors are of the same order of magnitude as the observed T₉₀ values for LGRBs.

Fig. C.1. Free-fall collapse timescales of the stable reverse-mass-transfer progenitors at carbon depletion. Each model is weighted by its contribution at z = 0. The red shading areas are stacked histograms for each metallicity with the solid black line indicating the total differential rate.

All Tables

All Figures

Fig. 1. Typical example of the evolution of a potential LGRB progenitor through the stable reverse-mass-transfer channel. The initial properties of the binary are Z = Z⊙, M1, ZAMS ≈ 25.1 M⊙, q = 0.95, and P ≈ 26.8 days. In all panels, red (blue) solid lines refer to the primary (secondary) star. The columns are split into five different evolutionary phases, focusing on the mass transfer phases. The second and fourth columns show the mass transfer and ‘reverse’ mass transfer phases, respectively. The first, third and last columns show the long, nuclear timescale, detached evolution of the system. Top two rows: Total mass (solid lines), He core mass (dashed lines) and CO core mass (dotted lines). The He (CO) core boundaries are set where the hydrogen (helium) fraction falls below 0.1. Third row: Mass transfer rate from Roche lobe overflow. The periods of Roche lobe overflow are marked by a red (blue) shading when the donor is the primary (secondary). Fourth row: Evolution of stellar radius (red: primary; blue: secondary) and separation of the system (black). Fifth row: Surface rotational angular frequency in units of the critical (i.e., mass shed) value. Sixth row: the dimensionless spin parameter cJ/GM2 of the star.

In the text

Fig. 2. Specific angular momentum distribution of the mass shells of the primary star in the example binary shown in Figure 1, at core carbon depletion in red. The kink at ≃7.8 M⊙ is the transition between the carbon-oxygen core and helium shell. It forms an accretion disk of 1.46 M⊙ during the collapse. For comparison, we plot the specific angular momentum of the ISCO, jISCO, for a non-rotating Schwarzshild BH (dashed black) and a maximally rotating Kerr BH (dotted black), as a function of BH mass. Additionally, we track the evolution of jISCO for the example primary during the collapse as its mass and spin change (dashed green). It tracks jISCO at the collapse of each mass shell and starts at 3 M⊙ with a BH of 2.5 M⊙ and 0.5 M⊙ neutrino mass loss.

In the text

Fig. 3. Four grid slices at Z⊙ with q = 0.85, 0.90, 0.95, and 0.99. Orange symbols represent systems where only the primary fills its Roche lobe throughout the evolution of the binary model. Green squares indicate systems where also the secondary fills its Roche lobe, but not at the same time as the primary. Purple symbols indicate systems where both stars fill their Roche lobe at the same time. The type of symbol shows the stability of the interaction: squares indicate stable mass transfer, while binary models reaching instability criteria are indicated with a diamond. Gray squares indicate no mass transfer, while black dots indicate a non-converged model or Roche lobe overflow when the model is initiated. For systems where the companion fills its Roche lobe and the interaction remains stable, we indicate if an accretion disk is formed during the collapse of the model. This is indicated by a coloured circle inside the square symbol, where the colour indicates the mass of the accretion disk.

In the text

	Fig. 4. Grids slices at q = 0.95 at Z = 10−3 Z⊙, 0.2 Z⊙, and 0.45 Z⊙. The same symbol and colours are used as in Figure 3. Less stable reverse mass transfer (green squares) occurs as the metallicity decreases and those that do, do not produce accretion disks upon their collapse.
In the text

Fig. 5. Redshift evolution of the stable reverse-mass-transfer LGRB event rate density for the IllustrisTNG metallicity and star formation history. The total potential population is indicated in black, while the dashed red line shows the intrinsic LGRB rate from Ghirlanda & Salvaterra (2022) with the best-fit parameters based on observations. The red shaded area is 10% to 20% of this rate. The five green shaded lines indicate the rate density for an energy limit of 1051 erg and different ηϕ.

In the text

	Fig. 6. Cumulative rate density (R[E > Ex]) within z < 2.5 from the IllustrisTNG simulation split up per metallicity with ηϕ = 10−3. The vertical line indicates a typical breakout energy of 1051 erg.
In the text

	Fig. 7. Total potential stable reverse-mass-transfer LGRB rate density over cosmic time from the IllustrisTNG simulation, Madau & Fragos (2017) with a log-normal metallicity distribution with σ = 0.5, and Neijssel et al. (2019) with a ln-normal distribution with σ = 0.39. No energy requirements are implemented.
In the text

	Fig. 8. Number of stable reverse-mass-transfer LGRBs per solar mass star formation (efficiency) as a function of metallicity. These metallicity bias functions are shown for the two Fryer et al. (2012) (rapid in purple; delayed in green) and the Patton & Sukhbold (2020) (orange) prescriptions.
In the text

Fig. A.1. Evolution of stable reverse-mass-transfer systems with MZAMS, 1 ≈ 39.87 and q = 0.95 and different initial periods, including the first mass transfer and reverse phase. The rows from top to bottom show: the total angular momentum in the primary star, the radius of the primary star, and finally the relative Roche lobe filling factor of the primary and secondary. The colour map indicates the initial period at the zero-age main-sequence, with blue indicating the tightest and yellow the widest systems. The horizontal dashed black lines indicate when the star fills its Roche lobe.

In the text

Fig. A.2. Stable-reverse-mass-transfer phase for the same models as in Figure A.1 with the same colour map, except for P = 26.83 days, which is not included. The models are aligned when starting the stable reverse-mass-transfer phase. The rows from top to bottom show the total angular momentum in the primary, the relative Roche lobe overflow factor of the secondary, the total hydrogen in the primary, and the fraction of critical rotation at the surface. In the top figure, the shaded areas indicate the reverse mass transfer.

In the text

Fig. A.3. Amount of angular momentum (AM) gained during the stable reverse-mass-transfer phase (blue; dots) and the amount of AM lost between the end of the mass transfer and core carbon depletion (red, cross) for the models in Figure A.2. The shortest-period systems gain the most angular momentum, but also lose more angular momentum due to their large radial extension compared to long-period binaries. The fractional loss is approximately 80% of the AM gained at each initial period.

In the text

	Fig. A.4. Angular momentum gain (blue) and loss (red) during the reverse mass transfer and post-reverse mass transfer, respectively, as a function of initial donor mass for PZAMS = 19.31 days and q = 0.95 at Z = Z⊙. Only the binaries undergoing stable reverse mass transfer at this period are shown.
In the text

	Fig. B.1. Similar to Figure 6 but for the cumulative rate density from the Neijssel et al. (2019) star formation rate.
In the text

	Fig. B.2. Similar to Figure 6 but for the cumulative rate density from the Madau & Fragos (2017) star formation rate.
In the text

	Fig. B.3. Stable reverse-mass-transfer LGRB rate density from the Madau & Fragos (2017) star formation history.
In the text

	Fig. B.4. Stable reverse-mass-transfer LGRB rate density from the Neijssel et al. (2019) star formation history.
In the text

	Fig. C.1. Free-fall collapse timescales of the stable reverse-mass-transfer progenitors at carbon depletion. Each model is weighted by its contribution at z = 0. The red shading areas are stacked histograms for each metallicity with the solid black line indicating the total differential rate.
In the text