Binary population synthesis models for core-collapse gamma-ray burst progenitors

ABSTRACT

Long-duration gamma-ray bursts (GRBs) are understood to be the final fate for a subset of massive, stripped envelope, rapidly rotating stars. Beyond this, our knowledge of the progenitor systems is limited. Using the Binary Population and Spectral Synthesis (bpass) stellar evolution models, we investigate the possibility that some massive stars in binaries can maintain the angular momentum required for jet production, while still loosing their outer envelope through winds or binary interactions. We find that a total hydrogen mass of M_H < 5 × 10⁻⁴ M_⊙ and a helium ejecta mass fraction of F_He < 0.20 provide the best thresholds for the supernova type II/Ibc and Ib/Ic divisions, respectively. Tidal interactions in binaries are accounted for by applying a tidal algorithm to post-process the stellar evolution models output by bpass. We show that the observed volumetric GRB rate evolution can be recreated using two distinct pathways and plausible distributions for burst parameters. In the first pathway, stars are spun up by mass accretion into a quasi-homogeneous state. In the second, tides maintain rotation where otherwise the star would spin-down. Both lead to type Ic supernova progenitors, and a metallicity distribution consistent with the GRB host galaxy population. The inferred core angular momentum threshold for jet production is consistent with theoretical requirements for collapsars, given the assumptions made in our model. We can therefore reproduce several aspects of core-collapse supernova/GRB observation and theory simultaneously. We discuss the predicted observable properties of GRB progenitors and their surviving companions.

methods: numerical, methods: statistical, stars: evolution, gamma-ray burst: general, stars: rotation, supernovae: general

1 INTRODUCTION

Long-duration gamma-ray bursts (GRBs) mark the evolutionary end points for a subset of massive, rapidly rotating, and stripped envelope (SE) stars. Their preference for low metallicity, star-forming host galaxies (Fruchter et al. 2006; Savaglio, Glazebrook & Le Borgne 2009; Levesque et al. 2010a; Perley et al. 2016b; Modjaz et al. 2019; Palmerio et al. 2019), association with energetic hydrogen and helium-poor type Ic-BL (broad line) supernovae (Galama et al. 1998; Iwamoto et al. 1998; Woosley, Eastman & Schmidt 1999; Nakamura et al. 2001; Stanek et al. 2003; Hjorth et al. 2003; Nomoto et al. 2004b, a; Modjaz et al. 2008, 2016; Cano et al. 2017), and theoretical considerations, all point to a core-collapse scenario. This could be a ‘collapsar’, (Woosley 1993; Woosley & MacFadyen 1999), in which a collapsing star powers a rapidly spinning black hole central engine, which launches relativistic jets. If the jets produced can escape the envelope (a requirement which might explain the association with SE supernovae, e.g. Modjaz et al. 2016), and we are aligned along the jet axis, the event is observable as a relativistically beamed GRB (MacFadyen & Woosley 1999; MacFadyen, Woosley & Heger 2001; Woosley & Bloom 2006; Hjorth 2013). Alternatively, the central engine could be a rapidly spinning neutron star with a strong magnetic field (a magnetar, Mazzali et al. 2014), which powers relativistic jets as it spins down.

Although the requirement for an SE progenitor with a large reserve of angular momentum at core-collapse is understood, the exact nature of the progenitors and the evolutionary pathways which lead to long-duration (core collapse) GRBs are not well constrained (see Levan et al. 2016; Fryer et al. 2019, and others for reviews). Models invoking single stars struggle to maintain enough angular momentum over the stellar lifetime, while simultaneously losing the outer envelope through winds (Vink, de Koter & Lamers 2001; Hirschi, Meynet & Maeder 2005; Vink & de Koter 2005; Woosley & Bloom 2006; Vink, Gräfener & Harries 2011; Modjaz et al. 2016).

One promising pathway is through chemically homogeneous evolution. In the later stages of binary evolution, the primary can expand and fill its Roche lobe, triggering accretion on to the secondary star. This can spin-up the secondary, and if sufficient angular momentum is transferred, rotational mixing can occur within the star (Cantiello et al. 2007; de Mink et al. 2009; Eldridge, Langer & Tout 2011; de Mink et al. 2013). If the star becomes fully mixed, it evolves chemically homogeneously, becoming smaller, hotter and hydrogen deficient, evolving blue-wards on the Hertzsprung–Russell diagram. If this process occurs at low metallicity, the secondary can retain the angular momentum gained from mass transfer, and end its life as a rapidly spinning type Ic progenitor – a good candidate for producing a GRB. Eldridge et al. (2019b) calculated the volumetric event rate of GRBs as a function of redshift, assuming that they arise solely from a quasi-homogeneous evolution (QHE) pathway (Yoon, Langer & Norman 2006; Cantiello et al. 2007), using the BPASS stellar evolution and population synthesis code (Eldridge et al. 2017). They found that although QHE stars could account for the observed GRB rate, there are significant uncertainties arising from luminosity function and beaming corrections. Additionally, BPASS only implements QHE pathways at metallicities Z < 0.004. GRBs are not exclusively found at such low metallicity, with observed examples arising in environments exceeding Solar abundance (e.g. Graham et al. 2009; Levesque et al. 2010b; Schady et al. 2015; Michałowski et al. 2018). QHE therefore cannot be the sole contributing pathway.

In this paper, we explore the possibility that massive stars in binaries, without undergoing QHE, can also produce GRBs given the right evolutionary history, even at moderate to high metallicities (see e.g. Podsiadlowski et al. 2004; Izzard, Ramirez-Ruiz & Tout 2004; Petrovic et al. 2005; Detmers et al. 2008; Trenti, Perna & Jimenez 2015). This route is particularly promising as binary interactions are likely responsible for the rapid rotation of massive stars (de Mink et al. 2013). To investigate, we introduce a tidal interaction post-processing algorithm to the bpass binary stellar evolution models in order to account for angular momentum transfer between orbital and stellar rotation. We investigate subsets of the resultant high mass, rapidly spinning, SE population, and explore the consequences of matching synthetic GRB rates arising from our selection criteria to observations. The paper is structured as follows: In Section 2, we outline the categorization of models into different classes of core-collapse progenitor, producing updated BPASS predictions for the properties of supernova progenitors. Section 3 describes our implementation of tides within BPASS, and the effect of their inclusion. A comparison to observations is presented in Section 4, with a Bayesian analysis used to infer the most likely parameter values in our two-progenitor pathway model. Section 5 details the results, including predictions for progenitor properties such as their position on the Hertzsprung–Russell diagram, and the interior angular momentum at core collapse. A discussion of the physical interpretation of these models follows in Section 6, with conclusions presented in Section 7. Where required, a flat ΛCDM cosmology with h  = 0.7, Ω_M  = 0.3, and Ω_Λ  = 0.7 is used.

2 CLASSIFICATION OF CORE-COLLAPSE PROGENITOR MODELS

GRBs, when close enough (typically within ∼1 Gpc) and followed up with prompt and deep optical–NIR observations, typically show evidence for an associated broad-lined type-Ic supernova (Ic-BL SN). The preference for stripped-envelope supernovae of this type is likely a jet escape requirement, as an envelope of helium or hydrogen around the star may prove an insurmountable barrier to the ejection of fast-moving material and hence stifle any incipient jet. It appears as though not every Ic-BL SN has an accompanying GRB, although in some cases the jet might be choked, as may have been the case in SN 2008D (Mazzali et al. 2008; Piran et al. 2013; Modjaz et al. 2016; Sobacchi et al. 2017; Barnes et al. 2018). Nevertheless, the minimum requirement for an SE is an observational constraint on the progenitor stars of GRBs. Using the binary stellar evolution models of bpass (v2.2.1, Stanway & Eldridge 2018), the first step in identifying GRB pathways is to find models which likely produce type Ic-SNe.

At each step in the BPASS models, stellar properties such as mass and chemical composition are calculated. When core carbon burning ends, the stellar model is assumed to progress rapidly through the final stages of core burning (if massive) before ending its life (Woosley, Heger & Weaver 2002). If a supernova occurs, the stellar model is then replaced by a neutron star or black hole based on the final stellar mass and the amount of mass that could be removed by energy injection from a supernova of 10⁵¹ erg. Because supernovae are observationally classified according to their spectral properties, which reflect the chemical composition of the ejecta and circumstellar medium (Filippenko 1997), we can predict which supernova type each model will produce based on its chemical properties immediately before core collapse.

Our first task is to identify stars which produce supernovae. The fate of stars with final masses |$\lt 2\, \mathrm{M}_{\odot }$| is uncertain. They might produce weak, faint, and fast core-collapse supernovae (Moriya & Eldridge 2016), or, depending on the central carbon abundance, type Ia supernovae. Alternatively, they could form white dwarfs without a supernova. Supernovae from these stars would in any case be rare events (Eldridge et al. 2017).

We find models that likely experience a supernova by two methods. First, if at the end point of our model the CO core mass is greater than |$3\, \mathrm{M}_{\odot }$|⁠, a supernova occurs. Secondly, at lower masses, due to the complexity of second dredge-up and carbon burning in super-AGB stars, we require that the final total mass exceeds |$2\, \mathrm{M}_{\odot }$|⁠, the CO core mass is greater than |$1.38\, \mathrm{M}_{\odot }$| and core carbon burning has occurred (when the ONe core is great than |$\gt 0.1\, \mathrm{M}_{\odot }$|⁠). Then, if the remnant mass is in the range ∼1.4–3|$\, \mathrm{M}_{\odot }$|⁠, a neutron star results. Above this remnant mass, black holes are created, and supernovae in this regime are uncertain, with ‘islands of explodability’ in mass (e.g. Sukhbold et al. 2016; Auchettl et al. 2019; Sukhbold & Adams 2019; Woosley 2019). It is believed that most stars above a zero-age main-sequence mass of |$20\, \mathrm{M}_{\odot }$| will create a black hole remnant and have their explosion engulfed within the Schwarzschild radius, thus producing no visible display. This picture is increasingly backed up by a lack of high-mass stars seen in pre-explosion imaging (Smartt 2009, 2015), possible tension between star formation and SN rates (Horiuchi et al. 2011), and the direct observation of a quietly disappearing massive star (Adams et al. 2017). As such, we initially disregard any model which produces a remnant mass greater than |$3\, \mathrm{M}_{\odot }$| (i.e. a black hole) as a vanishing event (see Eldridge & Tout 2004, for a description of how remnant and ejecta masses are calculated).

