The Stellar-mass Function of Long Gamma-Ray Burst Host Galaxies

As the extremely powerful explosions of massive stars (Woosley & Bloom, 2006), long-duration gamma-ray bursts (GRBs) are detectable out to very high redshifts (Salvaterra et al., 2009b; Tanvir et al., 2009; Cucchiara et al., 2011). The redshift distribution of long GRBs and the properties of their host galaxies can therefore be used to probe the evolution of the star formation rate (SFR) density with cosmic time (e.g., Totani 1997; Wijers et al. 1998; Lamb & Reichart 2000; Porciani & Madau 2001). Relevant studies in this field have been carried out for more than two decades, but it remains an open question whether long GRBs can be good tracers of cosmic star formation. By studying the GRB redshift distribution, there is a general agreement about the fact that the comoving GRB rate does not purely follow the cosmic SFR, but exhibits an enhancement of GRB detections at high redshifts ($z\gtrsim 2-3$) with respect to the expectation from the star formation history (e.g., Daigne et al. 2006; Guetta & Piran 2007; Le & Dermer 2007; Kistler et al. 2008; Li 2008; Qin et al. 2010; Wanderman & Piran 2010; Virgili et al. 2011; Lu et al. 2012; Robertson & Ellis 2012; Wang 2013; Lien et al. 2014; Wei et al. 2014; Perley et al. 2015, 2016c; Vergani et al. 2015; Japelj et al. 2016; Wei & Wu 2017; Palmerio et al. 2019; Lan et al. 2019, 2021; Ghirlanda & Salvaterra 2022). Several possible mechanisms have been proposed to reconcile the discrepancy between the GRB rate and the SFR, such as cosmic metallicity evolution (Langer & Norman, 2006; Li, 2008; Salvaterra et al., 2012; Lloyd-Ronning et al., 2019; Ghirlanda & Salvaterra, 2022), an evolving initial mass function of stars (Xu & Wei, 2009; Wang & Dai, 2011), and an evolving luminosity function of GRBs (e.g., Lloyd-Ronning et al. 2002; Firmani et al. 2004; Salvaterra & Chincarini 2007; Salvaterra et al. 2009a, 2012; Cao et al. 2011; Tan & Wang 2015; Pescalli et al. 2016; Lan et al. 2019, 2021).

Among the possible mechanisms, cosmic metallicity evolution appears to be a compelling one. The progenitors of long GRBs are expected to be preferentially localized in low-metallicity environments, where the loss of mass and angular momentum are lessened, allowing a relativistic jet to form (MacFadyen & Woosley, 1999; Hirschi et al., 2005; Yoon & Langer, 2005; Woosley & Heger, 2006; Yoon et al., 2006). And indeed it has been shown that long GRBs tend to occur in lower metallicity (as well as lower stellar-mass) galaxies, with metallicity $Z\leq 0.3Z_{\odot}$ (see Perley et al. 2016a and references therein). However, some works suggested that GRBs could occur in high-metallicity host environments (e.g., Levesque et al. 2010; Elliott et al. 2012, 2013; Savaglio et al. 2012; Hao & Yuan 2013). It is important to note that previous studies of the GRB redshift distribution are necessarily imprecise. Specifically, their derived GRB rate at any redshift is an average over all galaxies at that epoch (a diverse population covering orders of magnitude in metallicities, SFRs, etc.) and cannot disclose which types of galaxies contribute most to the GRB rate. To thoroughly understand the link between GRB production and environment, one also is required to characterize the population of GRB host galaxies directly (Perley et al., 2016c).

One of the difficulties in attempting to connect GRB observations and their environments to particular progenitor systems is obtaining a complete sample of GRB host galaxies. Recently, Perley et al. (2016b, c) carried out a new survey of the GRB host galaxy population. Their survey provides the first host galaxy sample that is unbiased and sufficiently large to statistically examine redshift evolution in the host galaxy population in detail. In this work, we make use of the complete sample presented in Perley et al. (2016b, c) to quantify the stellar-mass function (SMF) of GRB host galaxies. Meanwhile, we use the observed redshift distribution of the complete sample to probe and constrain the evolution of the GRB host population in redshift.

The structure of the paper is as follows: In Section 2, we introduce the GRB host sample used for our analysis. We describe our analysis method in Section 3, while the model results are presented in Section 4. Finally, in Section 5 we draw our conclusions. Throughout this work we adopt a flat $\Lambda$CDM model with cosmological parameters $H_{0}=70$ km s^-1 Mpc^-1, $\Omega_{\rm m}=0.3$, and $\Omega_{\Lambda}=0.7$.

