Properties of High-redshift Gamma-Ray Bursts

The star formation history, in particular at high redshift, is difficult to determine. Studies have been done using a wide range of probes (emission from star formation regions, supernovae, and GRBs) over a broad range of photon energies from the infrared to gamma-rays. In general, observations of the emission from star-forming regions (Madau et al. 1996; Hopkins & Beacom 2006; Bouwens et al. 2020) predict a much sharper drop-off at high redshifts than astrophysical transients (Gal-Yam & Maoz 2004; Wiggins et al. 2015). From dust obscuration to the evolution of the IMF, uncertainties in model assumptions make it difficult to ascertain a firm star formation history at high redshift.

For our models, we adopt the star formation rate from Madau & Dickinson (2014) as the benchmark star formation rate in our simulation. This star formation rate is described by the following equation,

$$\begin{eqnarray}&&\rho =0.01\displaystyle \frac{{(1+z)}^{2.7}}{1+{[(1+z)/2.9]}^{5.6}},\end{eqnarray} \tag{ 1 }$$

where ρ is in units of M_⊙ yr⁻¹ Mpc⁻³.

In addition, we also run the simulations of the main theoretical model using two other star formation rates from Madau & Fragos (2017) and Hopkins & Beacom (2006) in order to estimate how the predicted GRB rate might change due to different assumptions of the star formation rate. The star formation rate from Madau & Fragos (2017) is described as

$$\begin{eqnarray}&&\rho =0.01\displaystyle \frac{{(1+z)}^{2.6}}{1+{[(1+z)/3.2]}^{6.2}},\end{eqnarray} \tag{ 2 }$$

and the star formation rate from Hopkins & Beacom (2006) is