Otherwise, if the stellar structure at the end of the model still satisfies the carbon–oxygen, oxygen–neon, and total mass conditions listed above, we assume that a visible supernova will be produced. We then split these models into three categories, types Ib, Ic, and II, based on their chemical composition (Eldridge et al. 2011, 2013, 2017).

The type I/type II divide is defined by the total mass of hydrogen in the star; observationally this is probed by the presence of hydrogen lines in the SN spectrum, or lack thereof. Theoretically, a small residue of hydrogen may still be present without leaving observational evidence, as long as it is below some threshold. We vary this threshold, labelling models with a total hydrogen mass less than M_H,tot as type I and those with more as type II. The type Is are further divided into Ibs (which have visible helium lines) and Ics (with no helium or hydrogen in their spectra). This is decided by the fraction of helium in the ejecta, F_He,ejecta, and this threshold is also varied. The total hydrogen mass is used, rather than the fraction of hydrogen in the ejecta, because (unless QHE is underway) the majority of the hydrogen should be on the surface of the star and therefore mixed into the ejecta upon supernova.

The models in each category are then assigned a weighting, which corresponds to the estimated occurrence frequency of each model in a star formation episode, assuming |$10^{6}\, \mathrm{M}_{\odot }$| of metal enriched gas were allowed to collapse and form stars with 100 per cent efficiency. The weightings are informed by observations of stars in terms of binary fractions, mass distributions, binary orbital periods, and mass ratios (Moe & Di Stefano 2017) and correspond to the weightings from v2.2 of the bpass spectral synthesis models (Stanway & Eldridge 2018). We use the standard bpass initial mass function (IMF). This is parametrized as a power law in |$\frac{\mathrm{ d}N}{\mathrm{ d}m}$|⁠, and has a slope of −1 in the range 0.1–0.5 M_⊙, and −2.35 from 0.5 to a maximum of 300 M_⊙.

For a |$10^{6}\, \mathrm{M}_{\odot }$| stellar population, at each metallicity Z (12 over the range 10⁻⁴ to 0.04 by metal mass fraction) we calculate the relative contribution to the total number of supernovae from each type (II, Ib, and Ic), summed over 10 Gyr from t = 0. The average of the fractions over the Z = 0.008, 0.010, 0.014, and 0.020 metallicities are then compared to observed fractions from volume limited supernova studies. The metallicity range used for comparison is broadly representative of the range seen in the local Universe. For the type I/type II ratio we refer to Eldridge et al. (2013), which is a volume-limited sample and therefore representative of the intrinsic ratio. Shivvers et al. (2017) is used for the Ib/Ic ratio, due to their careful reclassification of a number of SE SNe from the Lick Observatory Supernova Search (LOSS).

In using the LOSS data, we have classified type IIb events as Hydrogen-rich (so that they are in the type II category), and we have grouped regular, peculiar, and broad-line Ic SNe together. We also subdivide the type IIs into IIP (plateau) and II-other classes, but do not use this for constraining the pre-core-collapse properties as we would have too many degrees of freedom. Instead, we use the selection criteria of Eldridge et al. (2017), where type IIPs are defined to have |$M_\mathrm{H,tot} \gt 1.5\, \mathrm{M}_{\odot }$| and a hydrogen to helium ratio greater than 1.05.

We use the difference between the bpass weighted fractions and the observed fractions to calculate χ², which we then minimize by varying M_H,tot and F_He,ejecta across parameter space. We show the results of this process in Fig. 1. We find that the best-fitting values are M_H,tot  = |$5\times 10^{-4}\, \mathrm{M}_{\odot }$| and F_He,ejecta = 0.20, with SN types II, Ib, and Ic making up 75 per cent, 14 per cent, and 10 per cent of the local Universe supernova rate, respectively. This compares to the Eldridge et al. (2013) and Shivvers et al. (2017) combined percentages of 74 ± 7.6, 14.5 ± 2.6, and 11.5 ± 2.0. The differences correspond to a χ², equivalent to a reduced χ² in this case, of 0.27. The large uncertainties on the observed fractions, taken together with the large number of free parameters, allow for a range of values above these best-fitting thresholds. Hydrogen masses up to 0.0022, and F_He,ejecta values in the range 0.17–0.26, are permitted within the 68 per cent confidence interval. None the less, there is a clear preference for lower values, with the upper end of the explored parameter space firmly excluded.

Figure 1.

The result of matching type II, Ib, and Ic SN relative fractions to observations, by minimizing χ² over M_H,tot and F_He,ejecta parameter space. The best-fitting classification thresholds, assuming that black hole producing events do not go supernova, are marked by a cross. Solid contours represent the 68 and 90 per cent confidence regions for one degree of freedom.

Dessart et al. (2011) predicted the supernova spectra from Wolf–Rayet stars, and found that hydrogen lines are initially visible if more than |$10^{-3}\, \mathrm{M}_{\odot }$| of hydrogen is present, likely making such events type IIb. This is consistent with our M_H,tot  = |$5\times 10^{-4}\, \mathrm{M}_{\odot }$| type I/II threshold. Dessart et al. (2012a) then showed that helium ejecta fractions as low as 0.2–0.3 can produce visible helium lines in the spectrum, again consistent with our results.

Using the selection criteria determined from comparison to the local Universe supernova rates, we extend these cuts over the entire BPASS metallicity range. The predicted SE to hydrogen-rich ratio as a function of metallicity is shown in Fig. 2, with the same ratio as reported by various observational studies also shown. The bpass metallicity scaling used in this paper includes modifications to the previous scale used (Xiao et al. 2018), to improve agreement with the data at low metallicities (see Fig. 44 of Eldridge et al. 2017, and Table 1). We note that the stellar evolution is driven primarily by iron abundance and bulk metallicity mass fraction, rather than the oxygen abundance. The bpass relative rates generally agree with observation, including with those studies listed in (Smartt 2009), which are not shown on the figure for clarity. However, at the highest metallicities, there are discrepancies. Our models appear to overpredict the rate of SE events relative to Hydrogen-rich events at significantly super-Solar metallicities. We caution that in this regime the remnant mass and ejecta composition are strongly sensitive to the assumed mass-loss rates on the asymptotic giant branch and immediately preceding the explosion. Small adjustments in the assumed wind prescription may have a significant effect on the supernova type ratio and further work is required to explore this.

Figure 2.

The ratio of SE to hydrogen-rich (type II) supernovae as a function of metallicity, using our best-fitting M_H,tot and F_He,ejecta thresholds (black line), assuming that supernovae only occur if a black hole is not produced. Metallicites in the range 0.008–0.02 by mass fraction are indicated by a grey shaded band. The red shaded region is a 0.1 dex uncertainty, intended to give an indication of the discrepancies between different metallicity scales (Eldridge et al. 2017). Where a sample metallicity is not explicitly defined, we assume that it samples the local metal mass fraction spread of 0.008–0.02.

Using our classification criteria, we place every pre-core-collapse model (in the metallicity range 0.008–0.020) on a Hertzsprung–Russell diagram in Fig. 3. We also show a range of known or candidate progenitors, identified from pre-explosion imaging. In most cases, there is good agreement. The only candidate which could be in tension is that of the type Ic progenitor, SN 2017ein (Kilpatrick et al. 2018; Van Dyk et al. 2018), although it is currently unclear how much of the pre-explosion stellar light was from the progenitor, rather than a possible binary companion, and so this data represent an upper limit. We note that increasing the remnant threshold for supernova production (i.e. assuming lower supernova mass-ejection energies or less efficient black hole stifling of any core-collapse event) would improve the agreement here. The physical properties of observed Wolf–Rayet stars are similarly more typical of the most luminous Ic progenitor models in our classification, with luminosities around log₁₀(L/L_⊙) ∼ 5.6 (Neugent & Massey 2019), although we have no constraint on whether these stars will go supernova or not.

Figure 3.

Following the categorization of bpass model such that the observed supernova relative rates are reproduced, we place these models on the HR diagram (in orange) according to core-collapse type. The shading represents the number density of models. Observed progenitors (in black) include the vanishing star N6946-BH1 (marked by a cross, Adams et al. 2017), for which we have adjusted the luminosity to reflect an updated distance estimate (Eldridge & Xiao 2019), SN 1987A (marked by a star in the IIP section, Walborn et al. 1987), several IIP SNe from Smartt (2015), type IIb and IIL SNe (both classified as II-other, and marked by dots and crosses, respectively, Smartt 2015), SN 1993J (also a IIb event, Aldering, Humphreys & Richmond 1994), the Ib progenitor of SN iPTF13bvn (Eldridge & Maund 2016) and finally the candidate Ic progenitor of SN 2017ein (Van Dyk et al. 2018). The marked luminosity for 2017ein is an upper estimate, assuming that it was a single star (Kilpatrick et al. 2018). This progenitor system is also indicated on the ‘vanishing’ panel (in purple) for reference.

A possible resolution to the type Ic tension is that 2017ein may have produced a supernova and a black hole remnant. Modelling of pre core-collapse stellar structure shows that some black hole producing events can successfully explode if the core has a low compactness parameter (Sukhbold et al. 2016), although this is just one of several relevant quantities (Mapelli et al. 2019). Furthermore, work in this area has typically focused on red supergiants, since these are the class of stars which we would have expected to see in pre-explosion imaging – SE supernovae are typically too distant to expect to see the progenitor. Therefore, we have few constraints on the explodability of such stars, and the generic vanishing threshold we have applied might not be applicable to SE stars. Another possibility is that supernovae can occur in stars that would otherwise implode thanks to energy injection from a rotationally powered central engine. This could be from disc winds (e.g. MacFadyen & Woosley 1999) or a jet (for example, type Ic-BL SNe might harbour choked jets, Mazzali et al. 2008; Piran et al. 2013; Modjaz et al. 2016; Izzo et al. 2019). We refer the reader to Piran et al. (2019) for a review of these ideas.