The redshifts of about 40% of all GRBs detected by the Swift satellite have been measured (Lan et al., 2021). Though this is a significant improvement compared to the pre-Swift era, from the point of view of characterizing the GRB host population or redshift distribution, the sample is still far from being considered unbiased and complete. The Swift GRB Host Galaxy Legacy Survey (SHOALS) built a new unbiased sample of long GRB host galaxies, which is composed of 119 GRBs at $0.03<z<6.29$ drawn from the Swift GRB catalog using a series of selection criteria and requiring the Swift/Burst Alert Telescope (BAT) fluence to exceed $S>10^{-6}$ erg $\rm cm^{-2}$ (Perley et al., 2016b). Their optimized observability cuts combined with host galaxy follow-up enabled redshift measurements for 110 of these bursts (representing a completeness of 92%). Subsequently, Perley et al. (2016c) presented rest-frame near-infrared (NIR) luminosities and stellar masses for this uniform sample of 119 GRB host galaxies through the Spitzer Infrared Array Camera (IRAC) observations. As the first large, highly complete, and multiband survey of an unbiased GRB host galaxy sample, SHOALS provides unprecedented insight into galaxy evolution and cosmic history (Perley et al., 2016b, c).

In this work, we use the SHOALS sample to investigate the stellar-mass distribution of GRB hosts across cosmic history and its possible implications for interpreting the observed GRB redshift distribution. We calculate the isotropic-equivalent energies $E_{\rm iso}$ in the 1–$10^{4}$ keV rest-frame energy range for these 110 GRBs having redshifts. Note that there is a low-luminosity/energy burst (GRB 060218) with $E_{\rm iso}\leq 10^{49}$ erg. As the low-luminosity/energy bursts have been claimed to be a distinct population (Soderberg et al., 2004; Cobb et al., 2006; Chapman et al., 2007; Liang et al., 2007), GRB 060218 is discarded in our analysis. The remaining sample includes 109 GRBs, whose stellar mass–redshift distribution is presented in the left panel of Figure 1. It is also noteworthy that 30 out of these 109 GRB host galaxies only have stellar-mass upper limits (Perley et al., 2016c). In principle, the galaxies with stellar-mass upper limits should be removed in the determination of the stellar-mass distribution of GRB hosts. However, the direct abandonment of these hosts would introduce a new bias on shaping the GRB redshift distribution. In the following, we will simply treat the upper mass limits as the stellar-mass measurements for these 30 GRB hosts.

In order to infer the model-free parameters, a maximum likelihood method firstly proposed by Marshall et al. (1983) is adopted. The likelihood function $\mathcal{L}$ is defined as (e.g., Chiang & Mukherjee 1998; Narumoto & Totani 2006; Ajello et al. 2009, 2012; Abdo et al. 2010; Zeng et al. 2014, 2016; Lan et al. 2019, 2021)

where $N_{\rm obs}$ is the number of the observed sample, $N_{\rm exp}$ is the expected number of GRB detections, and $\Phi(E_{\rm iso},\;M_{\ast},\;z,\;t)$ is the observed rate of GRBs per unit time at redshift $z\sim z+\mathrm{d}z$ with energy $E_{\rm iso}\sim E_{\rm iso}+\mathrm{d}E_{\rm iso}$ and stellar mass $M_{\ast}\sim M_{\ast}+\mathrm{d}M_{\ast}$, which can be expressed as

where $\Delta\Omega=1.33$ sr is the $Swift$/BAT field of view, $\psi(z)$ is the comoving formation rate of GRBs (in units of $\rm Mpc^{-3}$ $\rm yr^{-1}$), $(1+z)^{-1}$ accounts for the cosmological time dilation, $\phi(E_{\rm iso})$ is the normalized isotropic-equivalent energy function of GRBs, and $\varphi(M_{\ast},z)$ is the normalized SMF of GRB hosts, which evolves with redshift dependents on the assumed model (see more below). ${\rm d}V(z)/{\rm d}z=4\pi cD_{L}^{2}(z)/[H_{0}(1+z)^{2}\sqrt{\Omega_{\rm m}(1+z)^{3}+\Omega_{\Lambda}}]$ is the comoving volume element in a flat $\Lambda$CDM model, where $D_{L}(z)$ is the luminosity distance.

Since in the collapsar model the births of GRBs just mean the deaths of short-lived massive stars, the GRB formation rate $\psi(z)$ should in principle be related to the cosmic SFR $\psi_{\star}(z)$, i.e.,

where $\eta$ is the GRB formation efficiency in units of ${\rm M}_{\odot}^{-1}$. The SFR (in units of ${\rm M}_{\odot}$ $\rm yr^{-1}$ $\rm Mpc^{-3}$) can be parameterized as (Hopkins & Beacom, 2006; Li, 2008)

For the GRB energy function $\phi(E_{\rm iso})$, we employ a broken power-law form

