Long gamma-ray burst rate at very high redshift

TL;DR

This paper predicts the detection rates of high-redshift long gamma-ray bursts by future missions, highlighting the potential to observe Population III star-originated GRBs and their implications for early universe star formation.

Contribution

It introduces a model combining Pop I&II and Pop III star contributions to estimate high-redshift GRB detection rates for upcoming missions.

Findings

Pop I&II GRB rate decreases steeply with redshift.

Future missions could detect Pop III GRBs at z>8.

Detection of Pop III GRBs can reveal early star formation history.

Abstract

Future missions for long gammma-ray burst (GRB) observations at high redshift such as HiZ-GUNDAM and THESEUS will provide clue to the star formation history in our universe. In this paper focusing on high redshift (z>8) GRBs, we calculate the detection rate of long GRBs by future observations, considering both Population (Pop) I&II stars and Pop III stars as GRB progenitors. For the Pop I&II star formation rate (SFR), we adopt an up-to-date model of high-redshift SFR based on the halo mass function and dark matter accretion rate obtained from cosmological simulations. We show that the Pop I&II GRB rate steeply decreases with redshift. This would rather enable us to detect the different type of GRBs, Pop III GRBs, at very high redshift. If 10% or more Pop III stars die as an ultra-long GRB, the future missions would detect such GRBs in one year in spite of their low fluence. More luminous GRBs are expected from massive compact Pop III stars produced via the binary merger. In our conventional case, the detection rate of such luminous GRBs is 3-20 /yr (z>8). Those future observations contribute to revealing of the Pop III star formation history.

Tables

Table 1. Table 1 The initial distribution functions in This paper.

IMF	Initial Mass Ratio Function	Initial Period Function	Initial Eccentricity function
flat	flat	1/a	e
		$a_{\min }$*

Table 2. Table 2 The long GRB fraction of Pop III f GRB subscript 𝑓 GRB f_{\rm GRB} which consist of highly rotating Helium stars and highly rotating main-sequence stars

progenitor type	Highly rotating Helium stars	highly rotating main-sequence stars
${\alpha }_{\mathrm{CE}}\lambda =1$ model	1.1%	2.5%
${\alpha }_{\mathrm{CE}}\lambda =0.1$ model	1.6%	7.8%

Equations

$${\rho }_{\mathrm{SFR}}$$

$$=$$

$$SFR(t)=\frac{d{M}_{\ast }}{dt}=\frac{d{M}_{\ast }}{d{M}_{\mathrm{h}}}\dot{M_{\mathrm{h}}}=\frac{\alpha M_{\ast }}{M_{\mathrm{h}}}\dot{M_{\mathrm{h}}},$$

$$\alpha =\frac{d{M}_{\ast }}{d{M}_{\mathrm{h}}}\frac{M_{\mathrm{h}}}{M_{\ast }}=\frac{SFR/M_{\ast }}{\dot{M_{\mathrm{h}}}/M_{\mathrm{h}}}=\frac{\mathrm{sSFR}}{\mathrm{SMAR}}.$$

$$\alpha =\frac{\frac{d\log M_{\ast }}{dt}}{\frac{d\log M_{\mathrm{h}}}{dt}}=\frac{d\log M_{\ast }}{d\log M_{\mathrm{h}}}.$$

$$SFR=\frac{2\times 1.7\times 10^{-2}}{(M_{\mathrm{h}}/M_{\mathrm{br}}{)}^{-1.1}+(M_{\mathrm{h}}/M_{\mathrm{br}}{)}^{0.3}}\dot{M_{\mathrm{h}}},$$

$$\log ({\rho }_{\mathrm{SFR}}/[M_{\odot }\mathrm{y}{\mathrm{r}}^{-1}\mathrm{Mp}{\mathrm{c}}^{-3}])=-0.32z+\mathrm{0.47.}$$

$$\frac{d{N}_{\mathrm{GRB}}}{d{\Omega }_{\mathrm{obs}}d{t}_{\mathrm{obs}}}(z>z_0)={\int }_{z_0}^{\infty }\frac{{\dot{n}}_{\mathrm{GRB}}}{1+z}\frac{d{V}_{\mathrm{c}}}{dzd{\Omega }_{\mathrm{obs}}}dz,$$

$$\frac{d{V}_{\mathrm{c}}}{dzd{\Omega }_{\mathrm{obs}}}=\frac{D_{\mathrm{L}}^2}{(1+z{)}^2}\frac{c}{H_0}\frac{1}{\sqrt{{\Omega }_{\mathrm{m}}(1+z{)}^3+{\Omega }_{\Lambda }}},$$