Regardless of the mechanism, it is likely that some black hole producing events cause visible supernovae, and that these have been classified as vanishing in our analysis so far. One example are pair-instability supernovae (PISNe; Heger & Woosley 2002), which we classify as occurring in stars with helium cores between 64 and 133 |$\, \mathrm{M}_{\odot }$| at core collapse. However, these particular events likely disrupt the star leaving no remnant. They are also exceptionally rare and do not significantly impact the analysis thus far.

It is likely that GRBs do indeed require such a black hole central engine, but these events are also luminous and associated with supernovae. In order to decide which BPASS models end their lives as GRBs, we therefore need to consider the rotation of stars which are Ic-like in their chemical composition, but which also produce black holes at core collapse. We note that these additional black hole producing SNe are rare, given that they have higher mass progenitors. They make a small difference to the relative SN rates in the local Universe, therefore having a minimal impact on our supernova categorization – this is demonstrated later, see for example Fig. 10.

3 INTRODUCING COMPLEX TIDAL INTERACTIONS TO BPASS

Surface composition and structure of the pre-explosion progenitor allow us to identify the subset of BPASS models which produce Ic SNe and black holes. In order to produce a GRB, however, rapid rotation is also required. bpass accounts for the first-order effects of mass-loss and mass transfer on the spins of stars in binaries, but we must also consider the transfer of angular momentum between the orbit and stellar spins due to tidal interactions. bpass currently only invokes tidal interactions when Roche lobe overflow occurs (Eldridge et al. 2017), however, tides may have a substantial impact at all stages of stellar evolution. In this section we describe our first-order application of tides to the bpass output models.

3.1 Radiative damping of the dynamical tide

There are two mechanisms for tidal dissipation within stars. These are convective damping of the equilibrium tide, and radiative damping of the dynamical tide (Zahn 1975, 1977; Hut 1981; Goldreich & Nicholson 1989; Hurley, Tout & Pols 2002). Tidal forces are strongly dependent on radius, whatever the dissipation mechanism, and so the nature of the stellar envelope (rather than the core) determines which mechanism is dominant. Given that we are interested in stars with masses greater than |$2\, \mathrm{M}_{\odot }$|⁠, which have radiative envelopes for the majority of their lives, we employ radiative damping of the dynamical tide and assume that convective damping is negligible.

In the radiative regime, the mass and orbit of a companion star introduces a time-varying external gravitational potential. This variation couples to g modes (where buoyancy is the restoring force) within the radiative envelope. The density within this zone is not constant, and the frequencies corresponding to orbital motion preferentially drive g modes deeper within the radiative zone, near the convective core, because the Brunt–Väisälä oscillation frequency is density dependent. The induced g modes distort the star, allowing gravitational torques to act. The excited waves would be standing waves if they were reflected at the stellar surface, but because the radiative time-scale there is short, they are only partially reflected, undergoing a phase shift. Angular momentum is therefore transported from the core to the surface, or vice versa, by the induced g modes (Zahn 1975; Goldreich & Nicholson 1989).

3.2 Implementation of tides

We do not recalculate the bpass models to include tides at each time-step, as this would be extremely computationally expensive. Instead, our approach is analogous to the rapid population synthesis models of Hurley et al. (2002). We choose a subset of the bpass output models and implement a post-processing algorithm as described below. First, we identify every primary and secondary star model which produces a black hole remnant and has the chemical properties of a Ic progenitor immediately before core collapse. In this model subset, the lowest initial mass, at any metallicity, is |$10\, \mathrm{M}_{\odot }$|⁠. Using this to inform our strategy, we opt to consider every primary and secondary model in bpass with Zero Age Main Sequence (ZAMS) masses greater than |$7\, \mathrm{M}_{\odot }$|⁠, expecting that binary interaction, enhanced by tides, might cause some slightly lower mass stars to move on to black hole/type-Ic pathways.

For each model, over each time-step, we calculate the expected change in the orbital semimajor axis Δa_tides due to tides. The change in a and corresponding change in the stellar rotational angular velocity Ω are given by

$$\begin{eqnarray*} {\Delta }{a_\mathrm{tides}}=\frac{6ka}{t_\mathrm{damp}}q(1+q)\Big (\frac{R}{a}\Big)^{8}\Big (\frac{\Omega }{\omega }-1\Big), \end{eqnarray*}$$

(1)

and

$$\begin{eqnarray*} {\Delta }{\Omega }=-\frac{3k}{t_\mathrm{damp}}\frac{q^{2}}{r_{g}^{2}}\Big (\frac{R}{a}\Big)^{6}(\Omega -\omega) , \end{eqnarray*}$$

(2)

where k is an apsidal motion constant (typically in the range 0.01–0.1, we adopt 0.05, Zahn 1975), q is the ratio of secondary mass M to primary mass, R is the stellar radius, ω is the orbital angular velocity, and r_g = I^0.5M^−0.5R⁻¹ (Zahn 1977; Hut 1981; Hurley et al. 2002; Eldridge & Tout 2019). The moment of inertia I is assumed to be that of a solid sphere, |$\frac{2}{5}MR^{2}$|⁠, using the total stellar mass M (and assuming solid body rotation). This approximation is made for the current study as detailed calculations using stellar structure models would be computationally expensive. Most stars show only a small differential rotation gradient while still on the main sequence, although post-main-sequence stars may show a significant disconnect between envelope and core rotation (Heger, Langer & Woosley 2000). The final term left undefined is the damping time-scale, t_damp. This is the time-scale on which tidally induced distortions dissipate. For radiative damping, it can be shown that the damping time-scale is

$$\begin{eqnarray*} t_\mathrm{damp} = \frac{ka^{5}}{(1.9782\times 10^{4})\, MR^{2}\, (1+q)^{5/6}\, E_{2}}\, \mathrm{yr}, \end{eqnarray*}$$

(3)

where M, R, and a are in Solar units (Hurley et al. 2002), and E₂ is a second-order dimensionless coefficient, parametrizing the strength of the tidal interaction given the internal structure of the star. E₂ was calculated for a limited number of stellar models by Zahn (1975). Subsequently, Hurley et al. (2002) fitted a functional form to these values, producing a prescription for E₂ as a function of mass. We note that Kushnir et al. (2017) significantly improves upon this, however their implementation relies on knowledge of the radial structure of the star – information which is not readily available in the standard bpass output files, and would again make this process computationally expensive. A full implementation would require recalculating 250 000 individual detailed bpass stellar evolution models, which is far beyond the scope of this paper. Comparing a handful of E₂ parameters calculated by Kushnir et al. (2017) to the form assumed by Hurley et al. (2002) gives results that are in broad agreement.

Observations of massive stars, for example the VLT-FLAMES Tarantula survey (Ramírez-Agudelo et al. 2015, 2017), suggest that they are born spinning with initial stellar rotational velocities typically around 30–40 per cent of their critical break-up velocity. Dufton et al. (2019) construct a probability density function for the rotational velocities of (apparently) single O stars in the Tarantula Nebula (Large Magellanic Cloud, LMC) and NGC 346 (Small Magellanic Cloud) and again find a median rotation around 40 per cent of critical. By contrast, Stevance et al. (2018) used spectropolarimetry of galactic WO stars to show that their rotational velocities are less than 10 per cent of critical, however, these stars are unambiguously at a late stage of evolution, thus suggesting that the spins of these massive stars evolve significantly during their lifetime. However, in low metallicity environments, the spin-down is likely reduced. Vink et al. (2011) and Vink & Harries (2017) report the detection of young Galactic and LMC Wolf–Rayet stars which have surface rotations that could be conducive to GRB production.

For primary models, we adopt an initial equatorial rotational velocity that is 0.4 of the critical rotation. This is faster than previously assumed in bpass and its predecessors. The current prescription, first used by Hurley, Pols & Tout (2000), gives early-type stars (specifically, those that end their lives as black holes) initial rotational velocities which are 10–20 per cent of the equatorial break-up velocity. Our ansatz is much faster, but consistent with the observations described above.

We add the tidally driven change in orbital separation, Δa_tides, to the orbital change due to stellar evolution processes which have already been accounted for, such as mass-loss, so that Δa_total = Δa_tides + Δa_BPASS. The bpass and tidal changes can act against each other. The orbital evolution from this point forward may correspond more closely to that of a bpass model at the same stellar evolution stage which began life with a slightly different orbital separation. If the new orbit is now closer to a different bpass model (with the same stellar masses, at the closest matching time-step), we switch models and update the parameters, continuing from there. We let the semimajor axis a, orbital angular velocity ω, and stellar rotational angular velocity Ω vary smoothly and do not update these numbers when a model change occurs. Otherwise, we stay on the same model.

This process is continued until a significant event happens in the system – this is defined as either a common envelope phase or Roche lobe overflow. At this point, we synchronize the primary spin to the orbital period and maintain this synchronization until core collapse. Because tidal forces are strongly dependent on radius, any post-main-sequence expansion is likely to result in synchronization (Hurley et al. 2002). Since we assume that most high-mass systems are synchronized before the end of the primary model, we start each secondary model in a synchronized state, where the rotational angular velocity of the secondary equals the orbital angular velocity of the compact object companion. The percentage of model outcomes changed by the influence of tides is around ∼5 per cent, out of ∼7000 models considered per metallicity. Most of these changes do not alter the final event categorization, as many models are shifted on to a slightly different evolutionary track which still results in the same class of transient.