where $A$ is a normalization factor, and $a$ and $b$ are the power-law indices before and after the break energy $E_{b}$.

Empirically, the galaxy SMF can be well described by a Schechter form (Drory & Alvarez, 2008). Therefore, we use the same Schechter expression for the stellar-mass distribution of GRB host galaxies:

where $\varphi_{0}$ is a normalization factor, $\xi$ is the power-law index, and $M_{b,0}$ is the break mass.

Considering the fluence threshold that is used to define the SHOALS sample (i.e., $S_{\rm lim}=10^{-6}$ erg $\rm cm^{-2}$ in the 15–150 keV), the expected number of GRB detections should be

where $T\sim 7$ yr is the observational period of Swift that covers the SHOALS sample and $M_{\ast,\rm lim}(z)$ represents the stellar-mass threshold of Spitzer-IRAC at redshift $z$. Since $z<10$ and $M_{\ast}<10^{12}$ ${\rm M}_{\odot}$ for the current GRB host population, we adopt a maximum redshift $z_{\rm max}=10$ and a maximum stellar mass $M_{\ast,\rm max}=10^{12}$ ${\rm M}_{\odot}$. The energy function is assumed to extend between minimum and maximum energies $E_{\rm min}=10^{51}$ erg and $E_{\rm max}=10^{56}$ erg. The energy threshold appearing in Equation (7) can be computed by

where $N(E)$ is the GRB photon spectrum. To describe the burst spectrum, we employ a typical Band function with low- and high-energy spectral indices $-1$ and $-2.3$, respectively (Band et al., 1993; Preece et al., 2000; Kaneko et al., 2006). For a given isotropic-equivalent energy $E_{\rm iso}$, the peak energy of the spectrum $E_{p}$ is estimated through the empirical $E_{p}$–$E_{\rm iso}$ correlation (Amati et al., 2002; Nava et al., 2012), i.e., $\log\left[E_{p}(1+z)\right]=-29.60+0.61\log E_{\rm iso}$.

For each GRB host galaxy (with measured redshift) from the SHOALS sample, the absolute magnitude at a wavelength of $3.6/(1+z)$ $\mu$m can be calculated by

where $m_{\rm obs}$ is the observed IRAC aperture magnitude and $\mu(z)$ is the distance modulus at $z$. It has been suggested that the quantity $M_{3.6/(1+z)}$ alone can be used as a reasonable stellar-mass proxy (Perley et al., 2016c). As shown in the right panel of Figure 1, there is a strong correlation between the absolute magnitudes $M_{3.6/(1+z)}$ and the logarithmic stellar masses $\log M_{\ast}$ of GRB host galaxies. For simplicity, here we use a second-order polynomial to model the $\log M_{\ast}$–$M_{3.6/(1+z)}$ relation

where the best-fitting parameters are $p_{0}=3.59_{-1.13}^{+1.13}$ , $p_{1}=-0.11_{-0.11}^{+0.11}$, and $p_{2}=0.0087_{-0.0025}^{+0.0025}$ at 68% confidence level. The observed magnitude $m_{\rm obs}$ can then be translated to stellar-masses $M_{\ast}$ at different $z$ via Equations (9) and (10). We show that the stellar-mass threshold $M_{\ast,\rm lim}$ appearing in Equation (7) may be approximated well using a magnitude limit of Spitzer-IRAC, $m_{\rm obs,lim}=25.5$ mag, according to Equations (9) and (10) (see the dotted line in the left panel of Figure 1).

Considering that (i) GRBs may not be good tracers of star formation, (ii) the GRB energy function may evolve with redshift, and (iii) there is a strong evolution in the galaxy SMF (Fontana et al., 2006; Drory & Alvarez, 2008; Duncan et al., 2014), here we explore four different evolution scenarios. In the first no evolution model, the GRB formation rate strictly follows the cosmic SFR, and neither the GRB energy function nor the GRB host SMF evolves with redshift. In the energy evolution model, while the GRB rate is still proportional to the SFR and the GRB host SMF is still taken to be nonevolving, the break energy in the GRB energy function increases with redshift as $E_{b,z}=E_{b,0}(1+z)^{\delta}$. In the density evolution model, the GRB rate follows the SFR in conjunction with an extra evolving factor, i.e., $\psi(z)=\eta\psi_{\star}(z)(1+z)^{\delta}$, and both the GRB energy function and the GRB host SMF are taken to be nonevolving. In the last mass evolution model, while the GRB rate is still fully consistent with the SFR and the break energy in the GRB energy function is still a constant, the GRB host SMF evolves with redshift as the following relations (Drory & Alvarez, 2008):

where $\delta$ and $\gamma$ are two additional evolution parameters. When $\delta\rightarrow 0$ and $\gamma\rightarrow 0$, Equation (3) reduces to Equation (6).