$$\begin{eqnarray}&&\mathrm{log}(\rho )=\left\{\begin{array}{l}3.28\mathrm{log}(1+z)-1.82,\ \ \ \ z\leqslant 1.04\\ -0.26\mathrm{log}(1+z)-0.724,\ 1.04\lt z\leqslant 4.48\\ -{f}_{\mathrm{highz}}\mathrm{log}(1+z)+4.99,\ \ \ z\gt 4.48\end{array}\right.\end{eqnarray} \tag{ 3 }$$

where f_highz dictates the drop in star formation at high redshift. Most observations of star formation predict that f_highz lies between 5 and 8 whereas astrophysical transients predict f_highz values in the 0–4 range. These differences could be explained by observational biases in the different samples. They may also be explained by an evolution in the progenitors of these transients as a function of redshift. In the next two subsections, we study the effects of metallicity (and hence redshift) on GRB and supernova progenitors. Although these models produce very similar results at low redshift (with just different shapes to their fits to the data), the different dependencies on the redshift at high values of the redshift can lead to very different rate predictions.

Long GRBs are produced in the most massive stars, and as such, the rate of GRBs depends sensitively on the poorly understood initial mass distributions of massive stars, a.k.a. the IMF. A number of stellar population studies in the Local Group of galaxies have been used to infer observed IMFs (Salpeter 1955; Miller & Scalo 1979; Kroupa 2001; Chabrier 2003). These studies typically focus on stars formed at high metallicities, but the general consensus is that, as long as there are some metals to drive cooling, the IMF will not change dramatically with metallicity (however, see Chon et al. 2021, discussed below). All of these studies obtain similar results for the mass distribution of stars above 1 M_⊙:

$$\begin{eqnarray}&&f(M)\propto {M}^{-\alpha }{dM},\end{eqnarray} \tag{ 4 }$$

where f(M) is the fraction of stars produced within the range of M and M + dM and α is typically set to ∼2.3 (e.g., for Salpeter-predicted α = 2.35). For our study, we use the Kroupa distribution to include the flattening of the IMF at low masses:

$$\begin{eqnarray}{\alpha }_{\mathrm{IMF}}=\left\{\begin{array}{ll}0.3, & \mathrm{if}\ M\lt 0.08\\ 1.3, & \mathrm{if}\ 0.08\lt M\lt 0.5\\ 2.35, & \mathrm{if}\ M\gt 0.5\end{array}\right.,\end{eqnarray} \tag{ 5 }$$

where stellar masses are measured in solar mass units and the IMF is normalized to unity over the stellar mass range: 0.0001–500 M_⊙. The mean stellar mass is 〈M〉 ∼ 0.7 and the fraction of single stars that create a neutron star or black hole is F_cr = 0.007 if one assumes that all stars with a mass greater than 8.5 M_⊙ experience core collapse. This corresponds to a star formation rate of 1.5 stars per solar mass of star formation and, in the Milky Way, where the star formation rate is roughly 1.9 M_⊙ yr⁻¹ (Kennicutt & Evans 2012), a supernova rate of roughly 1.4 supernovae per century. We further assume that, barring considerable mass loss from winds, stars above 25 M_⊙ collapse to form black holes and, if rotating with sufficient angular momentum, will form GRBs.

Metallicity can alter the nature of star formation in a number of ways. For example, as the metallicity decreases, metals cannot efficiently cool protostellar clumps, and the critical mass needed for clump collapse that initiates star formation increases (for a review, see Bromm 2013; Rosen & Krumholz 2020). In general, it is believed that, until the metallicity is very low, metal-driven cooling is sufficiently efficient to minimize the effect on the IMF. However, more recent studies argue that there is significant flattening of the IMF at a metallicity between 0.01 and 0.1z_⊙ (Chon et al. 2021), corresponding to a redshift of 2.5–5 (see Section 2.3.3).

In addition, metals can alter the fragmentation of stars where less fragmentation occurs (fewer binaries are produced) at lower metallicities (Bromm et al. 2001). The former effect serves to produce more massive stars, flattening the IMF and producing more GRB progenitors. Less fragmentation also flattens the IMF, but it produces fewer binaries. If binary interactions are essential to the formation of GRBs, the lack of fragmentation at high redshift could reduce the overall GRB rate at low metallicities (and high redshifts). Our understanding of star formation is far from complete. Magnetic fields can play an important role in star formation, but the effect of magnetic fields on star formation and, in particular, with redshift evolution, remains an active area of research (Rosen & Krumholz 2020). At high metallicities, stellar winds are much stronger and can also alter star formation (Calderón et al. 2020).

For our models, we will focus on two effects: the flattening of the IMF at low metallicity (from the larger protostellar collapse criterion and less fragmentation) and a smaller binary fraction (from the reduced fragmentation). For the IMF evolution, we must assume both the dependence of this flattening on the metallicity and evolution of the metallicity with redshift (which depends on the mixing of metals).

Because we do not exactly know the evolution of the IMF with metallicity, we combine both the star formation and metallicity evolution into a single set of uncertain quantities (we do need a standalone metallicity/redshift prescription for our stellar evolution models; see Section 2). To vary a wide range of evolutionary parameters, we use a simple parameterized prescription for the IMF slope (α_IMF) for massive stars with masses above 0.5 M_⊙ above a critical redshift z₀:

$$\begin{eqnarray}&&{\alpha }_{\mathrm{IMF}}=2.35-{f}_{\mathrm{IMF}}(z-{z}_{0}),\end{eqnarray} \tag{ 6 }$$

where we vary f_IMF from 0 to 0.4 and vary z₀ from 2.5 to 9. A z₀ = 2.5–5 assumes a flattening of the IMF much earlier than expected in simulations but in line with more recent results (Chon et al. 2021). As we shall see in our comparison to GRB observations, either this is required or a drastic change in the currently predicted star history at high redshift (or both) is necessary.

Many of the progenitor scenarios for long-duration GRBs require binary interactions (for reviews, see Fryer et al. 1999, 2007, 2019). If less fragmentation occurs at low metallicities, one might expect fewer binaries and, hence, fewer GRBs (Belczynski et al. 2010). In addition, with less expansion in the giant phase for massive stars, even if binaries form, we expect fewer binary interactions. We implement an additional decrement to the GRB rate with redshift:

$$\begin{eqnarray}&&{\mathrm{Fraction}}_{\mathrm{Binaries}}=1.0-{f}_{\mathrm{bin}}(z-{z}_{0}),\end{eqnarray} \tag{ 7 }$$

where we vary f_bin from 0 to 0.1. This interacting binary fraction will be multiplied to our long-duration GRB rate. Figure 1 shows the ratio of long GRBs to the supernova rate for a range of values of f_IMF, f_bin, and z₀. Here we assume that, without mass loss, stars above 25 M_⊙ form black holes, and 10% of these stars have sufficient angular momentum to form a GRB. Because only these massive stars form GRBs, the GRB rate depends sensitively on the IMF. At high metallicity, winds reduce core masses, effectively raising this critical GRB-forming mass beyond 25 M_⊙. Our f_bin = 0.1 model assumes that binaries are critical in forming GRBs and that the binary fraction goes down with the decrease in metallicity. At this time, we do not know the source of fast-rotating stellar cores, and this example just shows one possible effect of binary evolution.

Zoom In Zoom Out Reset image size

2.3.1. Metallicity and Wind Mass Loss

The primary effect of metallicity on stellar evolution is its role in mass loss. Standard mass loss in massive stars is produced through line-driven mechanisms (Nugis & Lamers 2000; Kudritzki 2002). At solar metallicities, this mass loss can be quite extreme, drastically altering the mass and core-structure of massive stars. These mass-loss–metallicity effects have long been incorporated into studies of massive-star progenitors of supernovae and GRBs (Heger et al. 2003; Fryer et al. 2007; Young & Fryer 2007). These studies found that for the BHAD engine for GRBs, mass loss prevents the formation of GRBs at high metallicity, but once the metallicity drops below roughly 0.1 solar metallicity, mass-loss and metallicity effects are, relative to the extreme effects at high metallicity, minor. There is a maximum mass at which the Type I collapsar forms. Without mass loss, above about 130 M_⊙, pair instabilities cause a runaway implosion and then an explosion that disrupts the star. However, if mass loss is extensive, this limit can be set by mass loss. As the metallicity in the star increases, stars lose too much mass to collapse to form black holes (Heger et al. 2003). Mass loss is most effective for the most massive stars and, with increasing metallicity, mass loss will first affect the most massive stars, altering their final outcome from a collapsar to a normal supernova. With increasing metallicity, the mass of stars where mass loss alters the final outcome begins to drop. To incorporate this effect, we assume the maximum mass of GRB formation based on the mass loss:

$$\begin{eqnarray}&&{M}_{\max }^{\mathrm{GRB}}=\min (130\,{M}_{\odot },-{f}_{\mathrm{wind}}{\mathrm{log}}_{10}({Z}_{\mathrm{metal}}/{z}_{\odot }){M}_{\odot }),\end{eqnarray} \tag{ 8 }$$

where Z_metal is the metallicity, z_⊙ is solar metallicity, and f_wind dictates the dependence of the mass loss on the metallicity. We vary f_wind from 30 to 9000, where 9000 limits mass loss to stars with metallicities at or above solar. For lower values, the effects of mass on the fate of massive stars can be more extended. We use 130 M_⊙ for the absolute maximum mass of GRB formation assuming that stars undergo pair instability above this limit. It is worth noting that mass loss can allow these massive stars to avoid this pair-production instability. For the IMF at high metallicities, these stars do not contribute significantly to the GRB rate and our removal of these will not affect our rate. Beyond the pair-instability mass gap, the most massive stars (ZAMS star mass above ∼300M_⊙) can form black holes. We will discuss these Population III collapsars below (Section 3.1).

Mass loss from winds also removes angular momentum from stars. For the fastest-rotating stars, winds are efficient at carrying away angular momentum, and stars with high mass loss will spin slower than those without this mass loss. Figure 2 shows the angular momentum profiles of two massive stars (40 and 65 M_⊙) at two metallicities (solar and 0.02 solar) from the GENEC code. These calculations assume no strong coupling between layers from a magnetic dynamo, producing the high-spin progenitors needed to produce BHADs. Even without this coupling, mass loss removes enough angular momentum to make much slower progenitors at high metallicities. The angular momentum at lower metallicities is a factor of 4 higher than the models at solar metallicity. If these models are correct, the GRB rate could increase as the metallicity drops from solar to 0.02 solar. However, we note that there is growing evidence that there is some coupling between burning layers, and this mass-loss effect on angular momentum could be muted (Belczynski et al. 2020).

Zoom In Zoom Out Reset image size

2.3.2. Other Metallicity Effects on Stellar Evolution

While small compared to the effects at high metallicity, metallicity effects below 0.1 z_solar can lead to detectable differences in the properties of a GRB. This is because the structure of the stellar cores continues to change with metallicity, producing more compact cores that, as we shall discuss in Section 3, alter both the power in the GRB jet and the nature of the supernova produced along with this GRB. For our study, we use the evolution of the progenitors from these GENEC models (Belczynski et al. 2020), using both the structure and angular momentum of these models as a function of metallicity. Figure 3 shows the radii and masses of the helium and carbon/oxygen cores as a function of ZAMS mass for a set of KEPLER models at solar and 10⁻⁴ solar metallicities (Woosley et al. 2002). Note that although the core masses are roughly the same (or less massive at solar metallicity), the radii of the cores are typically higher at solar metallicity.

Zoom In Zoom Out Reset image size

For our calculations, we have developed rough fits (as a function of mass and metallicity for these core masses and radii). For the CO core mass, we use

$$\begin{eqnarray}{M}_{\mathrm{CO}}=\left\{\begin{array}{ll}2.0+0.375({M}_{\mathrm{ZAMS}}-10), & \mathrm{if}\ {M}_{\mathrm{ZAMS}}\lt {M}_{\mathrm{crit}}\\ \max (2.0+0.375({M}_{\mathrm{ZAMS}}-10\,{M}_{\odot }) & \\ -0.1({M}_{\mathrm{ZAMS}}-{M}_{\mathrm{crit}}),2.0), & \mathrm{otherwise}\ \end{array}\right.,\end{eqnarray} \tag{ 9 }$$

where M_ZAMS is the ZAMS mass of the star in M_⊙ and ${M}_{\mathrm{crit}}=180/(5+2{\mathrm{log}}_{10}{Z}_{\mathrm{metal}})$. For the corresponding radius of these CO cores, we use

$$\begin{eqnarray}{\mathrm{log}}_{10}(r)=\left\{\begin{array}{ll}9.5+{M}_{\mathrm{ZAMS}}/100, & \mathrm{if}\ {M}_{\mathrm{ZAMS}}\lt {M}_{\mathrm{crit}}\\ \max (11,9.5+0.01{M}_{\mathrm{crit}} & \\ +0.35({M}_{\mathrm{ZAMS}}-{M}_{\mathrm{crit}})), & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 10 }$$

2.3.3. Metallicity Distributions with Redshift

For our metallicity–redshift relation, we use the same prescriptions used by Young & Fryer (2007) and Dominik et al. (2013). For the mean metallicity with redshift, we use the value described by Pei et al. (1999):

$$\begin{eqnarray}{Z}_{\mathrm{metal}}=\left\{\begin{array}{ll}{10}^{-{a}_{2}z}, & \mathrm{if}\ z\lt 3.2\\ {10}^{-{b}_{1}-{b}_{2}z}, & \mathrm{if}\ 3.2\lt z\lt 5\\ {10}^{-{c}_{1}-{c}_{2}z}, & \mathrm{if}\ z\gt 5\end{array}\right.,\end{eqnarray} \tag{ 11 }$$

where we use the fitting parameters set by Pei et al. (1999): a₂ = 0.5, b₁ = 0.8, b₂ = 0.25, c₁ = 0.2, c₂ = 0.4.

We allow for a distribution of metallicities as a function of redshift by assuming different galaxy masses have different metallicities (Tremonti et al. 2004):

$$\begin{eqnarray}&&12+\mathrm{log}(O/H)=-1.492+1.847\mathrm{log}{M}_{\mathrm{galaxy}}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&-0.08026{(\mathrm{log}{M}_{\mathrm{galaxy}})}^{2},\end{eqnarray} \tag{ 13 }$$

where the metallicity is set by the relative abundance of oxygen to hydrogen. For the galaxy mass distribution with redshift, we again use the formalism by Young & Fryer (2007) and Dominik et al. (2013) using the fitting formulae from Fontana et al. (2006):

$$\begin{eqnarray}&&f({M}_{\mathrm{galaxy}})\propto \mathrm{log}{\left({10}^{\mathrm{log}(M)-{M}_{* }(z)}\right)}^{{\alpha }_{* }(z)}\exp ({10}^{\mathrm{log}(M)-{M}_{* }(z)}),\end{eqnarray} \tag{ 14 }$$

where log(M) is the logarithm of the galaxy mass in units of solar masses, M_*(z) = 11.16 + 0.17z − 0.07z² and α_*(z) =−0.18 – 0.082z. Above a redshift of 4, we set the mass function to the value at a redshift of 4.

By assuming that GRB observables trace the single-star evolution, we are effectively studying collapsar progenitors for different collapsar types (Heger et al. 2003):

• 1.  
  Type I collapsars are normal GRBs produced by massive stars that form direct-collapse black holes. These systems do form a proto-neutron star, but the standard convective engine is unable to produce a supernova explosion. The GRB-associated supernova is produced by the GRB itself. This is our standard GRB model.
• 2.  
  Type II collapsars are GRBs produced by stars who first produced weak supernova explosions. The subsequent fallback produces a black hole. We will assume these systems make dim GRBs because we expect their accretion rates to be lower than the accretion rates of Type I collapsars. These bursts will only be detected at low redshift, and we do not include these systems in the study of high-redshift GRBs.
• 3.  
  Type III collapsars are GRBs produced in the most massive stars (above the pair-instability gap). As we shall discuss below, these collapsars may be much dimmer than normal GRBs because these systems collapse from large proto-black holes (Fryer et al. 2001), and the engine weakens significantly with increasing black hole mass.

Popham et al. (1999) estimated the strength of a Blandford–Znajek jet based on the properties of the black hole and accretion disk. In this paradigm, the magnetic field develops in the disk and its strength is relative to the disk energy. This is directly proportional to the accretion rate but also depends on the spin and mass of the black hole. We use the following fit to the data from Heger et al. (2003):

$$\begin{eqnarray}\begin{array}{rcl}{L}_{\mathrm{BZ}} & = & {f}_{\mathrm{Lpeak}}{10}^{50}{a}_{\mathrm{spin}}^{2}{10}^{0.1/(1-{a}_{\mathrm{spin}})-0.1}\\ & & \times {\left(\displaystyle \frac{{M}_{\mathrm{BH}}}{3{M}_{\odot }}\right)}^{-3}\displaystyle \frac{{dM}/{dt}}{0.1{M}_{\odot }\,s-1}\,\mathrm{erg}\,{{\rm{s}}}^{-1},\end{array}\end{eqnarray} \tag{ 15 }$$

where M_BH and a_spin are the black hole mass and spin, respectively, dM/dt is the accretion rate, and f_Lpeak is a parameter that covers the efficiency at which disk energy is converted to the jet power and, ultimately, the GRB luminosity.

With this formula, we can derive the power in the jet and, if we assume there is a one-to-one correspondence between jet power and GRB luminosity (assuming some efficiency less than 1), the GRB luminosity as a function of time (Figure 4). There is a generic trend in the luminosity where the power in the jet decreases with time as the black hole accretes mass and the event horizon moves outward. Another feature of the models is that, at high metallicities, the more massive stars cannot drive strong jets. This is because mass loss removes the angular momentum of the star and the low black hole spin limits the power from the Blandford–Znajek mechanism.

Zoom In Zoom Out Reset image size

As we move to high redshift, the metallicity decreases (more massive stars do not lose their angular momentum) and the IMF flattens out. If we include both effects, using our metallicity versus redshift relationship, we can study the evolution of the distribution of GRB luminosities as a function of redshift. Figure 5 shows these peak luminosity distributions as a function of redshift. In our f_IMF = 0.4, z₀ = 5, f_bin = 0 model, the IMF flattens, producing a larger number of bright GRBs.

Zoom In Zoom Out Reset image size

Variations in the spin rate and initial seed magnetic fields between different stars should broaden this distribution. Without a better theoretical understanding of stellar evolution, this is difficult to predict. We study these effects (see Section 4.1) by assuming that our estimates produce average luminosities and using observed distributions as a guide (Salvaterra et al. 2012; Pescalli et al. 2015). These observational studies suggest a distribution $f{(L)\propto (L/{L}_{\mathrm{break}})}^{-a}$, where a = 1.5 and 2.3 below and above the break luminosity (L _break). In this distribution model, our nominal luminosity corresponds to the "break" luminosity in these studies, and the variation in the luminosity with redshift becomes even less pronounced. Although these observations also predict a large number of dim bursts (which may be caused by lower-mass progenitors with fallback black holes), our focus in our study is on bright GRBs, and we limit our distribution to these bright bursts.

In principle, a changing IMF may also change the level of beaming, if the level of beaming of GRBs depends on the mass of the progenitor. While at this time there is no clear theoretical argument to suggest a variation in beaming with stellar mass, some observations suggest a mild beaming dependence on the GRB luminosity: ${{\rm{\Theta }}}_{\mathrm{beaming}}\propto {L}_{\mathrm{GRB}}^{0.2}$ (Ghirlanda et al. 2007), which corresponds to a small difference in opening angle with redshift. There is also empirical evidence for a stronger decrease in jet-opening angle with increasing redshift, possibly arising from the higher pressure in the envelopes of lower-metallicity stars that is expected due to their lower opacity and corresponding higher density (Lloyd-Ronning et al. 2020). If GRBs do, in fact, become more beamed with redshift, we would expect a lower GRB rate with redshift than the results from our models without beaming evolution.

We can estimate the duration as well based on the time it takes for the GRB jet luminosity to drop one order of magnitude. The primary constraint on the duration of the GRB is the accretion rate. If the black hole accretes mass quickly, its event horizon will move outward, reducing the energy in the jet under this disk-driven jet paradigm. What this means is that, even if the mass accretion rate remains steady for long periods of time, the GRB duration will drop off. As with the GRB luminosity, the difficulty in calculating the duration is tying the power source of this jet to a GRB duration. Figure 6 shows the distribution of the duration of the power source (timescale for the luminosity to drop by an order of magnitude) as a function of redshift. Because of the strong dependence of the luminosity on the black hole mass, there is very little evolution in this duration. This duration corresponds to the duration of the energy source near the black hole. Although it is likely that there exists some correlation between this energy source and burst duration, its relationship with the burst duration, e.g., T₉₀ is not well understood at this time.

Zoom In Zoom Out Reset image size

A supernova-like transient is also produced when the GRB jet disrupts the star. It drives a strong, asymmetric explosion called a hypernova (Iwamoto et al. 1998). Just as the structure of the star alters the conditions around a BHAD, the compactness of the core can alter the transient produced when the jet disrupts the star. To first order, Type Ic hypernovae rely on the same physics (transport of energy from a ⁵⁶Ni power source—but note that shock heating may also be important) as Type Ia supernovae and can be understood from the physical understanding from studies of those supernovae. The Arnett law describes the peak luminosity as (Arnett 1980)

$$\begin{eqnarray}&&L=2.055\times {10}^{10}{L}_{\odot }{M}_{\mathrm{Ni}}{\rm{\Lambda }}(x,y),\end{eqnarray} \tag{ 16 }$$

where M_Ni is the nickel mass and Λ(x, y) is a dimensionless function that depends upon a shaping parameter $y={(\kappa {M}_{\mathrm{ejecta}}/{v}_{\mathrm{ejecta}})}^{1/2}$, where κ, M_ejecta, and v_ejecta are the opacity, mass, and velocity of the ejecta. Under this formalism, the peak luminosity depends on the following properties of the explosion: ejecta mass, the nickel mass, the composition (for the opacity), and the explosion energy (${E}_{\exp }$).

To understand these dependencies, we have run a series of calculations of supernovae varying the ejecta mass, nickel mass, explosion energy, and radius. For these calculations, we use the gray diffusion supernova code described in De La Rosa et al. (2017), designed for large grids of supernova calculations. It includes a flexible framework for initial conditions (to probe the dependencies on the progenitor and explosion properties) and a recipe to include shock heating. Figure 7 presents the results of some of these explosions, showing trends in nickel and ejecta masses as well as explosion energy. From these calculations, we have derived fits to the peak luminosity and duration (defined as the time spent within one e-folding of the peak luminosity):

$$\begin{eqnarray}&&{L}_{\mathrm{peak}}\propto {M}_{\mathrm{ejecta}}^{-0.43}{M}_{\mathrm{Ni}}{E}_{\exp }^{0.024},\end{eqnarray} \tag{ 17 }$$

where M_ejecta is the ejecta mass, M_Ni is the nickel mass, and ${E}_{\exp }$ is the explosion energy. The corresponding

$$\begin{eqnarray}&&{t}_{\mathrm{dur}}\propto {M}_{\mathrm{ejecta}}^{0.68}{M}_{\mathrm{Ni}}^{-1.25}{E}_{\exp }^{0.032}.\end{eqnarray} \tag{ 18 }$$

Zoom In Zoom Out Reset image size

The nickel production in a collapsar is either produced in the disk wind (Surman & McLaughlin 2004) or in the shock as ejecta plows through the star (Nakamura et al. 2001). In both cases, it depends on the density, and hence mass, of the core. For a zeroth-order approximation of this dependence, we have assumed that the nickel mass, M_Ni, is proportional to CO core mass, M_CO. By combining this assumption with our fitted formulae for the supernova peak luminosity and duration (Equations (17) and (18)), our supernova properties depend only on the CO core mass and the explosion energy. Using our stellar and explosion properties discussed earlier in the paper (Equations (9) and (15)), the metallicity dependence on redshift (Equation (11)), and the IMF dependence on redshift (Equation (6)), we can model the distribution of both peak luminosities and supernova durations. Figure 8 shows the distribution of peak supernova luminosities as a function of redshift for two different models from the evolution of the IMF. The corresponding distributions for the supernova duration (defined by the time when the luminosity is within an e-folding of the peak flux) are shown in Figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although this discussion focused on GRB-associated hypernovae, we note that hypernovae may also be produced in systems where the jet fails to produce a strong burst of gamma-rays (i.e., the jet decelerates before exiting the star). The rate of hypernovae is difficult to determine. Guetta & Valle (2007) found that roughly 7% of all Type Ib/c supernovae exhibited broad lines characteristic of hypernovae and that the observations of hypernovae suggested a relative rate of normal GRBs to hypernovae is somewhere between 4% and 30%. This difference could be explained simply by the beaming fraction of GRBs and Guetta & Valle (2007) argued that the hypernova rate was consistent with the normal GRB rate. For our models, we assume there are twice as many hypernovae as GRBs (presumably formed from dim GRBs or failed-jet hypernovae).

Other types of massive-star explosions also evolve with a changing IMF, and by combining GRB results with observations of other powerful explosions, we can further constrain the IMF Lazar & Bromm (2022). Most notable are the pair-instability supernovae produced in a runaway thermonuclear burning phase instigated by electron/positron pair production in a hot core (Barkat et al. 1967). Although it appears the exact stellar masses that undergo pair instability is sensitive to the ¹² C(α, γ) rate (Costa et al. 2021; Farmer et al. 2020), we use the limits and metallicity dependencies set by Fryer et al. (2001) and Heger et al. (2003). At high metallicities, winds reduce the stellar mass, preventing the high temperatures needed to make pairs. We include these mass-loss effects, reducing the rate of pair-instability supernovae, but we assume that, if an explosion is produced, its yields and luminosity do not change.

Our standard models do not include the potentially large fraction of dim GRBs. At very low redshift, these dim GRBs can dominate the observed GRB sample. We have augmented are standard model by extending this dim portion of the GRBs. The luminosity function for dim GRBs is quite uncertain from both the theoretical and observational points of view. It is difficult to constrain the nature of these dim bursts. Dim GRBs could be produced for a wide range of reasons (e.g., off-angle jets or failed jets), and with our lack of understanding of jet physics, it is difficult to theoretically determine both the rate and properties of these dim GRBs. Observational constraints are also difficult because these dim GRBs are difficult to constrain because their flux can be mostly below the detection thresholds of the available observatories. To have a crude estimation of the dim GRB contribution in the BAT sample, we explore a dimmer luminosity function by lowering the fiducial luminosity to f_fiducial = 10⁵⁰ erg s⁻¹ for dim GRBs.

Figure 14 shows the expected detections of dim GRBs in the BAT sample based on these two luminosity functions. Based on our current assumptions, results show that these dim GRBs should consist of no more than a few detections per year in the BAT sample (i.e., ≲1%). Dim GRBs are difficult to detect but should also produce hypernovae. This places some limits on the relative fraction of dim GRBs to normal GRBs at low redshift. With these constraints, we do not expect too many dim GRBs in the observed GRB sample.

Zoom In Zoom Out Reset image size

Using the calibrated values of L_fiducial and f_GRB, we estimate detection rates of high-redshift GRBs from a next-generation detector. For this study, we use the characteristics of the newly proposed Gamow Explorer (White et al. 2021), a proposed multiwavelength space telescope aiming to explore the early universe with GRBs. The instrumental design provided by the Gamow team includes the effective area as a function of energy, the background count rate of $1.632\times {10}^{-5}\ \mathrm{count}\ {{\rm{s}}}^{-1}\ {\mathrm{arcmin}}^{-2}$, and the PSF area of $721\ {\mathrm{arcmin}}^{2}$. The Gamow photon counts for each burst C_GRB is estimated using the Gamow effective area and the burst energy spectrum. The total background count during a burst duration is estimated as C_bkg = 1.632 × 10⁻⁵ × 721 × T₉₀, where T₉₀ is randomly selected from the BAT-detected GRBs. A burst would be considered detected if the source count C_GRB is larger than the photon count that corresponds to a probability of <10⁻¹⁰ from Poisson fluctuations of the total background count C_bkg. Moreover, a minimum source count of 5 is required in order to take into account additional systematic noise.

Figure 15 shows the estimated Gamow GRB detection rate for these models, and Table 2 summarizes the accumulated number of GRB detections for different redshifts. The differences between these different models show the true strength of high-redshift observations. For example, for a fiducial luminosity of 10^51.5 erg, the detection rate above a redshift of 5 ranges from 4 to 40. Above a redshift of 8, this rate ranges from a few tenths to 10. Even a null detection by Gamow will constrain our understanding of our models and the evolution of the IMF.

Zoom In Zoom Out Reset image size

Figure 16 shows an example of the expected constraints on parameters related to IMF evolution that could be placed by the Gamow high-z GRB survey. We assume values for the evolution of the IMF based on the current best estimates of the redshift evolution (Pescalli et al. 2016). A good fit to that redshift evolution finds redshift parameters of z₀ = 2.5 and f_IMF = 0.15 from Equation (6). With the current limited high-redshift GRB observations, it is difficult to truly measure this number. To determine how the added statistics from a mission like Gamow can constrain these IMF evolution parameters, we use a suite of simulations with gradual changes of f_IMF from 0.0 to 1.0 and z₀ from 2.0 to 6.0. For these additional simulations, all other parameters are fixed to the same values. That is, we use the star formation rate from Madau & Fragos (2017), a fiducial luminosity of 10^51.5 erg s⁻¹ cm², binary fraction f_bin = 0.01, and wind parameter f_wind = 45. For each of these models, we estimate the Gamow GRB detections at different redshifts and compare that with the expected Gamow GRB detection rates based on current GRB observations (i.e., the blue regions in Figure 15). The f_IMF and z₀ values that give a GRB detection rate within the blue region are included in the constraint area shown in Figure 16.

Zoom In Zoom Out Reset image size

In addition to predicting the evolution (both the rate and properties) of GRBs, our models predict the evolution with redshift (due to mass loss and IMF evolution) of other transients produced by massive stars, e.g., pair-instability supernovae and hypernovae. Figure 17 shows the rate of pair-instability supernovae and hypernovae with respect to the supernova rate as a function of redshift. As with GRBs, both of these transients are reduced at high metallicity because of mass loss from metallicity-driven winds. At low redshifts, there is an increase in the rates of these systems because of the change in winds. But at high redshifts, the rate can increase even further as the IMF flattens. In our most extreme case, the pair-instability rate ultimately becomes more than 10 times higher than the supernova rate. In this section, we discuss this impact and the potential for a broader range of observational constraints on this redshift evolution. These transients are among the brightest astrophysical outbursts, allowing them to be observed to large distances (Whalen et al. 2014; Wiggins et al. 2015). It is possible that the James Webb Space Telescope (Gardner et al. 2006) will be able to observe pair-instability supernovae out to redshift 10 (Whalen et al. 2014).

Zoom In Zoom Out Reset image size

Hypernovae and pair-instability supernovae produce different nucleosynthetic yields than normal supernovae, and, if the IMF flattens, these unique nucleosynthetic signatures can be compared to observations (e.g., in metal-poor stars) to constrain the flattening of the IMF at high redshift. A number of studies have calculated the yields of hypernovae and pair-instability supernovae (Nakamura et al. 2001; Heger & Woosley 2002; Umeda & Nomoto 2002; Chen et al. 2014; Kozyreva et al. 2014). Using the yields from Umeda & Nomoto (2002), we can calculate the nucleosynthetic yields as a function of redshift for our different models. Figure 18 shows the yields (abundance divided by iron production: [X/Fe]) for four different IMF-evolving models divided by the flat [X/Fe] ratios. Per star, more iron-peak elements are produced for these models with flattened IMFs.

Zoom In Zoom Out Reset image size

Other studies have attempted or identified approaches to use metal-poor stars to search for evidence of hypernovae and pair-instability supernova explosions (e.g., Umeda & Nomoto 2002; Tumlinson et al. 2004; Ren et al. 2012). In some cases, scientists have argued that the yields of these stars matched hypernovae with little indication for evidence of signatures from pair-instability supernovae (Umeda & Nomoto 2002; Tumlinson et al. 2004; Takahashi et al. 2018). The fact that, thus far, no clear signature of pair-instability supernovae exists either suggests that we do not know the exact range of masses for pair-instability supernovae or that the IMF does not flatten as much as suggested in our most extreme model. But much more work must be done to move from finding evidence of these extreme explosions to probing the IMF at high redshifts from these nucleosynthesis observations.

Another probe of the nuclear yields at high redshift is to identify high-redshift gas (e.g., through active galactic nuclei or GRBs) and, through absorption-line spectroscopy, probe the nucleosynthetic yields. Missions like Gamow are designed to obtain these absorption-line spectra to constrain these yields.

An additional diagnostic of an evolving IMF is the number of ionizing photons produced by stars. Stellar model predictions of the number of ionizing photons produced during the lifetime of massive stars depend on rotation and metallicity, but these differences are actually less than differences produced by different stellar codes (Topping & Shull 2015). However, with all models, there is a clear trend where the number of ionizing photons produced increases with stellar mass. To study the effect of our evolving IMFs, we fit the models in Topping & Shull (2015):

$$\begin{eqnarray}&&{Q}_{0}={Q}_{\mathrm{base}}[1+{f}_{\mathrm{ion}}({M}_{\mathrm{star}}-20)],\end{eqnarray} \tag{ 19 }$$

where Q₀ is the lifetime integrated H-ionizing photons for a star of mass (M_star), Q_base is roughly 10⁵³ photons, and f_ion is determined by fitting to the data in Topping & Shull (2015) (f_ion = 0.07 − 0.15). Figure 19 shows the ratio of this number of photons with redshift for four of our evolutionary models.

Zoom In Zoom Out Reset image size

An increase in ionizing photon production from an evolving IMF could relieve tension in reionization models where large escape fractions of ionizing photons from galaxies (>10%) are required at high redshift to complete reionization by z ≈ 5.5–6. Observations of galaxies at z < 4, where escape fractions of ionizing photons can be measured, predominantly report low escape fractions (<10%) (Siana et al. 2010; Sandberg et al. 2015; Grazian et al. 2017; Rutkowski et al. 2017; Pahl et al. 2021; Saxena et al. 2022), though recently higher escape fractions have been measured in some galaxies (Steidel et al. 2018; Fletcher et al. 2019; Izotov et al. 2021; Marques-Chaves et al. 2021; Saxena et al. 2022). Higher Lyman-continuum photon production efficiencies in galaxies at high redshift are favored in Finkelstein et al. (2019) where they consider scenarios for reionization with empirically supported low galaxy ionizing photon escape fractions.

On the other hand, if Lyman-continuum photons are too numerous at high redshift, then reionization will complete sooner than observed quasar sight lines suggest and predicted opacities in the IGM may be too low. This may put constraints on the more extreme models shown in Figure 19, barring a lower escape fraction. However, because the uncertainties in reionization, such as the escape fraction and the emissivity of ionizing photons, are degenerate the incidence of GRBs at high redshift may help constrain each of these properties.

Finally, energetic explosions such as hypernovae may have an enhanced dust production rate at high redshift (Brooker et al. 2021). There is some evidence of this increased dust production at low metallicities (Nanni et al. 2020). Figure 20 shows the total dust (plus two specific dust grains) production using the models from Brooker et al. (2021) and our IMF/metallicity evolution. Initially, the dust production increases before decreasing as we make increasingly massive stars. This production does not include any contribution from pair-instability supernovae, which may be a major contributor to high-redshift dust production (Schneider et al. 2004). If they produce large amounts of dust, these results would have to be revised at high redshift.

Zoom In Zoom Out Reset image size

Coupled with a larger sample of GRB observations probing the same redshift space, this broad suite of observations will be able to set firm limits on the evolution of the IMF at high redshift.