An example of the methodology is shown in Fig. 4, demonstrating the process of jumping between bpass models when tidal interactions significantly alter the binary orbit. The synchronization of the binary is also shown when the common envelope phase occurs at t ∼ 35 Myr.

Figure 4.

An example of our methodology. The blue dashed line indicates the orbital period evolution of the initial model (labelled Model A). The red dashed line shows the orbital period of the model at the time of core collapse (Model B). Other models can be passed through in between these. The solid green line is the primary star’s rotational angular velocity, which starts at 40 per cent of critical, and the solid black line is the orbital angular velocity. A vertical grey line indicates when the model switch occurs. The system in this example is at a metallicity of Z = 0.014 and consists of a M_ZAMS = 8 |$\, \mathrm{M}_{\odot }$| star with a 7.2 |$\, \mathrm{M}_{\odot }$| companion, and an initial orbital period of log₁₀(P/d) = 0.2. The system starts in a state where the primary is spinning faster than the orbit of the companion, and so the system moves apart as angular momentum is transferred. Tidal forces get stronger as the primary expands, leading to synchronization at ∼35 Myr. In this case, the added orbital angular momentum gained from tidal interactions is sufficient to eject the common envelope, and consequently the binary moves apart. Without tides, the system would retain more mass, shortening the primary lifetime, such that it explodes in <40 Myr without significant evolution in orbital period.

3.3 The effect of tides on high-mass stars

We compare the number of models which now produce either Ic SNe with a neutron star remnant, or direct collapse to a black hole, to the previous numbers from before tides were considered. We find that models are both added to, and removed from, the pool of models which produce Ic SNe and black holes. The overall effect is shown in Fig. 5. For a more complete breakdown of the changes due to tides, refer to Fig. A1 of the online Appendix, which show the number of high-mass stellar death events (Ib, Ic, PISN, and vanishing) per |$10^{6}\, \mathrm{M}_{\odot }$| of star formation for each metallicity and bpass model type, both before and after tides. The effect of tides is complex and depends on the specific system in question. In most systems, given our initial rotational velocities, the binary is pushed apart, leading to reduced envelope mass-loss from the primary (and a shorter main-sequence lifetime). However, this same transfer of angular momentum to the orbit can make the ejection of a common envelope more efficient. These effects are roughly balanced, such that the overall change in the number of SE progenitors is small, as Fig. 5 shows.

Figure 5.

The change in the number of Ic SNe and black holes produced per |$10^{6}\, \mathrm{M}_{\odot }$| of star formation, as a function of metallicity, due to the inclusion of tidal interactions.

In order for a stripped-envelope star to produce a GRB, it must also launch jets. A key parameter of interest in jet production is the angular momentum of the star at core collapse. Woosley (1993) and Woosley & MacFadyen (1999) calculated that the specific angular momentum of material just outside the newly formed black hole should be >∼10¹⁶ cm²s⁻¹ in order for accretion to occur, otherwise material directly infalls and energy extraction is inefficient.

In Fig. 6, we visualize the effect of tides on the final angular momentum of Ic-SN and black hole progenitors. For the evolution without tides, the spin of the star in question is synchronized to the orbit in the final time-step in order to obtain an estimate of the specific angular momentum j. The resultant distributions of j, with and without tidal evolution over the stellar lifetime, are binned, and the difference in the normalized fraction contributing to each j bin is shown. Notably, it is the high metallicity models which are predominantly affected by the inclusion of tides. The evolution of these systems was previously dominated by wind-driven mass-loss, however, tidal interactions are acting to maintain angular momentum in the primary when it would otherwise be lost. Although tides typically push binaries apart over the main-sequence lifetime, they tend to end their lives spinning more rapidly due to synchronization in the mass transfer and common envelope phases.

Figure 6.

The difference in each specific angular momentum bin between the number of stellar models producing Ic-SNe and black holes at core collapse, before and after tides are considered. Metallicity is indicated by the colour gradient. The y-axis represents the bin heights for the number of models after tides, subtracted from the bins heights before tides were considered, where the histograms are normalized so that the enclosed area of each is equal to unity. This therefore provides an indication of changes in the j distribution due to tides.

Fig. 6 also shows a high-j spike due to tides, which is mostly populated by low-metallicity systems. Inspecting the mass distribution of these high angular momentum stars, they are typically very high mass. Our assumed IMF in this work allows for ZAMS masses up to |$300\, \mathrm{M}_{\odot }$|⁠. Because tidal forces are extremely radius sensitive (proportional to the eighth power in R in the prescription used), and because low-Z stars are more likely to maintain very high masses all the way to core collapse, the highest j bins are naturally populated by the lowest Z stars.

Another effect of introducing tides to BPASS is that the delay times, the intervals between star formation and core-collapse, are changed. In general, over the main-sequence lifetime, tides act to push binaries apart, therefore decreasing the mass lost to a companion. This means that stars leave the main sequence with higher masses than before tides were considered, which typically shortens lifetimes and therefore delay times. However, if the stars are further apart, the orbital angular momentum is greater and common envelopes are more efficiently ejected, therefore reducing the system mass (and increasing delay times). These two effects, which act to add and remove mass, respectively, dominate equally as often. In Fig. 7, as in Fig. 6, we show the change of the contribution to each delay time bin from the Ic-SN and black hole progenitor population, before and after tides are included. The overall effect of tides on the delay time distributions (DTDs) is a small shift to shorter delays.

Figure 7.

The change in the number of models falling in each delay time bin, for black hole and Ic-SN producing models. Delay times are measured from the zero age main sequence. The overall effect of tides on the DTD of massive stars is small, since model lifetimes are both shortened and lengthened by decreased mass transfer and more efficient common envelope ejection, respectively. The y-axis is calculated in the same way as Fig. 6.

Because the application of this tidal algorithm changes the frequency with which primary models end in different configurations, the weightings of the secondary models are changed. These are calculated based on the primary model endpoints and follow the evolution of post-supernova systems to later time-steps. Using the IMF and initial binary population parameters of the primary models, and the associated new end points for these binaries, the secondary models are re-weighted to account for the tidal evolution. In subsequent analysis, these new weightings are used.

4 CALCULATION OF LONG GAMMA-RAY BURST VOLUMETRIC EVENT RATES AND PROGENITOR PROPERTIES

4.1 Transient population synthesis

The search for additional, tidally induced GRB pathways (referred to hereafter as ‘tidal GRBs’) is motivated by the observed metallicity distribution of GRB progenitor environments (see e.g. Trenti et al. 2015; Palmerio et al. 2019), but the addition of any such pathway must also be consistent with observed long GRB volumetric event rates. For the existing QHE pathway, any stars with a remnant mass |$\gt 3\, \mathrm{M}_{\odot }$|⁠, which accrete >5 per cent of their initial mass at Z < 0.004, fulfill our composition and spin requirements for a GRB (Eldridge et al. 2019b). We have now further identified the bpass models which are rapidly spinning at core-collapse, at all metallicities, and which would produce a type Ic-SN if collapse to a black hole did not hamper the escape of emission from an explosion. We expect a currently undefined subset of these to also launch GRBs and hypothesize that the principle determinant of whether they do so is whether they exceed a threshold in specific angular momentum, j_cut. We will assume that these two pathways are the only contributors to the GRB population and now turn our attention to deciding which of the tidally spun models produce GRBs, and whether we can re-create the observed GRB rate and its evolution over cosmic history using plausible selection criteria.

We vary the specific angular momentum threshold to define sets of plausible GRB progenitors. The model weightings are used to determine, for a given mass of stars formed at metallicity Z, how many GRBs we would expect. We also know the delay times of these models. In order to construct the GRB rate as a function of redshift, we follow the methodology of Eldridge et al. (2019b) to apply a metallicity-dependent cosmic star formation rate history. The star formation rate density (SFRD; Madau & Dickinson 2014) as a function of redshift is given by

$$\begin{eqnarray*} \psi (z) = 0.015\frac{(1+z)^{2.7}}{1+((1+z)/2.9)^{5.6}} \, \mathrm{M}_{\odot }\,\mathrm{yr}^{-1}\,\mathrm{Mpc}^{-3}, \end{eqnarray*}$$

(4)

and this used in conjunction with the formalism of Langer & Norman (2006), which decomposes the SFRD into contributions from different metallicities as a function of redshift. The SFRD at or below a metallicity Z, at redshift z, is given by

$$\begin{eqnarray*} \Psi \left(z,\frac{Z}{\mathrm{Z}_{\odot }}\right) = \psi (z) \frac{\hat{\Gamma }[0.84,(\mathrm{\mathit{ Z}/Z}_{\odot })^{2}10^{0.3z}]}{\Gamma (0.84)}, \end{eqnarray*}$$

(5)

where |$\hat{\Gamma }$| and Γ are the incomplete and complete Gamma functions. We use 0.020 for the Solar metallicity |$\, \mathrm{Z}_{\odot }$|⁠, which corresponds to 12 + log(O/H) = 8.93 in the abundance distribution pattern adopted by bpass (see Eldridge et al. 2017).

Low metallicities dominate star formation only in the very early Universe, with the peak of Solar-metallicity star formation at around z ∼ 1, and the overall peak in star formation at z ∼ 2. By applying the GRB progenitor model rate per |$10^{6}\, \mathrm{M}_{\odot }$| to the SFR at each redshift and metallicity, and accounting for the delay times, we can construct the intrinsic (and hence estimate the observed) GRB rate. This requires, however, the application of an angular momentum cut-off j_cut for the tidal GRB pathways, and rate corrections due to jet beaming and the GRB luminosity function.

4.2 Bayesian parameter estimation

In order to investigate whether our two-pathway model is plausible, we allow four model parameters to vary, and infer their values by comparison to observational data. The parameters are,

• the jet half-opening angle θ,

• the lower limit on the isotropic equivalent energy of GRBs, E_low,

• the specific angular momentum j_cut,⊙, at Solar metallicity, above which GRBs can occur in our black hole producing, SE, rapidly rotating models,