We optimize the model-free parameters by maximizing the likelihood function (Equation (1)). Since there are several parameters in our models, the Markov Chain Monte Carlo (MCMC) technology is applied to give multidimensional parameter constraints from the observational data. Here we employ the MCMC code from GWMCMC,¹¹1https://github.com/grinsted/gwmcmc which is an implementation of the Goodman & Weare 2010 Affine invariant ensemble MCMC sampler (Goodman & Weare, 2010; Foreman-Mackey et al., 2013).

Using the above analysis method, we can infer the values of each model’s free parameters, including the GRB formation efficiency $\eta$, the GRB energy function, the host galaxy SMF, and the evolution parameters. The best-fitting parameter values together with their $1\sigma$ uncertainties for different models are listed in Table 1. The last two columns report the natural-logarithmic-likelihood value $\ln\mathcal{L}$ and the Akaike information criterion (AIC) score, respectively. The AIC score of each fitted model is evaluated as (Akaike, 1974; Liddle, 2007)

where $f$ is the number of model-free parameters. If there are two or more models for the same data, $\mathcal{M}_{\rm 1},\mathcal{M}_{\rm 2},...,\mathcal{M}_{\rm N}$, and they have been separately fitted, the one with the smallest AIC score is the one most preferred by this criterion. With the${\rm AIC}_{i}$ characterizing model $\mathcal{M}_{i}$, the unnormalized confidence that this model is true is the Akaike weight $\exp(-{\rm AIC}_{i}/2)$. The relative probability of $\mathcal{M}_{i}$ being the correct model is then

Whether long GRBs are biased tracers of the star formation activity remains an open question. Studies of the GRB redshift distribution seem to indicate that the comoving GRB rate density exhibits similar behavior to the star formation history at low redshifts ($z\lesssim 2$), but shows an enhancement of GRBs at high redshifts ($z\gtrsim 2-3$) compared to what would be expected if the GRB rate followed the cosmic SFR strictly (e.g., Daigne et al. 2006; Kistler et al. 2008; Perley et al. 2015). However, in the technique of star formation models, the GRB rate measured at any epoch is an average over all galaxies at that time (a diverse galaxy population spanning orders of magnitude in metallicities, SFRs, etc.). Thus, it is difficult to inspect which types of galaxies contribute most to the GRB formation rate. To get a better understanding of the connection between GRB production and environment, one has to characterize the population of GRB host galaxies directly (Perley et al., 2016c).

Based on the complete GRB host sample presented in Perley et al. (2016b, c), we construct the SMF of GRB host galaxies as well as the isotropic-equivalent energy function and redshift distribution of GRBs in the frameworks of different evolution models. The maximum likelihood method is applied to perform a joint analysis for the distributions of redshift, isotropic-equivalent energy, and host galaxy mass of the GRB sample, and the MCMC algorithm is adopted to provide the best-fitting parameters in each model. Our results confirm that GRBs must have experienced some kind of redshift evolution, being more luminous or more population in the past than present day. In order to match the observed distributions, the typical GRB energy should increase to $(1+z)^{2.47^{+0.73}_{-0.89}}$ or the GRB rate density should rise to $(1+z)^{1.82^{+0.22}_{-0.59}}$ on top of the known cosmic evolution of the SFR. We find that the GRB host SMF can be well described by the Schechter function with a power-law index $\xi\approx-1.10$ and a break mass $M_{b,0}\approx 4.9\times 10^{10}$ ${\rm M}_{\odot}$, independent of the assumed evolutionary effects.

However, since the SMF of all galaxies is suggested to be evolving with redshift, the intrinsic form of the stellar-mass distribution of GRB host galaxies may also have a redshift dependence. We also explore models in which the GRB host SMF evolves with redshift. In this case, we find that the redshift distribution of GRBs with measured redshift, energy, and stellar mass can be successfully fitted by the star formation history, which involves with an evolving SMF of GRB hosts. In order to reproduce the observed redshift distribution well, the normalization coefficient of the SMF should increase to $(1+z)^{1.96^{+0.44}_{-0.47}}$ and the break mass should evolve with redshift according to $\log M_{b}(z)=11.24^{+0.25}_{-0.43}-0.38^{+0.28}_{-0.19}\ln(1+z)$. Our results suggest that long GRBs may be biased tracers of star formation and the evolving SMF may also explain the discrepancy between the GRB rate and the SFR.

The current sample size of GRB hosts is still too small to exhibit an obvious redshift evolution in the stellar mass distribution. Ongoing large-area surveys will be able to increase the host galaxy sample size. With more abundant observational information in the future, the properties of the population of GRB host galaxies and their possible evolution with redshift will be better characterized, and determinations on the SMF of GRB hosts, as discussed in this work, will be considerably improved.