$$\displaystyle\frac{dn_{\rm GRB}}{dL}\propto\left\{\begin{array}[]{ll}\left(\frac{L}{L_{\star}}\right)^{-a},&\mbox{for}~{}L\leq L_{\star}\\ \left(\frac{L}{L_{\star}}\right)^{-b},&\mbox{for}~{}L>L_{\star}.\end{array}\right.$$

$${\epsilon }_{\mathrm{pk}}=1.8\times {\left(\frac{1}{2.34\times 10^{-5}}\times \frac{L}{10^{52}\text{mbox}erg\text{mbox}{s}^{-1}}\right)}^{0.5}\text{mbox}keV.$$

$$\frac{d{N}_{\text{mbox}III}}{dM}=\text{mbox}const.$$

$${\alpha }_{\mathrm{CE}}\left(\frac{G{M}_{\mathrm{c},1}{M}_2}{2a_{\mathrm{f}}}-\frac{G{M}_1{M}_2}{2a_{\mathrm{i}}}\right)=\frac{G{M}_1{M}_{\mathrm{env},1}}{\lambda R_1},$$

$${\alpha }_{\mathrm{CE}}\left(\frac{G{M}_{\mathrm{c},1}{M}_{c,2}}{2a_{\mathrm{f}}}-\frac{G{M}_1{M}_2}{2a_{\mathrm{i}}}\right)=$$

$$+\frac{G{M}_2{M}_{\mathrm{env},2}}{\lambda R_2},$$

$$\dot{M}=\min \left(\frac{3}{10}\frac{M}{{\tau }_{\mathrm{KH}}},10^{-14}{\left(1-\frac{\Omega }{{\Omega }_{\mathrm{K}}}\right)}^{-0.43}\right)$$

$$\dot{J}=\frac{2}{3}\dot{M}{R}^2\Omega .$$

$$R_{\mathrm{GRB}}=f_{\mathrm{B}}\cdot f_{\mathrm{GRB}}\cdot \left(\frac{f_{\mathrm{b}}}{1+f_{\mathrm{b}}}\right)\cdot \frac{SFR}{M_{\mathrm{III}}}.$$

Full text

1. introduction

The gravitational collapse of massive stars is considered as the trigger of the long gamma-ray burst (GRB), which is confirmed by observing long GRBs associated with hypernovae (Galama et al., 1998; Hjorth et al., 2003; Melandri et al., 2014). However, the long GRB rate is not simply proportional to the star formation rate (SFR) (e.g. Wanderman & Piran, 2010; Lien et al., 2014, 2015). Long GRBs, which are the most luminous astrophysical phenomenon, are a powerful tool to reveal the high redshift star formation especially at $z\gtrsim 10$. The long GRB at the highest redshift ever observed is GRB 090429B with a photometric redshift $z\sim 9.4$ (Cucchiara et al., 2011). There are future plans of long GRB observations at high redshift such as HiZ-GUNDAM (Yoshida et al., 2016) and THESEUS (Yuan et al., 2016; Amati et al., 2018; Rossi, Stratta, Maiorano, Amati, Nicastro & Palazzi, 2018). Those observational missions can probe the high redshift universe. We have to identify what type of stars become long GRB progenitors to calculate the long GRB rate at high redshift. At present, the prediction of the GRB rate at $z\gtrsim 10$ is difficult because of the lack of the observational knowledge of GRBs and the SFR at $z\gtrsim 10$.

As the long GRB progenitor at high redshift, many authors have considered not only Population I and II (Pop I&II) stars, but also Population III (Pop III) stars (e.g. Lloyd-Ronning et al., 2002; Bromm & Loeb, 2006; Belczynski et al., 2007; Campisi et al., 2011; de Souza et al., 2011; Toma et al., 2011, 2016; Ghirlanda et al., 2015; Burlon et al., 2016). Pop III stars are first stars formed from the primordial gas with no metal. Pop III stars are more massive stars than Pop I&II stars due to a lack of effective coolant such as metal and dusts (Omukai et al., 2005; Dayal & Ferrara, 2018). It has been suggested that two different Pop III star formation modes. The first generation of Pop III stars (Pop III.1) are formed from the primordial gas unaffected by the previous star formation, where the main coolant is H2 molecule (Tegmark et al., 1997; Abel et al., 2002; Bromm et al., 2002; Yoshida et al., 2006). Recent studies suggest that the radiation feedback from the massive protostar leads to the typical mass of Pop III.1 of $\sim 40~{}{M}_{\odot}$ (Hosokawa et al., 2011). The second generation Pop III stars (Pop III.2) are formed from the no metal gas ionized by the radiation from the previous star formation (Johnson & Bromm, 2006; McKee & Tan, 2008). In such ionized gases, the hydrogen dueteride (HD) cooling is more effective than H2 cooling, so that the typical mass ($\sim 20~{}{M}_{\odot}$) is slightly less than the mass of Pop III.1 (Hosokawa et al., 2012).

Some simulations show that disk fragmentations frequently occur, which implies the existence of binary Pop III stars (Saigo et al., 2004; Machida et al., 2008; Stacy & Bromm, 2013; Susa et al., 2014). Kinugawa et al. (2014) predicted detection of gravitational waves (GWs) from binary black hole (BH) mergers originated from Pop III stars. The first GW detection with LIGO, GW150914 (Abbott et al., 2016a), was a 30+30 ${M}_{\odot}$ binary black hole merger, which supports the BH binary formation from Pop III stars.

The Pop III star formation rate has been studied with semi-analytical method or numerical simulations (e.g. de Souza et al., 2011; Johnson et al., 2013). The Thomson scattering optical depth for cosmic microwave background photons measured with Planck (Planck Collaboration et al., 2014, 2016a, 2016b) is lower than the previous values measured with WMAP (Dunkley et al., 2009). This provides tight constraints on the star formation history at high redshift (Visbal et al., 2015; Hartwig et al., 2016; Inayoshi et al., 2016).

Since the mass distribution of Pop III stars is biased to heavier range than those for Pop I&II stars, we can expect that Pop III stars tend to launch long GRBs easily. de Souza et al. (2011) calculated the Pop III SFR using a semi-analytical method and estimated the Pop III GRB rate. Yoon et al. (2012) calculated the rotating Pop III stellar evolution as a GRB progenitor model and showed that rotating massive Pop III stars experience the chemically homogeneous evolution and can launch GRB jets at the final stage of their evolution. On the other hand, Nakauchi et al. (2012) showed that Pop III blue supergiant stars, which hold a massive hydrogen envelope, may give rise to a GRB with a duration of $\sim 10^{5}$ s in the observer frame because of the long mass accretion phase. The peak luminosity of the Pop III ultra long GRBs was estimated as $\sim 5\times 10^{50}~{}\rm erg~{}s^{-1}$. For binary Pop III stars, Belczynski et al. (2007) discussed the tidal spin-up and the envelope ejection by the binary interaction and calculated the number of GRB progenitors based on a criterion from the angular momentum and amount of envelope.

In this paper focusing on high redshift ($>8$) GRBs, we calculate the detection rate of long GRBs by future observations, considering both Pop I&II stars and Pop III stars as GRB progenitors. The future detections of GRBs at very high redshift ($z\gtrsim 10$) will unveil the star formation history in the very early era. For Pop I&II stars, we adopt an up-to-date model of high-redshift SFR proposed by Harikane et al. (2018) with the halo mass function and dark matter accretion rate obtained from cosmological simulations in Ishiyama et al. (2015).

On the other hand, in the case of Pop III stars, we consider both ultra-long GRBs from massive stars with heavy envelope and classical GRBs from massive compact stars, which experienced the binary merger.

The long GRB rate from Pop I&II stars at high redshift is discussed in §2. In §3, we consider the long GRB rate from Pop III stars considering two cases: the ultra-long GRB and classical GRB progenitors. Our results are summarized in §4.

2. GRB from Pop I&II stars

2.1. Pop I&II star formation rate

In order to calculate the long GRB rate from Pop I&II stars at high redshift, first we estimate the cosmic star formation rate densities (SFRDs) of Pop I & II stars. Since the cosmic SFRDs at $z>10$ is poorly constrained from observations, we consider two models in Behroozi & Silk (2015) and Harikane et al. (2018) for SFRDs. Below we shortly review their calculations. See Behroozi & Silk (2015) and Harikane et al. (2018) for more details.

In both the models, the SFRD $\rho_{\rm SFR}$ can be calculated based on the following equation:

[TABLE]

where $M_{\mathrm{h}}$, $\frac{dn}{dM_{\mathrm{h}}}$, and $\dot{M_{\mathrm{h}}}$ are dark matter halo mass, halo mass function, and dark matter accretion rate, respectively. In Behroozi & Silk (2015), they obtain halo mass function and the dark matter accretion rate at $z=5-15$ from the Bolshoi N-body simulation (Klypin et al., 2011). The Bolshoi simulation is calculated in the redshift range of $z=0-80$ in a 250 $h^{-1}$ Mpc with the mass resolution of $1.9\times 10^{8}\ M_{\odot}$. The SFR per halo $SFR$ in Behroozi & Silk (2015) is expressed as

[TABLE]

where $M_{*}$ is the stellar mass, and $\alpha$ is the ratio of the specific star formation rate (sSFR) and the halo specific mass accretion rate (SMAR), expressed as,

[TABLE]

In the Behroozi & Silk (2015), they assume that $\alpha$ remains constant at $z\geq 5$ over the galaxy’s star formation history. The constant $\alpha$ is motivated because the galaxy’s history would be dominated by a single feedback mode; they consider only the stellar feedback (supernovae and reionization) neglecting the AGN feedback. This assumption also implies a relation $M_{*}\propto M^{\alpha}_{\mathrm{h}}$, because equation (4) is interpreted as

[TABLE]

Based on this assumption and simulation results for halo mass function and dark matter accretion rate, Behroozi & Silk (2015) start their calculation from $z=5$ using the abundance matching results for $\alpha$ and the $M_{*}$-$M_{\mathrm{h}}$ relation at $z=5$. Figure 1 shows their SFRDs for $z=5-15$. Their SFRDs agree well with the observations at $z=5-9$, while they are slightly higher than the recent estimates at $z=10$ based on UV luminosity function observations (Bouwens et al. 2015; Oesch et al. 2018).

On the other hand, the clustering analysis of $z\sim 4-7$ Lyman-break galaxies in Harikane et al. (2018) and Harikane et al. (2016) provides an empirical equation for the SFR per halo as

[TABLE]

where $M_{\rm br}=10^{11.35}M_{\odot}$. We use the halo mass function, $\frac{dn}{dM_{\mathrm{h}}}$, and the dark matter accretion rate, $\dot{M_{\mathrm{h}}}$, at $z=5-30$ from the Phi-1 simulation in Ishiyama et al. (2015) 111http://hpc.imit.chiba-u.jp/~ishiymtm/db.html, which is calculated in the redshift range of $z=0-30$ in a $32\ h^{-1}\mathrm{Mpc}$ box with the mass resolution of $3.28\times 10^{5}\ M_{\odot}$. Hereafter, we call this model HPhi-1.

Figure 1 shows the cosmic SFRDs at $z=5-30$ in the HPhi-1model. In the redshift range of $z=7-25$, the model is well fitted by a simple power law function,

[TABLE]

As shown in Figure 1, the HPhi-1 model predicts lower SFRDs than previous models, but the results agree with UV luminosity function observations at $z=5-10$ (Bouwens et al., 2015; Oesch et al., 2018). For comparison, we also plot the extrapolated fitting formula of Madau & Dickinson (2014) in Figure 1. This fitting formula was derived from the observation data at $0<z<8$. Since the extrapolated fitting formula for $z>8$ is not consistent with the observation results, hereafter, we focus on HPhi-1 model, and Behroozi & Silk (2015) which are based on the halo mass function and dark matter accretion rate obtained from cosmological simulations.

2.2. GRB rate from Pop I&II stars

The redshift evolution of the GRB rate does not follow the star formation rate (e.g. Wanderman & Piran, 2010; Lien et al., 2014, 2015). This discrepancy may be explained by the metallicity effect on the progenitor formation (Yoon & Langer, 2005; Hirschi et al., 2005; Woosley & Heger, 2006; Yoon et al., 2006; Kinugawa & Asano, 2017). However, we can expect significantly low metallicity at high redshifts, where its effect may be not important any longer. Then, the GRB rate will roughly follow the star formation rate differently from the GRB rate at lower redshift.. Wanderman & Piran (2010) estimated the comoving GRB rate at $z=8$ as $\dot{n}_{\rm GRB}=10^{+12}_{-8.5}~{}\mbox{Gpc}^{-3}~{}\mbox{yr}^{-1}$ above the luminosity of $L=10^{50}~{}\mbox{erg}~{}\mbox{s}^{-1}$. Assuming a broken power-law for the redshift evolution, more recent analysis by Lien et al. (2014, 2015) provided the rate at $z=8$ as $\dot{n}_{\rm GRB}=6.2^{+3.8}_{-5.3}~{}\mbox{Gpc}^{-3}~{}\mbox{yr}^{-1}$.

Assuming that the GRB rate is proportional to the star formation rate, we extrapolate those rates from $z=8$. Given the field of view $\Omega_{\rm obs}$, the GRB occurnece probability at $z>z_{0}$ is

[TABLE]

where $(1+z)$ in the denominator is the time dilation effect due to the cosmological expansion. Using the luminosity distance $D_{\rm L}$, the differential of the comoving volume $V_{\rm c}$ is written as

[TABLE]

where $H_{0}$, $\Omega_{\rm m}$, and $\Omega_{\Lambda}$ are the standard cosmological parameters.

The detection rate depends on the GRB spectrum, luminosity function, and energy band and sensitivity of instruments. The luminosity functions assumed in Wanderman & Piran (2010) and Lien et al. (2014, 2015) are broken power-law,

[TABLE]

The parameters in Wanderman & Piran (2010) are $L_{\star}=10^{52.5}~{}\mbox{erg}~{}\mbox{s}^{-1}$, $a=1.2$, and $b=2.4$, which implies that lower luminosity GRBs dominate the GRB number. The lower limit in luminosity is taken as $10^{50}~{}\mbox{erg}~{}\mbox{s}^{-1}$. Lien et al. (2014, 2015) adopts $L_{\star}=10^{52.05}~{}\mbox{erg}~{}\mbox{s}^{-1}$, $a=0.65$, and $b=3.0$. In this case, GRBs at $L\sim L_{\star}$ dominate the GRB number.

The GRB spectra are well described by the Band function (Band et al., 1993), whose parameters are the peak energy $\varepsilon_{\rm pk}$, low-energy index $\alpha$, and high-energy index $\beta$. While the peak energy in the rest frame is assumed as a constant $511$ keV in the analysis of Wanderman & Piran (2010), Lien et al. (2014, 2015) adopts a modified Yonetoku relation

[TABLE]

In this paper, as a typical value, we fix the indices as $\alpha=-1$ and $\beta=-2.25$.

As a future observation mission, we consider wide field X-ray monitor with Lobster Eye optics, which may be adopted by the missions in HiZ-GUNDAM (Yoshida et al., 2016) or THESEUS (Yuan et al., 2016; Amati et al., 2018; Rossi, Stratta, Maiorano, Amati, Nicastro & Palazzi, 2018). With such an instrument, we can expect a sensitivity of $\sim 10^{-10}~{}\mbox{erg}~{}\mbox{cm}^{-2}~{}\mbox{s}^{-1}$ for 100 s exposure, and a field of view $\sim 0.2$ str 222Private communication with HiZ-GUNDAM working group. See also Yuan et al. (2016). Multiple Lobster Eye systems can enlarge the field of view depending on the budget in future plans. We conservatively assume single Lobster Eye system. . In Figure 2, we plot the expectation of the detection rate adopting the two models: the models in Wanderman & Piran (2010) (circles) and Lien et al. (2014, 2015) (triangles) for the GRB rate, luminosity function and spectral peak energy. If we adopt the star formation rate in Behroozi & Silk (2015) (black and green), both the models suggest a few GRB detection per year for $z>10$. However, the Hphi-1 SFR, which seems consistent with the observed rate at $z\simeq 10$, leads to a detection rate $\ll 1$ for $z>12$ (see red and blue symbols). In addition, we should take into account the efficiency of the redshift confirmation, which depends on the performance of the follow-up infrared telescope onboard HiZ-GUNDAM or THESEUS. The confirmation of GRBs at $z>12$ seems not easy. However, in other words, this provides the opportunity to detect other types of trangient phenomena at high redshifts, such as GRBs from Pop III stars.

3. GRB from Pop III stars

3.1. Pop III star formation rate

Pop III stars are first stars formed from metal-free gases. Pop III stars are massive (e.g. McKee & Tan, 2008; Hosokawa et al., 2011), and have no stellar wind mass loss(Krtička & Kubát, 2006). Since Pop III stars are formed in anomalous circumstance, their formation history may be different from Pop I&II stars.

At present, we have not significant constraint on the Pop III star formation rate from observations. However, the Pop III star formation rate has been estimated using the cosmological simulation. We consider two Pop III SFRDs in de Souza et al. (2011) and Inayoshi et al. (2016).

The SFRD of de Souza et al. (2011) is calculated by a semi-analytical approach, in which they assume that Pop III stars are formed in dark matter haloes at their collapse. They adopt the Sheth-Tormen mass function (Sheth & Tormen, 1999) to estimate the number of dark matter halos at given redshift. They divide the populations into the Pop III.1 and Pop III.2 stars. Pop III.1 stars are the first generation Pop III stars formed in the dark matter ’minihaloes’, where only H2 molecular cooling is the dominant cooling process. Pop III.2 stars are the second generation Pop III stars formed from ionized gases that are in the HII region made by previous Pop III stars or a virialization shock in the halo with the virial temperature $T_{\rm Vir}\gtrsim 10^{4}$ K. The hydrogen deuteride (HD) cooling is efficient below 200 K in this region due to an enhanced free electron fraction. As a result, the typical mass of Pop III.2 is slightly smaller than the Pop III.1 (Hosokawa et al., 2012) but still massive enough to form a BH at their end of life, which may launch a GRB jet.

In order to determine how many Pop III.1 and Pop III.2 stars are formed in collapsed haloes, de Souza et al. (2011) introduce a parameter, the star formation efficiency $f_{*}$. The original SFRD in de Souza et al. (2011) are obtained with $f_{*}=0.1$, and 0.01 for Pop III.1, and Pop III.2, respectively. In de Souza et al (2011), using a criteria based on the virial mass and reionization history, the evolutions of SFRD of Pop III.1 and PopIII.2 are calculated. They also follows the progress of the metal enrichment by protogalactic wind, which prevents the Pop III star formation.

The SFRD in Inayoshi et al. (2016) is also calculated by a semi-analytical approach using the same Sheth-Tormen mass function (Sheth & Tormen, 1999) based on the model described in Visbal et al. (2015). The treatments for the progresses of ionization and metal enrichment are different from the model in de Souza et al. (2011). Inayoshi et al. (2016) take into account the constraint of the star formation from the Thomson scattering optical depth $\tau_{\rm e}=0.066+1\sigma$ where $\sigma=0.016$, which is measured by Planck Collaboration et al. (2016a). These changes generate significant constraint on models of Pop III star formation (Visbal et al., 2015; Inayoshi et al., 2016), although this constraint depends on some parameters such as the escape fraction of photon, the initial mass function (IMF), and so on. Inayoshi et al. (2016) shows that the constraint of the total Pop III star formation density is $\rho_{*,III}\lesssim 6\times 10^{5}~{}{M}_{\odot}\rm~{}Mpc^{-3}$ for $\tau_{\rm e}=0.066+1\sigma$, the escape fraction of ionizing photons from mini-halos $f_{\rm esc}=0.1$, and the flat IMF ($10~{}{M}_{\odot}<M<100~{}{M}_{\odot}$).

The SFRD in de Souza et al. (2011) already conflicts with the limit by the Planck observation so that we adopt the SFRD in de Souza et al. (2011) decreasing by a factor of 0.3, which corresponds to the upper-limit of $\rho_{*,III}$ given by Inayoshi et al. (2016). Figure 3 shows the SFRDs we adopt in this paper based on de Souza et al. (2011) and Inayoshi et al. (2016).

3.2. Ultra-long GRB rate from Pop III

Hereafter, we assume that the GRB rate is proportional to SFRDs in Figure 3. The number of stars are calculated from the initial mass function, which is assumed as the flat

[TABLE]

between $10M_{\odot}$ and $100M_{\odot}$. This implies that the average mass of Pop III stars, $M_{\mbox{III}}$, is $55M_{\odot}$.

In the most optimistic scenario, all such heavy stars give rise to a GRB. Then, the apparent GRB rate is simply estimated as $f_{\rm B}SFR/M_{\mbox{III}}$, where $f_{\rm B}$ is the beaming factor. For metal free stars like Pop III stars, however, the stellar wind is suppressed so that a massive envelope remains at the collapse (Krtička & Kubát, 2006). In such cases, its long free fall time leads to ultra-long GRBs (Nakauchi et al., 2012, 2013), whose duration is $\sim 10^{4}$ s. The jet opening angles in ultra-long GRBs are estimated as wider than $10^{\circ}$ (Levan et al., 2014). Here, we adopt an optimistic opening angle $\theta_{\rm j}=20^{\circ}$, which implies the beaming factor $f_{\rm B}=0.06$. Finally we obtain the comoving GRB rates at $z=8$ as $\dot{n}_{\rm GRB}=1000~{}\mbox{Gpc}^{-3}~{}\mbox{yr}^{-1}$ and $440~{}\mbox{Gpc}^{-3}~{}\mbox{yr}^{-1}$ for the SFRs in de Souza et al. (2011) and Inayoshi et al. (2016), respectively.

The observed typical luminosity of ultra-long GRBs is $\sim 10^{49}~{}\mbox{erg}~{}\mbox{s}^{-1}$ (Gendre et al., 2013; Peng et al., 2013; Levan et al., 2014). As the luminosity function in equation (12), we adopt the same parameters with those in Lien et al. (2014, 2015) but $L_{\star}=10^{49}~{}\mbox{erg}~{}\mbox{s}^{-1}$ with the lower and upper limits $10^{47}~{}\mbox{erg}~{}\mbox{s}^{-1}$ and $10^{51}~{}\mbox{erg}~{}\mbox{s}^{-1}$, respectively. The GRB spectra are assumed as the same Band function with the modified Yonetoku relation. The integration time for an instrument like HiZ-GUNDAM or THESEUS on geocentric orbits would be limitted below $\sim 1000$ s. Assuming a sensitivity of $2\times 10^{-11}~{}\mbox{erg}~{}\mbox{cm}^{-2}~{}\mbox{s}^{-1}$ for 1000 s exposure between 0.5 and 4 keV, and a field of view $\sim 0.2$ str333Private communication with HiZ-GUNDAM working group. See also Yuan et al. (2016). , we plot the expectation of the ultra-long GRB detection in Figure 4.

The estimated detection rates of ultra-long GRBs are significantly higher than the usual GRB rate. The dominant sources in the HiZ-GUNDAM/THESEUS era may be ultra-long GRBs originated from Pop III stars. However, the assumption that all Pop III stars die as a GRB may be too optimistic. Practically the rates in Figure 4 are upper limits of Pop III ultra-long GRB rates. In this optimistic case, ultra-long GRBs from Pop III stars may be detected even at $z=6-8$ especially for the modified de Souza SFRD model. If more than 10% of Pop III stars induce an ultra-long GRB, $\sim$ one detection per year is expected. de Souza et al. (2011) assumed that only 1% of Pop III stars cause a GRB at their end of life. Under such a conservative assumption, detection of an ultra-long GRB seems very difficult.

3.3. Pop III binary population synthesis & binary merger model

Classical long GRBs, whose typical luminosity is more than $10^{52}~{}\mbox{erg}~{}\mbox{s}^{-1}$ with the typical duration of $10~{}\mbox{s}$, are favorable to detect rather than ultra-long GRBs. The binary interaction may produce ideal progenitors to cause classical GRBs.

In order for long GRBs to occur, the progenitors need high angular momentum. Although the angular momentum of Pop III stars is unknown, the remnant of the binary merger during a common envelope (CE) phase possibly have a high angular momentum. When the radius of primary giant suddenly becomes larger or a radical mass transfer shrinks the orbit, the secondary star sometime plunges into the primary envelope. The secondary star spirals in and the envelope of primary will be evaporated. After the CE phase, the binary becomes a close binary which consists of the secondary and the core of the primary giant or the two stars merges during a CE phase. In the latter case, the envelope evaporated, and a highly rotating helium star would remain (Fryer & Heger, 2005). Furthermore, the highly spinning progenitors evolve as chemically homogeneous stars (Yoon et al., 2012). Since such highly rotating stars have small radius, the jet can break out the stellar surface with a high accretion rate like Pop I&II GRBs. For such idealized progenitors, GRBs can be as luminous as usual observed GRBs with duration of $\sim$ 10 s (e.g. Suwa & Ioka, 2011). Thus, we focus on the binary merger model (Fryer & Heger, 2005) as the Pop III GRB progenitor. We consider two channels for Pop III GRB progenitors: (1) highly rotating helium stars and (2) highly rotating main-sequence stars.

The highly rotating helium stars are made by the binary mergers during a CE phase that contains only post main sequence stars. On the other hand, highly rotating main-sequence stars are made by the binary mergers during a CE phase that contains a Giant star and a main sequence star. Using the population synthesis method, we calculate these binary merger fraction of Pop III stars and estimate the Pop III long GRB rate. According to the binary population synthesis method, we set the zero age main sequence binary parameters, such as primary mass $M_{1}$, mass ratio $M_{2}/M_{1}$, separation $a$, and eccentricity $e$, using the initial distribution functions, and calculate each stellar evolution. The numerical code judges whether stars experience the binary interactions (BIs) such as the tidal friction, the mass transfer, the CE phase, and so on and updates the parameters $M_{1},M_{2},a$, and $e$ in each time step. We repeat this calculation using different initial binary parameters chosen by the Monte-Carlo method (Kinugawa et al., 2014). We use the flat IMF that is suggested by some simulations (Hirano et al., 2014; Susa et al., 2014). We assume the other initial distributions are the same as those of Pop I binaries (Heggie, 1975; Abt, 1983; Kobulnicky & Fryer, 2007) as summarized in Table 1.. Using the Monte Carlo method with those initial distribution functions, we calculate the Pop III evolutions of the stellar radius and the core mass and check whether the binary interaction occurs or not. We calculate $10^{6}$ binaries for each models.

The calculation code is the same as the Pop III binary population synthesis code in Kinugawa et al. (2014) and Kinugawa et al. (2017). This code was used to calculate the binary black hole merger rate and the detection rate of LIGO gravitational wave observations. The binary black hole merger rate calculated by this code (Kinugawa et al., 2014, 2016) with a similar initial parameter set to that in this paper is consistent with the LIGO result (Abbott et al., 2016b; The LIGO Scientific Collaboration & the Virgo Collaboration, 2018).

We use the following formalism for the CE phase. The criterion of the mass transfer leading to a CE phase is the same as that of our previous paper (Kinugawa et al., 2014). In order to calculate the separation just after the CE phase $a_{\rm f}$, we use the energy balance prescription (Webbink, 1984)

[TABLE]

for a binary of a giant star and a main sequence star, where $a_{\rm i}$, $R_{1}$, $M_{1}$, $M_{\rm c,1}$, $M_{\rm env,1}=M_{1}-M_{\rm c,1}$, and $M_{2}$ are the binary separation just before the CE phase, the radius, the mass, the core mass and the envelope mass of the giant, and the mass of the companion star, respectively. The value $\alpha_{\rm CE}$ is the efficiency parameter how much the orbital energy can be used in ejecting the envelope. The parameter $\lambda$ is for the envelope binding energy. If the companion star is also a giant, Equation (15) changes into

[TABLE]

where $M_{\rm c,2}$, $M_{\rm env,2}=M_{2}-M_{\rm c,2}$, and $R_{2}$ are the core mass, the envelope mass, and the radius of the companion star, respectively (Dewi et al., 2006). The CE parameters $\alpha_{\rm CE}$ and $\lambda$ are not well understood (Ivanova et al., 2013). We adopt the typical CE parameter values adopted in the previous binary population studies ($\alpha_{\rm CE}\lambda=1$ and $0.1$) (Belczynski et al., 2007; Kinugawa et al., 2014). The simulation of the CE phase is so difficult that the CE parameter are theoretically uncertain (e.g. Ivanova et al., 2013), but there are some observation constraints. Those values ($\alpha_{\rm CE}\lambda=1$ and $0.1$) can reproduce the observation results such as the separation distribution of observed white dwarf binaries, and the period-eccentricity distribution of observed binary pulsars (e.g. Zorotovic, Schreiber, Gänsicke & Nebot Gómez-Morán, 2010; Hijikawa et al., 2019). The Pop III binary black hole merger rates using those values are consistent to the LIGO’s result (Kinugawa et al., 2016; Abbott et al., 2016b). A smaller $\alpha_{\rm CE}\lambda$ leads to a closer separation after the CE phase, and vice versa. Thus, the smaller $\alpha_{\rm CE}\lambda$ implies an efficient stellar merger.Figure 5 shows the $\alpha_{\rm CE}\lambda$ dependence of the binary separation. This figure demonstrates that a smaller $\alpha_{\rm CE}\lambda$ makes binaries easier to merge.

Just after the CE phase, if $a_{\rm f}$ is less than the sum of the radius of the giant’s core and the radius of the companion star (or the sum of radii of the two giants’ cores if the companion star is also a giant), we assume that the binary stars merge.

If two giants merge during the CE phase, the merged helium star obtains a large angular momentum from the orbital angular momentum of the binary. Thus, we assume that the helium star has a spin angular momentum of the Kepler velocity. The mass of the highly rotating helium star is the sum of the primary helium core and the secondary helium core. The radius $R_{\rm rem}$ is calculated by Equation (81) in Hurley et al. (2000). The spin angular momentum is calculated by $kMR^{2}\Omega_{k}$, where $k$ is a parameter for the momentum of inertia which is calculated from the density profile. We assume that $k$ of the highly rotating helium star is the same as the value adopted in Hurley et al. (2000) for the helium dense convective core of giants ($k=0.21$). After the merger, the highly rotating helium star loses the angular momentum by the stellar wind mass loss. The stellar wind mass lose rate of Pop III stars is very weak ($\sim 10^{-14}{M}_{\odot}~{}\rm{yr}^{-1}$). But, in the case of highly rotating stars, the rotation effect enhances the mass loss rate. We use the following formula as the mass loss rate for the rotating helium stars,

[TABLE]

(Yoon et al., 2012), where $\tau_{\rm KH}=GM^{2}/RL$, $M$, $R$, $L$, $\Omega$, and $\Omega_{\rm K}$ are the Kelvin-Helmholtz timescale, the stellar mass, the stellar radius, the stellar luminosity, the angular velocity of the star, and the angular velocity of the Kepler rotation, respectively. The angular momentum loss due to the wind mass loss is written as

[TABLE]

We assume that the highly rotating helium stars evolve as the chemically homogeneous stars (Maeder, 1987) and will change into a CO star. In the case of the chemically homogeneous stellar evolution, the heavy elements are possibly carried to the stellar surface. But, we do not consider the effect of the surface heavy elements on the mass loss enhancement, because Krtička & Kubát (2009) and Muijres et al. (2012) show that such effect is rather moderate. When stars collapse, we treat them as a direct collapse. If the mass of the star after the collapse is larger than 3 $M_{\odot}$, the star is regarded as a BH. If the mass of helium stars is more massive than 60 ${M}_{\odot}$, the stars possibly cause pair instability supernovae (Woosley et al., 2007). Thus, we assume highly rotating helium stars whose masses are $3~{}{M}_{\odot}<M<60~{}{M}_{\odot}$ as GRB progenitors.

On the other hand, if a giant and a main sequence star merge during the CE phase, the merged remnant becomes a highly rotating main-sequence star. The mass of the highly rotating main-sequence star is the sum of the primary giant’s core and the secondary main-sequence star. We use the results of Yoon et al. (2012) to determine the fate of highly rotating main-sequence stars. Yoon et al. (2012) shows that the highly rotating main-sequence star whose mass is $13~{}{M}_{\odot}\lesssim M$ can evolve as the chemically homogeneous, and if their mass is $\lesssim 84~{}{M}_{\odot}$, they do not become a pair instability supernova, and the inner cores of those stars have a significantly high angular momentum to launch a GRB jet. Thus, we assume highly rotating main-sequence stars whose mass range is $13~{}{M}_{\odot}<M<84~{}{M}_{\odot}$ become long GRB progenitors (Yoon et al., 2012).

3.4. Classical long GRB rate from Pop III

Table 2 shows the numbers of the long GRB progenitors for $10^{6}$ binaries, obtained from the calculation shown in the section 3.3. A few percents of Pop III binaries can cause classical GRBs. We calculate the long GRB rate $R_{\rm GRB}$ from Pop III stars, using the beaming factor $f_{\rm B}=0.01$, the binary fraction $f_{\rm b}=0.5$, Pop III SFRs, and the long GRB fraction of Pop III $f_{\rm GRB}$ which consist of highly rotating Helium stars and highly rotating main-sequence stars as

[TABLE]

The beaming factor $f_{\rm B}=0.01$ is chosen to make the opening angle the same order of the typical value for the Pop I&II case (Liang et al., 2008). The binary fraction $f_{\rm b}=0.5$ is also the same as those in the Pop I&II case (Sana et al., 2012, 2013), which is consistent with the BH-BHS merger rate Kinugawa et al. (2016); Belczynski et al. (2016).

Figure 6 shows the classical long GRB rate from Pop III stars. Figure 7 shows the expectation of the classical GRB detection originated from Pop III stars based on the binary interaction model with an instrument with a sensitivity of $10^{-10}~{}\mbox{erg}~{}\mbox{cm}^{2}~{}\mbox{s}^{-1}$ in 0.5–4 keV, and a field of view 0.2 str. The same parameters as those in Lien et al. (2014) are adopted as the classical GRB luminosity function, and the modified Yonetoku relation is used for the spectral peak energy. The black and red lines are detection rates with SFR in de Souza et al. (2011) and Inayoshi et al. (2016), respectively. The upper and lower lines correspond to the parameter of $\alpha_{\rm CE}\lambda=0.1$ and $1$, respectively. The blue dashed line is the fiducial detection rate for Pop I&II GRBs based on the HPhi-1 model taken from Figure 2.

Although the GRB fraction is small compared to the assumption for ultra-long GRBs in section 3.2, the brighter luminosity function provides higher detection rates as shown in Figure 5. The SFRDs for Pop III stars assumed in this paper, which do not violate the constraints given by the Planck observation, imply higher detection rates for Pop III GRBs than the rate for Pop I&II GRBs. This is encouraging for the future observational missions such as HiZ-GUNDAM or THESEUS.

4. Conclusion and Discussion

The Hphi-1 model suggests that the SFRD calculated from clustering analysis of galaxies and UV luminosity function observations steeply decreases at high redshift compared to the extrapolated SFRD of the Madau & Dickinson (2014) model. This result shows that the Pop I&II stars are hard to contribute for long GRBs at high redshift. However, the SFRD of Pop III stars can be higher than Pop I&II SFRD. At $z=8-9$, the detectable Pop I&II GRB rate based on HPhi-1+Lien et al. (2014) model is 136 $\rm yr^{-1}$ in the whole sky. On the other hand, the Pop III GRB rate using Inayoshi et al. (2016) SFRD is 50 $\rm yr^{-1}$. At $z=9-10$, the Pop III GRB rate (63 $\rm yr^{-1}$) is almost the same as the Pop I&II GRB rate (51 $\rm yr^{-1}$). GRB events at $z=8-9$ like GRB090423 (Chandra et al., 2010) and GRB090429B (Cucchiara et al., 2011) would be Pop III GRBs with a probability of a few tens of percent. However, we have not found Pop III-like signature for GRB090423 and GRB090429B at present. Note that our binary merger model for Pop III GRBs yields classical GRBs, whose characteristic may be similar to other usual long GRBs.

In this paper, we consider the GRB from Pop III stars, using two SFRDs of Pop IIII considering the constraint from the Planck observation. We calculate the detection rate of Pop III GRBs by future observations such as HiZ-GUNDAM and THESEUS. In the pessimistic model, since the Pop III stars hold the hydrogen envelope because of the weak stellar wind, the Pop III stars are hard to launch a classical long GRB. In this case, Pop III stars might launch an ultra-long GRB. Only if more than 10% of Pop III stars launch a GRB jet, the future missions can detect such an ultra-long GRB per year.

However, many massive binary black holes confirmed by gravitational waves (The LIGO Scientific Collaboration & the Virgo Collaboration, 2018) might be remnants of Pop III binaries (Kinugawa et al., 2016). If a significant fraction of Pop III stars are formed as a binary, we expect that highly rotating helium stars and highly rotating main-sequence stars are formed via the binary merger, and they evolve as a chemically homogeneous star. Our population synthesis calculation shows that several % of Pop III binaries become such highly rotating stars which possibly launch a long GRB. On the other hand, Belczynski et al. (2007) consider Pop III GRB progenitors which lose the envelope and spin up due to tidal spin up and show that such progenitors may be a very small fraction ($\lesssim 1\%$). Thus, the binary merger is more effective process to make a Pop III GRB progenitor than the tidal spin up. If such highly rotating stars launch a classical GRB resembling low redshift long GRBs, they can be detected by HiZ-GUNDAM and THESEUS. Those future observations help us reveal the Pop III SFRD.

Acknowledgment

We appreciate D. Yonetoku and T. Ishiyama for the information on the instruments for HiZ-GUNDAM and for the data of N-body simulations, respectively. This work was supported by JSPS KAKENHI Grant Number 18J00558(TK), 16J03329(YH), 16K05291, and 18K03665 (KA). This work is carried out by the joint research program of the Institute for Cosmic Ray Research (ICRR), The University of Tokyo.