• the index n, which allows for a metallicity dependence of this angular momentum cut (as j_cut ∝ j_cut,⊙(Z/Z_⊙)ⁿ).

The first two parameters determine whether a GRB is likely to be seen in the observed sample, while the latter two determine whether an event is likely to launch a jet at all. For observed GRB rates, we use the Swift Gamma-Ray Burst Host Galaxy Legacy Survey (SHOALS) sample (Perley et al. 2016a). This is the largest unbiased sample of GRBs with identified hosts (and hence redshift and metallicity estimates) consisting of bursts with isotropic equivalent energies, E_iso, greater than 10⁵¹ erg.

Of the four model parameters, two are corrective and change the number density of events equally across redshift (θ, E_low). The other two (j_cut, n) may also affect how the rate varies with redshift. The half-opening angle θ is used to account for the fact that GRBs are strongly beamed, and therefore most events are seen off-axis and not detected, so that the observed rate is much lower than the intrinsic one. This correction is given by [1 − cos θ]⁻¹ for bipolar jets, and assumes that within a viewing angle θ, we are equally likely to detect the burst whatever the orientation. The second factor corrects for the GRB luminosity function (actually an isotopic-equivalent energy E_iso function). This is required because the SHOALS comparison data are comprized exclusively of high-energy bursts, above 10⁵¹ erg, which creates a luminosity-unbiased sample over a wide redshift range. However, we want to know the total number of bursts which are occurring. The assumed GRB luminosity function is taken from Pescalli et al. (2016). This has a power-law slope of −1.2 below, and −1.92 above, a break energy of 5 × 10⁵⁰ erg. The SHOALS data give us the number of events per redshift bin above 10⁵¹ erg, to obtain the total number at each redshift for comparison with modelled rates we integrate over the function down to a lower limit E_low. Finally, models are selected as GRB candidates if they have a specific angular momentum at core collapse greater than j_cut, which is proportional to (Z/Z_⊙)ⁿ.

We perform a Markov Chain Monte Carlo (MCMC) analysis to infer the probability density functions of θ, E_low, j_cut, and n. The python package emcee is used to perform our MCMC sampling (Goodman & Weare 2010; Foreman-Mackey et al. 2013). The posterior distribution of the parameters, P(θ, L_E, j_cut, n|D, M) is given by the product of the log likelihood and log prior

$$\begin{eqnarray*} P({\theta },L_\mathrm{E},j_\mathrm{cut},n|D,M) &\propto & P(D|{\theta },E_\mathrm{low},j_\mathrm{cut},n,M) \nonumber\\ &&\times \, P({\theta },E_\mathrm{low},j_\mathrm{cut},n), \end{eqnarray*}$$

(6)

where D denotes the data (in this case the SHOALS rates and their quoted uncertainties), and M is the two-pathway model as previously described. The priors assumed, P(θ, E_low, j_cut, n), are as follows:

• For the half-opening angle, we simply limit θ to the range 0 < θ < 22.5, the least informative prior we can use whilst ruling out very weakly beamed (total opening angle >45 deg) or isotropic emission.

• We have limited prior knowledge on the lowest possible GRB luminosity. We therefore restrict E_low to the range 45 < log₁₀(E_low/erg) < 50.7, where the upper bound is the break in the assumed luminosity function, and the lower bound is arbitrarily low. Although events are seldom seen at log₁₀(E_iso) < 50.7, we are trying to apply a minimal prior and allow for a small number of outliers in E_iso (Ajello et al. 2019).

• We use a top-hat prior which covers the range 16 < log₁₀(j_cut/cm² s⁻¹) < 19.3. The lower bound corresponds to the theoretical minimum j required for GRB production by a black hole central engine (e.g. Woosley 1993; MacFadyen & Woosley 1999). The upper bound is simply the maximum value of j obtained from our tidal calculations.

• For the index n, we again use a flat prior covering all physically realistic values (0 < n < 15). GRBs are rarer at high metallicity, and high stellar envelope opacity could help to impede jet propagation. Because opacity scales approximately as Z (e.g. Vink et al. 2001), we might expect n = 1 to be most probable, however, this is not favoured in our prior.

Initializing the MCMC with 100 walkers, 2000 steps, and a burn-in of 50 steps (checked by visual inspection), we obtain the joint marginalized distributions and correlations shown in the corner plot of Fig. 8.

Figure 8.

Covariances between the four fitted parameters in our two-pathway GRB model, and the marginalized posterior probability density distributions for each. Vertical dashed lines mark the 16th and 84th percentiles, and the orange crosshairs indicate the distributions medians. 1D histograms at the top of each column represent the marginalized distribution of the parameter indicated on the x-axis for that column.

5 PREDICTED LONG GAMMA-RAY BURST RATES AND PROGENITORS

5.1 MCMC results

The marginalized probability distribution in the half-opening angle θ is similar in form to those found by Racusin et al. (2009), and the well-fitted sample of Ryan et al. (2015). The posterior distribution in E_low favours a critical isotropic equivalent energy of 10^48.1 erg. This result arises despite a minimal prior which allows for all physically reasonable values below the break in the assumed luminosity function. The cut-off is in good agreement with observations, which have yielded bursts with log₁₀(E_iso) < 48 on only a few occasions, despite such events being theoretically detectable at low redshift (e.g. Racusin et al. 2009; Ajello et al. 2019). The posterior distribution of the angular momentum cut-off log₁₀(j_cut,⊙) has its median at 18.74. Finally, the distribution of the power-law index n drops away at very high values, with a peak and median at ∼4. The mean, median, and percentiles of all four posterior distributions are listed in Table 2.

Table 1.

The bpass metallicity scaling used in this work, adapted from Xiao, Stanway & Eldridge (2018). For a discussion of metallicity scales and uncertainties, we refer the reader to Eldridge et al. (2017).

Metal fraction	Previous	New
by mass	12 + log(O/H)	12 + log(O/H)	12 + log(Fe/H)
0.0001	6.60	7.00	5.23
0.001	7.61	8.00	6.24
0.002	7.91	8.31	6.54
0.003	8.09	8.43	6.72
0.004	8.21	8.51	6.84
0.006	8.39	8.61	7.02
0.008	8.52	8.69	7.15
0.010	8.62	8.75	7.25
0.014	8.77	8.77	7.40
0.020	8.93	8.93	7.57
0.030	9.13	9.05	7.76
0.040	9.27	9.13	7.90

Table 2.

Properties of the parameter posterior distributions, obtained from the MCMC analysis carried out in Section 5.

Parameter	Mean	16th Percentile	Median	84th Percentile
θ/deg	10.62	3.12	9.88	18.75
log10(Elow/erg)	48.31	47.42	48.11	49.30
log10(jcut,⊙/cm2 s−1)	18.64	18.14	18.74	19.13
Znindex	3.87	1.39	3.81	5.94

In Fig. 9, we show the posterior probability density for our model, compared with the observed E > 10⁵¹ erg SHOALS rates. The model (as a function of j_cut and n) predicts an intrinsic rate, which is converted to an observed rate by the θ and E_low parameters as described above. The probability density shown in Fig. 9 therefore represents the posterior distribution for the correction to the intrinsic rate, which is a convolution of the corrections arising from the opening angle and lower luminosity limit.

Figure 9.

The distribution of observed GRB rates, as a function of redshift (or lookback time). Darker shading represents higher probability density. The observed SHOALS rates are shown with their uncertainties. bpass produces intrinsic rates, these have been corrected using the distributions of θ and E_low which are output as posteriors from the MCMC run.

For best estimates of the parameters values, we use the posterior medians and 68 per cent credible intervals given by the 16th and 84th percentiles. This gives |${\theta }=9.9^{+8.9}_{-6.8}$| deg, log₁₀(E_low/erg|$)=48.1^{+1.2}_{-0.7}$|⁠, log₁₀(j_cut,⊙/cm²  s⁻¹) = |$18.7^{+0.4}_{-0.6}$| and |$n=3.8^{+2.1}_{-2.4}$|⁠. Using the median values gives the fit shown in Fig. 10 (we do not fit to the highest redshift SHOALS point, which gives an event rate consistent with zero). We note that the fitting just the QHE pathway can produce similarly good results given different assumptions for θ and E_low, as demonstrated by Eldridge et al. (2019b), but we know that GRBs occur above |$0.2\, \mathrm{Z}_{\odot }$| metallicity, and therefore the QHE pathway cannot be the sole contributor.

Figure 10.

The bpass prediction for the intrinsic volumetric rate of GRBs arising from both QHE and tidal pathways, shown by the solid purple line. The dashed dark blue line below is the contribution from QHE progenitors, and the dashed orange line represents bursts from the tidal influence pathway. These rates were obtained by selecting black hole-producing, SE progenitors with a specific angular momentum cut that best reproduces the evolution of the rate over cosmic time. The SHOALS rates have then been corrected for beaming and a luminosity function, using the medians of the θ and E_low MCMC posterior distributions.

5.2 Metallicity distribution

An independent test is to compare the metallicity distribution predicted for our tidal and QHE GRB progenitors, to that of observed host galaxies. In Fig. 11, we show the synthetic metallicity distribution of our GRB-producing stars at two redshifts, z = 0.2 and z = 1.5. Also shown are host galaxy metallicity distributions from Japelj et al. (2018), Modjaz et al. (2019), Graham, Schady & Fruchter (2019), and Palmerio et al. (2019).

Figure 11.

Cumulative metallicity distributions of our GRB progenitors at z = 0.2 and z = 1.5, marked by shaded grey and green lines, compared with observed host galaxy distributions from Japelj et al. (2018), Modjaz et al. (2019), Palmerio et al. (2019), and Graham et al. (2019). The comparison data has been shifted so that at the mass fraction used to define Solar metallicity in their scale, 12 + log(O/H) is the same as the bpass value at that mass fraction. Marked on the plot are |$0.2\, \mathrm{Z}_{\odot }$| and |$\, \mathrm{Z}_{\odot }$| metallicities in the (modified) bpass scaling.

To draw comparisons between metallicity distributions, we need to ensure that the scales being used do not have significant offsets. Graham et al. (2019) used the metallicity diagnostic and scaling of Kobulnicky & Kewley (2004), with Solar metallicity defined at 12 + log(O/H) = 8.69 (Allende Prieto, Lambert & Asplund 2001). This corresponds in their scale to a metal mass fraction of 0.014. In our bpass scaling, a mass of fraction of 0.014 corresponds to 12 + log(O/H) = 8.76. To reconcile this with Kobulnicky & Kewley (2004), 0.07 dex is added to each value in the Graham et al. (2019) distribution.

The other three comparison data sets use a Maiolino et al. (2008) scaling (Modjaz et al. 2019 provide a variety, we choose the same scaling for consistency), where Solar is again at 12 + log(O/H) = 8.69, but this now corresponds to 0.0134 by mass fraction (Asplund et al. 2009). Again, these data sets are shifted by 0.07, which brings the 0.0134 Solar value into agreement with bpass at 12+log(O/H) = 8.76. For a discussion of these issues, and their impact within bpass, we refer the reader to Eldridge et al. (2017) and Xiao et al. (2018).

Interestingly, there appears to be little evolution of the metallicity distribution in GRBs predicted by the bpass models with redshift, with an overall shift to lower values of only ∼0.2 dex between redshifts 0 and 5. The data similarly show a lack of variation in the observed fraction of high metallicity bursts out to redshift 2.5 (Graham et al. 2019). The observed samples none the less span a wide redshift range, and we compare each to the closer of the two redshift curves shown on Fig. 11. Anderson–Darling tests between the bpass results and the data fail to reject the null hypothesis (i.e. p > 0.05) that they are drawn from the same distribution, in every case except for Graham et al. (2019).

5.3 Delay-time distribution

In Fig. 12, the delay-time distribution for the QHE and tidal progenitors are shown. In the tidal case this is for all metallicities considered, whereas the QHE distribution is limited to |$Z\lt 0.2\, \mathrm{Z}_{\odot }$| by construction. The new tidal pathways have shorter delay times, decreasing the mean temporal offset between star formation and GRB events. We note that these times are technically only until the end of core carbon burning. However, the final stages of core burning before core collapse occupy <<1 Myr (e.g. Groh et al. 2014). The GRB progenitors have among the shortest delay times of any stars, particularly so for the tidal GRBs. This implies that they will be also be among the most luminous main-sequence stars. Unlike the QHE pathway, in which it is lower mass secondary stars that produce GRBs, in the tidal pathway the progenitor will usually not have been subject to a supernova kick, or have received a smaller kick. The progenitor stars would therefore be preferentially formed in the brightest regions of their host galaxies, and stay there, leading to a GRB distribution within their hosts that traces the host light – a trend which has been previously been observed (e.g. Fruchter et al. 2006; Eldridge et al. 2011; Lyman et al. 2017).

Figure 12.

DTDs for QHE and tidal GRB progenitors, across all metallicities considered in this work. QHE pathways are limited to |$Z\lt 0.2\, \mathrm{Z}_{\odot }$| by construction. The contribution from each model is weighted according to the (tide adjusted) bpass weightings, and delays are measured from zero-age main sequence. The total area of each distribution is normalized to unity.

5.4 Progenitor systems

Table 3 shows distribution statistics for the tidal progenitor systems in mass ratio, orbital period, initial mass, final mass, and delay time parameter space, over the metallicity range Z = 0.008–0.020. Stars with final masses of |${\sim }20\, \mathrm{M}_{\odot }$|⁠, in tight binaries with large mass ratios, are the most frequent progenitors. A corner plot showing this information, and covariances between parameters, is available in Fig. A2 of the online Appendix. There are signs of an excess in the number of twin systems (i.e. systems with mass-ratio near unity) which give rise to GRBs, however this is simply reflecting the increased likelihood of twin systems in all massive binaries, and is not specific to GRB progenitors.

Table 3.

Properties of the tidal GRB progenitors over the metallicity range 0.008–0.020. For GRBs arising from the primary star (the majority), mass ratio q, initial mass, and log(P) are all given at ZAMS. For GRBs arising from the secondary star in a binary, values are given immediately after the supernova of its primary companion. We list the minimum, maximum, mean, and standard deviation σ of the parameter distributions for the population.

Property	Min	Max	Mean	σ
Mass ratio q	0.03	0.90	0.48	0.28
log10(P/d)	0.00	2.20	0.45	0.50
Initial mass/|$\, \mathrm{M}_{\odot }$|	15.0	300.0	84.7	67.6
Final mass/|$\, \mathrm{M}_{\odot }$|	8.5	46.6	18.7	8.6
Delay time/Myr	2.5	16.0	4.5	2.1

5.5 Core angular momentum

Finally, we turn our attention to the angular momentum distribution of the tidally spun models selected as GRB progenitors. The angular momentum cut applies to j = Ω × R², where R is the radius of the star, however the quantity of interest for GRBs is the interior specific angular momentum. Theoretical modelling (Woosley 1993; MacFadyen & Woosley 1999) suggests that to launch a jet, |${\gtrsim}16\, \mathrm{cm}^{2}\, \mathrm{s}^{-1}$| is required at the innermost stable orbit around the newly formed black hole. To estimate these j values in our tidal GRB models, we make two assumptions. First, that the star has a constant rotational angular velocity Ω throughout its structure at the end point of our models. Secondly, that initial radius of the material which will form the accretion disc of the nascent black hole is at R_boundary – the radius in the pre-collapse star which encloses the post-collapse remnant mass. This material is assumed to retain its specific angular momentum during core collapse.

To calculate R_boundary, we use the files output directly from the bpass version of the stars code, which includes information on the radial structure of the stellar models.¹ The boundary radius is calculated by summing shells of mass until the remnant mass is enclosed, from r = 0 to r = R_boundary. The specific angular momentum at this radius is then given by

$$\begin{eqnarray*} j_\mathrm{boundary} = \frac{J_\mathrm{shell}}{M_\mathrm{shell}} = \frac{2}{3} \Omega R_\mathrm{boundary}^{2} , \end{eqnarray*}$$

(7)

where the angular momentum and mass of the shell just outside the boundary radius are J_shell and M_shell.

In Fig. 13, we show the specific angular momenta, evaluated at the remnant-ejecta boundary r = R_boundary, for the tidal GRB models immediately before core collapse in three metallicity bins. The angular velocity assumed in all cases is j_cut/R², the minimum value allowed by our metallicity-dependent surface momentum cut. We show two cases, the first has a Z^3.8 dependence on this cut (the favoured value from our MCMC run), the other has n = 1.4 (this is the lower, 1σ equivalent bound). We might expect no strong metallicity dependence on the angular momentum that a new-forming black hole requires to launch a jet, although an index of ∼1 may be expected for jet escape if envelope opacity is solely responsible for suppressing jets and hence dominating any Z dependence (Vink et al. 2001; Vink & de Koter 2005).

Figure 13.

The minimum specific angular momenta required for a GRB, measured at the remnant-ejecta boundary of each model at the point of collapse. These correspond to collapsing stars with minimum surface specific angular momenta j_cut. The upper panel assumes a Z^3.8 metallicity dependence, and the lower assumes that j_cut ∝ Z^1.4. Detailed stellar interior models are used to calculate j_boundary as discussed in the text. The models, which include primary and secondary stars, are binned into three metallicity ranges: low (10 × 10⁻⁴ ≤ Z ≤ 0.004), moderate (0.006 ≤ Z ≤ 0.008), and high (0.010 ≤ Z ≤ 0.040). The typical boundary momenta required are in the range log₁₀(j_boundary/cm² s⁻¹) ∼ 13 − 17 for n = 3.8, and 14.5–16.5 for n = 1.4, although these are minimum values and could easily be ∼100 times greater (see discussion in Section 6.1).

The distributions shown in Fig. 13 represent lower limits in that the angular momentum of any specific GRB progenitor may exceed the fitted cut level for the population. The distributions show scatter around 10¹³–10¹⁷ cm s⁻¹ (a wide spread, n = 3.8) and 10^14.5–10^16.5 cm s⁻¹ (more peaked, n = 1.4). We note that using an n ∼ 1 metallicity dependence preferentially shifts the low Z models to greater specific angular momenta, and brings the three metallicity bins into good agreement.

6 DISCUSSION

6.1 The production of GRBs

We have identified the subset of probability weighted stellar evolution models which are likely to generate a GRB through either quasi-homogenous evolution or the results of tidal interactions modifying the angular momentum of the progenitor star. The location of the tidal GRB progenitors on the Hertzsprung–Russell diagram is shown in Fig. 14, for all models in the metallicity range Z = 0.008 to Z = 0.020. Purple stars indicate the predicted total optical light from the binary (i.e. both primary and secondary star, or secondary plus remnant in rare cases) immediately before the supernova explosion, where the more luminous component is assumed to dominate the temperature measurement. The yellow stars indicate the properties of the surviving binary companion expected to be observable after the GRB has faded. Grey circles represent the properties of individual progenitor stars immediately before core collapse. In some cases, this progenitor is the secondary in the original ZAMS binary.

Figure 14.

Every model (primary and secondary) which produces a GRB via the tidal pathway, in the metallicity range Z = 0.008 to Z = 0.020, shown on the HR diagram. The model values for the individual pre core-collapse stars that go GRB are marked by grey circles. Also marked are the secondary stars left behind after a primary goes GRB (orange stars) and the unresolved systems that would be seen in pre-primary explosion imaging (purple stars). The shading represents the number density of models. Two example evolutionary tracks are overlaid. Track (A) is for a 50 M_⊙ primary with a 45 M_⊙ companion, starting with an orbital period of 0.2 d at a metallicity of 0.5 Z_⊙. Track (B) follows a 150 M_⊙ secondary star with an 11 M_⊙ black hole companion, starting with an orbital period of 1.4 d, at 0.4 Z_⊙.

The GRB progenitor binary component is often not responsible for all the light coming from an observed progenitor system; this is particularly true for the twin systems identified above, in which the secondary is likely to be very nearly as bright and evolved as the primary. Immediately prior to core collapse, 5 per cent of progenitor systems between |$0.4\, \mathrm{Z}_{\odot }$| and |$\, \mathrm{Z}_{\odot }$| metallicity have a secondary star that is more luminous than the pre-explosion primary.

Progenitors end in the hot and bright region on the upper left, as luminous or more luminous than typical Wolf–Rayet stars seen in the Local Group (Neugent & Massey 2019). Fig. 14 and Table 3 indicate that main-sequence stars are the most frequent companions left behind after a primary star goes GRB in the models considered here. This is consistent with earlier findings. Zapartas et al. (2017b) predict the companions expected for SE (type IIb, Ib, and Ic) supernovae using the binary_c rapid population synthesis code. They found that at |$0.3\, \mathrm{Z}_{\odot }$| metallicity, given their assumed IMF and binary parameters, ∼68 per cent of the progenitors should have a main-sequence companion at the point of explosion, and the remainder have compact object companions.

The modelled intrinsic GRB rate shown in Fig. 10 is around ∼10 per cent of the type Ic SN rate expected in the local Universe (z < 1). This is consistent with previous estimates of this fraction (Fryer et al. 2007) from observational data, corrected for selection effects.

The key property in our model determining whether or not a massive SE star produces a GRB is its internal specific angular momentum. Our estimates (using n = 3.8) for the minimum boundary specific angular momentum required to launch a GRB, log₁₀(j_boundary/cm² s⁻¹) ∼ 13–17, drops well below that expected from theory. This is also lower than the values found from detailed modelling of massive star interiors (Heger, Langer & Woosley 2000; Yoon, Dierks & Langer 2012; Fryer et al. 2019). The bpass models rapidly evolve through their final stages, and do not track the end phases of evolution in detail. The core properties we have used are therefore more representative of stellar structure at the end of core carbon burning, rather than at collapse, and stellar cores can contract further by around two orders of magnitude in density after carbon burning (Eldridge et al. 2019a), potentially raising the specific angular momentum of a given shell. During this time, the assumption of solid body rotation almost certainly breaks down, but circumstances in which the envelope rotates faster than the core are highly unlikely. The j values shown in Fig. 13 are already lower limits since we have adopted our log₁₀(j_cut) minimum value. Actual model values extend up to log₁₀(j) ∼ 19. If a further core contraction occurs, as expected, the true values may be higher still by several orders of magnitude (e.g. Heger et al. 2000). Even a modest cumulative increase of 2 dex would push the distribution towards log₁₀(j_boundary/cm² s⁻¹) ∼15–19, in better agreement with theoretical predictions (Woosley 1993; Woosley & MacFadyen 1999). This is the first time this critical threshold has been derived from observational data, albeit through fitting with stellar models. It arises naturally, without fine tuning, from our two-pathway model and associated assumptions.

6.2 The metallicity dependence of GRBs

GRBs show a deficit with respect to the cosmic volume-averaged star formation rate SFRD at low redshift, but occur in proportion to it at z ∼ 3 and above (Eldridge et al. 2019a). This is reflecting the well documented high-metallicity aversion of GRBs (e.g. Fruchter et al. 2006; Greiner et al. 2015; Modjaz et al. 2019), which is demonstrated further in Fig. 15. Here we show the GRB progenitor production efficiency, defined as GRB progenitor fraction per unit star formation at each metallicity, given our assumed IMF and binary parameters. There is a clear aversion to high metallicity. The QHE pathway becomes increasingly scarce as metallicity increases, until the Z = 0.004 mass fraction cut-off implemented by bpass is reached. Tidal pathways, on the other hand, peak at around 20 per cent of Solar. This is likely due to a trade-off between tides being more effective when the stellar envelope is extended, and increasing angular momentum loss through winds.

Figure 15.

The efficiency of GRB progenitor production across metallicity. For a given mass of gas at a single metallicity, we would expect N_GRB progenitors to be formed per Solar mass of gas. We show this number for all GRB progenitors, and separately for the QHE pathway. QHE progenitors show a smooth decline in occurrence rate as metallicity increases, before the 0.2 Z_⊙ cut-off is reached. The tidally induced pathways, however, peak around this metallicity.

We can further test this by examining the metallicity distribution of predicted progenitors. In Fig. 11, we compared our synthetic GRB metallicity distributions to those from observed host galaxies, making the assumption that the metallicity of the host stellar population changes by a negligible amount over the lifetime of a GRB progenitor. The resultant distribution accounts for the cosmic star formation history and model weightings, and is broadly consistent with the observed host galaxy population.

Given that our predicted metallicity distribution and volumetric rate evolution is in good agreement with observations, using our best-fitting Bayesian analysis values, thought must be given to the strong metallicity dependence suggested by that analysis, Z ∝ Zⁿ where |$n=3.8^{+2.1}_{-2.4}$|⁠. The uncertainties on this parameter are large, and at their upper end suggest a far stronger metallicity dependence for GRBs than seen in earlier work (e.g. Trenti et al. 2015).

We have chosen in this analysis to place the metallicity dependence as a modifier of the angular momentum threshold. However, we expect that the minimum angular momentum threshold at the remnant-ejecta boundary (i.e. at the point at which a jet is launched) is independent of metallicity. Thus the effective threshold encodes information not on jet launching but rather on the metallicity dependence of whether the jet, once launched, can escape the stellar envelope or is stifled before producing a GRB.

Successful jet break-out relies on a number of different properties of the collapsing star. In principle, it will depend on the column density of the material through which the jet must tunnel (i.e. the density and thickness of the stellar envelope at point of collapse), the probability of photons interacting with that material (i.e. the envelope opacity), and also the time-scale for jet escape (i.e. if the central engine deactivates before breakout is achieved there will be no visible event).

We calculate the envelope column density for each GRB progenitor model by summing the density of mass shells from the core–envelope boundary out to the surface, multiplying by the shell thickness at each stage. The detailed stars outputs are again used. We find that although high metallicity stars do have some of the highest columns, there is no distinct trend across the full metallicity range. We therefore cannot confidently attribute the high Z dependence suggested by our model fitting analysis to a column density effect.

The opacity of the ejecta which must be tunnelled through is also dependent on the metal content of the stellar material. For each photon, its probability of interaction scales with the number of possible electron energy level transitions which it may be able to excite. Heavier metallic elements, with their extended electron shells, dominate this probability, and so the opacity scales broadly linearly with the abundance of iron group elements. For higher opacities, more energy is dissipated from the jet in exciting electrons and so an initially more relativistic, more collimated jet is required. This corresponds to a greater reservoir of angular momentum also required to successfully tunnel through the envelope. Opacity is roughly proportional to metal mass fraction Z (Vink et al. 2001; Vink & de Koter 2005), and would lead to n ∼ 1 in our formalism. It therefore cannot reproduce the n = 3.8 dependence that results as a best fit from our Bayesian inference.

If instead we assume the lower 1σ bound from our analysis (n = 1.4), the resultant j distribution in Fig. 13 is good agreement with theory and removes any clear metallicity dependence at the point where the jet launches. A lower n index also introduces bursts at |$\, \mathrm{Z}_{\odot }$| and above, improving the metallicity distribution agreement with observations. The question then is why Z^3.8 was favoured by the Bayesian inference. In Fig. 10, a model which tracks the assumed cosmic star formation history overestimates the GRB rate at low redshift and underestimates it at high redshift. Larger powers of Z rectify this, allowing more bursts to occur at low metallicity (fewer at high values) – boosting the high-redshift rate (diminishing it at low redshift). However, a Zⁿ dependence that is too strong begins to remove too many bursts and makes it difficult to reproduce the observed population numbers with plausible E_iso and θ corrections. The |$n=3.8^{+2.1}_{-2.4}$| best fit arises from the combination of these factors, and inherits all their uncertainties. In particular we have not considered the substantial uncertainty on the cosmic volume-averaged metallicity evolution and its scatter as a function of redshift, or uncertainties in the somewhat better constrained cosmic SFRD history. The inclusion of other progenitor pathways, which could have a different occurrence rate over cosmic history, may also affect the best-fitting index n. Thus while we clearly favour a non-zero metallicity dependence, we caution against over interpretation of the best-fitting value, and note that an n ∼ 1 dependence is both permitted and explained by a physically plausible mechanism.

We note that the derived metallicity dependence arises from a fit to the evolution in the inferred GRB production rate. The fraction of very low metallicity (i.e Z < 0.001) models contributing to the volume averaged rate is very low and so the derived Zⁿ dependence is likely poorly constrained or not applicable at these very low metallicities. Indeed, it is possible that below a minimum metallicity threshold non-iron-dominated opacities begin to dominate and the metallicity dependence will plateau. There is very little difference in the stellar atmosphere opacity tables used by bpass in this regime (Eldridge et al. 2017).

6.3 Uncertainties in the stellar modelling

Throughout this paper, we have made assumptions and simplifications with regards to the tidal evolution of binary systems, and the structure of high-mass stars. For each model an initial weighting was selected from the IMF and binary parameter distributions, and assigned a value for an initial spin. The default bpass IMF incorporates stars up to |$300\, \mathrm{M}_{\odot }$| at ZAMS, in order to accommodate the rare but important Very Massive Star population seen in the Local Group. However, the occurrence of such stars and the distribution of binary parameters may well be metallicity dependent. This is unconstrained by observations and so not accounted for in the current version of bpass (v2.2.1, Stanway & Eldridge 2018, 2019).

The initial stellar angular velocity chosen for every primary model was 40 per cent of the Keplarian value at the equator (a typical value for OB stars Dufton et al. 2019). While a distribution of velocities could be sampled, the exact choice of Ω does not have a significant impact on tidal evolution (also seen by Zapartas et al. 2017a), as there is typically much more angular momentum stored in the orbit than in stellar spins. We also assume that the binaries start in circular orbits with their spin vectors aligned, so that only |$\dot{a}$| and |$\dot{{\Omega }}$| need be considered.

Given that most significant changes due to tides occur with the onset of mass transfer (Hurley et al. 2002), changes to evolution before that stage of factor a few are unlikely to have a major impact on our results. This is relevant to our assumptions concerning stellar structure. We take the value of the apsidal motion constant, k (which dependes on internal structure) to be 0.05 for every model. The theoretical modelled range of values is 0.01–0.1 (Zahn 1975). ΔΩ and Δa both scale with k⁻¹. Varying k by a factor of a few will have similarly small effects on these parameters. For the system shown in Fig. 4, varying k from 0.01 to 0.1 makes no difference to the size or timing of the model jump made by our tidal algorithm, and the system ends in the same state independent of k.

We note that the secondary star in our tidal prescription is treated as a point mass. For evolved binaries, where this object is a compact remnant, the approximation is reasonable. Otherwise, tidal distortion of the secondary may play a role in the system evolution, even if it remains inside its Roche lobe. One possible outcome of spin–orbit interactions is a merger. For this to occur, the total system angular momentum – the sum of the orbital and spin components – must be less than a critical value. Angular momentum is then transferred from the orbit to the spin, but the stars merge before a synchronization can occur (Eldridge & Tout 2019). Given our starting velocity of 0.4Ω_crit, mergers solely due to tides are rare in our model set, only occurring in very tight binaries.

Finally, a major assumption in this analysis is that most black hole producing core-collapse events do not produce a visible supernova. The picture is more complex than there being a simple mass cut-off, as Dessart, O’Connor & Ott (2012b) discuss in their work on black hole production in GRBs, with other factors including the rotationally driven magnetic field and envelope structure. Core compactness is also important criteria for successful supernovae (Sukhbold et al. 2016; Ertl et al. 2019), leading to so-called ‘islands of explodability’. We address the core compactness by assuming that the radius of interest within the pre-collapse star is that which encloses the eventual remnant mass. This is reasonable approximation, provided that there is negligible fallback accretion. We also assume that the core and envelope co-rotate, with a flat angular velocity profile throughout the entire star. As previously discussed, any deviation from this will likely result in the core spinning faster than assumed. Such an increase in core spin would increase the number of rapidly spinning cores, and therefore increase the inferred j_cut,⊙ threshold, improving agreement with collapsar theory.

6.4 Magnetars as GRB central engines

Although we have assumed in this paper that GRBs operate under the collapsar mechanism, which requires black hole formation upon core collapse, it is possible that newly formed magnetars might be able to launch jets too. Indeed, there is some evidence from GRB-SN energetics that this is the case (Mazzali et al. 2014). Therefore, we relaxed our remnant mass constraint to 1.4 M_⊙, allowing for the possibility that neutron stars are formed in core-collapse GRBs, and re-performed the MCMC analysis as described in Section 4.2. We can see from Fig. 5 that while black holes form readily at low Z, type Ic SNe with neutron star remnants are rare. By including neutron star producing Ic events, we are preferentially adding models at high metallicity. Because these end their lives spinning slower, most of those added are then rejected by the angular momentum cut. Therefore, neutron star-forming events can contribute, but are not an important pathway if black hole forming GRBs are also considered. The small difference their inclusion makes is demonstrated in Fig. A3 of the online Appendix.

We then took the extreme case that GRBs can only occur if a neutron star is formed in core collapse. An MCMC run under this assumption, as described in Section 4.2, yielded |${\theta }=12.8^{+6.4}_{-4.6}$| deg, log₁₀(E_low/erg|$)=49.6^{+0.6}_{-0.5}$|⁠, log₁₀(j_cut,⊙/cm² s⁻¹) = |$18.8^{+0.3}_{-0.4}$| and |$n=2.2^{+1.0}_{-1.3}$|⁠. Because SE, fast-spinning, neutron star-forming progenitors are rarer (than the equivalent black hole forming events), the best-fitting parameter values are naturally quite different, as shown in Fig. A3 of the Appendix. The most notable among them is E_low – the minimum isotropic equivalent energy is nearly two orders of magnitude higher than the faintest GRBs observed. Therefore, barring strong co-variances between the four model parameters which have not been considered, or a very different luminosity function than that assumed, this appears to disfavour magnetars as the sole central engines capable of launching GRBs. Furthermore, the metallicity distribution assuming just magnetar engines is inconsistent with the host galaxy data (see Fig. A4 of the online Appendix). There is narrow range of metallicity in which stars lose enough mass to produce neutron star remnants, but retain enough to not fully spin-down. Anderson–Darling p-values between the synthetic distribution and the hosts data sets are all <0.05, rejecting the null hypothesis that they are drawn from the same distribution at the 2σ level.

6.5 Future possibilities

In Fig. 14, we predicted that the progenitors of GRBs induced through tidal evolution are amongst the most luminous SE stars. They are also extremely hot. Thus despite their high luminosities they may be challenging to observe, requiring ultraviolet, or very blue-sensitive optical, detectors. Supernova 2017ein, the only SN Ic to have a progenitor system identified so far, occurred at a redshift ∼0.0027, with the candidate identified in archival Hubble Space Telescope (HST) imaging (Kilpatrick et al. 2018; Van Dyk et al. 2018). The star, if single, had an absolute F555W magnitude of −7.5. The nearest GRB confirmed so far is GRB 980425, associated with SN 1998bw. This occurred at z ∼ 0.0085 (Galama et al. 1998), corresponding to a luminosity distance which is approximately three times greater than SN 2017ein. If SN 2017ein had occurred at the distance of SN 1998bw, the progenitor (or unresolved progenitor system) would have had an apparent magnitude of ∼25.5 in this band.

It is also interesting to consider the possibility of surviving companions being seen once the supernova and GRB afterglow have faded. We have identified these companions as most-likely main-sequence stars, with a slightly preference for large mass ratios or twin systems. Such stars are faint, but not beyond the realms of possibility for HST or upcoming facilities such as the European Extremely Large Telescope (E-ELT) or the James Webb Space Telescope (JWST). The low metallicity QHE pathway is restricted to high redshifts where these kind of observations are unfeasible. However, the tidal GRB pathways we have identified, where tides maintain spin in an SE progenitor, are a plausible target.

Although we can reproduce the observed GRB rate using mass-transfer QHE and tidally spun-up pathways, it is worth considering other, rarer possibilities. A number of such pathways have been identified in the literature, including the merger of helium stars with a companion (Fryer & Heger 2005; Detmers et al. 2008; de Mink et al. 2013, 2014), binary-driven hypernovae (Rueda, Ruffini & Wang 2019), and white dwarf-compact object mergers (Fryer et al. 1999; Caito et al. 2009).

Another scenario we have not considered is that tidally spun-up stars could enter the QHE regime (de Mink et al. 2013; de Mink & Mandel 2016; Mandel & de Mink 2016; Song et al. 2016; Song et al. 2016; Marchant et al. 2016, 2017) – in our formalism the only way to induce QHE is through mass transfer at very low metallicity. Tidal interactions may in some cases lead to increased mass transfer and QHE, and there is an overall slight increase in the frequency of QHE occurring when tides are considered (see Fig. A1 of the online Appendix). This pathway has received significant attention in the literature as a route to producing massive black hole binaries. However, because the low metallicity, low mass-loss GRB requirement holds even if QHE is tidally induced, such events will not significantly contribute to the overall observed GRB rate. We note that our analysis above has produced a good match to the observed volumetric rates and properties of known GRBs without assuming exotic or extremely rare progenitor pathways. We therefore leave a full bpass study of the viability and occurrence rates of rarer alternative pathways for future consideration.

Finally, an extension of the bpass modelling project is the synthesis of SN light curves from the stellar structure models at the point of explosion (CURVEPOPS; Eldridge et al. 2018, 2019a). In a similar fashion, the light curves arising from the type Ic supernovae of our candidate GRB progenitors could be synthesized, and a comparison between the CURVEPOPS light curves and GRB-SN light curves made. This will be explored in future work.

7 CONCLUSIONS

We have used the binary stellar evolution models of bpass to investigate the nature of supernova and core-collapse GRB progenitors. By considering tidal interactions, and applying a prescription for the cosmic star formation rate history, we have shown that two binary pathways can explain the bulk of the observed GRB population. The first involves secondary stars spun up by accretion into a quasi-homogeneous state, and the second occurs when massive stars in binaries have their envelope removed while tides maintain the required angular momentum for jet production. This two pathway model can reproduce the rates as a function of redshift and the observed host metallicity spread, using reasonable parameter distributions for the jet opening angle and isotropic equivalent energy. The model predicts minimum angular momentum requirements which are in general agreement with collapsar GRB theory. Finally, predictions are made for the nature of the progenitor systems and their companions, which may be testable with next generation facilities. The exploratory work presented here does not represent a definitive or unique solution to the GRB progenitor problem, and there are several areas which need to be further explored or developed in future work. None the less, our approach can simultaneously reproduce several aspects of long GRB theory and observation, demonstrating the utility of population synthesis models in transient progenitor studies.

SUPPORTING INFORMATION

Appendix A.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

AAC is supported by Science and Technology Facilities Council (STFC) grant 1763016. AAC also thanks the William Edwards educational charity. ERS is supported by STFC consolidated grant ST/P000495/1. JJE acknowledges support from the University of Auckland and also the Royal Society Te Apārangi of New Zealand under Marsden Fund grant UOA1818.

We acknowledge use of Python algorithms including the corner module (Foreman-Mackey 2016), the excellent Python-based introduction to Bayesianism by VanderPlas (2014), and the use of Ned Wright and James Schombert’s Python cosmology calculator (Wright 2006).

Footnotes

1

Note: these detailed output files do not form part of the standard bpass stellar model data release due to their data volume and technical complexity.