Constraints on warm dark matter models from high-redshift long gamma-ray bursts

Abstract

Structures in warm dark matter (WDM) models are exponentially suppressed below a certain scale, characterized by the dark matter particle mass, m_x. Since structures form hierarchically, the presence of collapsed objects at high redshifts can set strong lower limits on m_x. We place robust constraints on m_x using recent results from the Swift data base of high-redshift gamma-ray bursts (GRBs). We parametrize the redshift evolution of the ratio between the cosmic GRB rate and star formation rate (SFR) as ∝(1 + z)^α, thereby allowing astrophysical uncertainties to partially mimic the cosmological suppression of structures in WDM models. Using a maximum-likelihood estimator on two different z > 4 GRB subsamples (including two bursts at z > 8), we constrain m_x ≳ 1.6–1.8 keV at 95 per cent CL, when marginalized over a flat prior in α. We further estimate that 5 years of a Sino-French space-based multi-band astronomical variable objects monitor like mission would tighten these constraints to m_x ≳ 2.3 keV. Our results show that GRBs are a powerful probe of high-redshift structures, providing robust and competitive constraints on m_x.

methods: statistical, gamma-ray burst: general, dark matter

INTRODUCTION

The current concordance cosmology, in which structure formation proceeds in a hierarchal manner driven by pressureless cold dark matter (CDM), has been remarkably successful in explaining the observed properties of large-scale structures in the Universe (e.g. Tegmark et al. 2006; Benson 2010) and the cosmic microwave background (e.g. Komatsu et al. 2011). Such observables probe scales in the range ∼1 Gpc down to ∼10 Mpc. On smaller scales, ≲1 Mpc, there are still some discrepancies between standard ΛCDM and observations (e.g. Menci, Fiore & Lamastra 2012). For instance, N-body simulations predict more satellite galaxies than are observed both around our galaxy (the so-called ‘missing satellite problem’; e.g. Klypin et al. 1999; Moore et al. 1999), and in the field as recently noted by the Arecibo Legacy Fast ALFA (ALFALFA) survey (e.g. Papastergis et al. 2011; Ferrero et al. 2012). Furthermore, simulations of the most massive Galactic CDM subhaloes are too centrally condensed to be consistent with the kinematic data of the bright Milky Way satellites (e.g. Boylan-Kolchin, Bullock & Kaplinghat 2011). Moreover, observations of small galaxies show that their central density profile is shallower than that predicted by CDM N-body simulations (e.g. Moore 1994; de Blok et al. 2001; Donato et al. 2009; Governato et al. 2012; Macciò et al. 2012).

Baryonic feedback is a popular prescription for resolving such discrepancies. Feedback caused by supernova explosions and heating due to the UV background (UVB) may suppress the baryonic content of low-mass haloes (e.g. Governato et al. 2007; Mashchenko, Wadsley & Couchman 2008; Busha et al. 2010; Sobacchi & Mesinger 2013b) and make their inner density profile shallower (e.g. de Souza & Ishida 2010; de Souza et al. 2011a). However, accurately matching observations is still difficult even when tuning feedback recipes (e.g. Boylan-Kolchin, Bullock & Kaplinghat 2012).

An alternative explanation would be found if dark matter (DM) consisted of lower mass (∼keV) particles, and thus was ‘warm’ (WDM; e.g. Bode, Ostriker & Turok 2001; Khlopov & Kouvaris 2008; de Vega & Sanchez 2012; de Vega, Salucci & Sanchez 2012; Destri, de Vega & Sanchez 2013; Kamada et al. 2013; Kang, Macciò & Dutton 2013). The resulting effective pressure and free streaming would decrease structure on small scales, though again fine tuning might be required to fully match all the observations (e.g. Boylan-Kolchin et al. 2011; Borriello et al. 2012; Macciò et al. 2012).

The most powerful testbed for these scenarios is the high-redshift Universe. Structure formation in WDM models (or in any cosmological model with an equivalent power-spectrum cutoff) is exponentially suppressed on small scales (e.g. Schneider et al. 2012; Schneider, Smith & Reed 2013). Since structures form hierarchically, these small haloes are expected to host the first galaxies. If indeed dark matter were sufficiently ‘warm’, the high-redshift Universe would be empty. Therefore, the mere presence of a galaxy at high redshift can set strong lower limits on the WDM particle mass.

Due to their high luminosity, gamma-ray bursts (GRBs) constitute a remarkable tool to probe the high-z Universe and small-scale structures. They provide a glimpse of the first generations of stars (e.g. de Souza, Yoshida & Ioka 2011b; de Souza et al. 2012), as well as constraints on primordial non-Gaussianity (Maio et al. 2012). As pointed out by Mesinger, Perna & Haiman (2005), the detection of a single GRB at z > 10 would provide very strong constraints on WDM models.

Here, we extend the work of Mesinger et al. (2005) by presenting robust lower limits on WDM particle masses, using the latest Swift GRB data. The current data, including many redshift measurements, allow us to perform an improved statistical analysis by directly comparing the distribution of bursts in various models as a function of redshift. Furthermore, we make more conservative¹ assumptions throughout the analysis, such as normalizing the star formation rate (SFR)-to-GRB ratio at high redshifts (thereby using a shorter, more accurate lever arm which minimizes modelling uncertainty), using an unbiased luminosity function (LF) and allowing the SFR-to-GRB ratio to evolve with redshift. Finally, we study the effectiveness of future observations in improving the current constraints.

Current limits on DM masses, m_x, are motivated by several observations. The Lyman α forest implies m_x ≳ 1 keV (e.g. Viel et al. 2008) and m_νs > 8 keV for sterile neutrinos (Seljak et al. 2006; Boyarsky et al. 2009a). Likewise, WDM models with a too warm candidate (m_x < 0.75 keV) cannot simultaneously reproduce the stellar mass function and the Tully–Fisher relation (Kang et al. 2013). Also, the fact that reionization occurred at z ≳ 6 implies m_x ≳ 0.75 keV (Barkana, Haiman & Ostriker 2001). However, all of these limits are strongly affected by a degeneracy between astrophysical (i.e. baryonic) processes and the DM mass. Our approach in this work is more robust, driven only by the shape of the redshift evolution of the z > 4 SFR. Furthermore, it is important to note that the SFR is exponentially attenuated at high redshifts in WDM models. Since astrophysical uncertainties are unable to mimic such a rapid suppression, probes at high redshifts (such as GRBs and reionization) are powerful in constraining WDM cosmologies.

The outline of this paper is as follows. In Section 2, we discuss how we derive the DM halo mass function and SFR in WDM and CDM models. In Section 3, we derive the corresponding GRB redshift distribution. In Section 4, we discuss the adopted observed GRB sample. In Section 5, we present our analysis and main results. In Section 6, we discuss possible future constraints using a theoretical mock sample. Finally, in Section 7, we present our conclusions.²

STRUCTURE FORMATION IN A WDM-DOMINATED UNIVERSE

Massive neutrinos from the standard model (SM) of particle physics were one of the first DM candidates. However, structures formed in this paradigm are incompatible with observations. Other alternative DM candidates usually imply an extension of the SM. The DM particle candidates span several orders of magnitude in mass (Boyarsky, Ruchayskiy & Shaposhnikov 2009b): axions with a mass of ∼10⁻⁶ eV, first introduced to solve the problem of charge parity (CP) violation in particle physics, supersymmetric particles (gravitinos, neutralinos and axinos) with mass in the range ∼eV–GeV, superheavy DM, also called Wimpzillas [also considered as a possible solution to the problem of cosmic rays observed above the Greisen–Zatsepin–Kuzmin (GZK) cutoff], Q-balls and sterile neutrinos with mass in the ∼ keV range, just to cite a few. For a review about DM candidates, see Bertone, Hooper & Silk (2005). Two promising candidates for WDM are the sterile neutrino (Dodelson & Widrow 1994; Shaposhnikov & Tkachev 2006) and gravitino (Ellis et al. 1984; Moroi, Murayama & Yamaguchi 1993; Kawasaki, Sugiyama & Yanagida 1997; Primack 2003; Gorbunov, Khmelnitsky & Rubakov 2008).

In WDM models, the growth of density perturbations is suppressed on scales smaller than the free-streaming length. The lighter the WDM particle, the larger the scale below which the power spectrum is suppressed. In addition to this power-spectrum cutoff, one must also consider the residual particle velocities. As described in Barkana et al. (2001), these act as an effective pressure, slowing the early growth of perturbations. Below we describe how we include these two effects in our analysis (for more information, see Barkana et al. 2001; Mesinger et al. 2005).

Power-spectrum cutoff

The free-streaming scale, below which the linear perturbation amplitude is suppressed (e.g. Colombi, Dodelson & Widrow 1996; Bode et al. 2001; Viel et al. 2005), is given by the comoving scale

\begin{equation} R_{\rm fs} \approx 0.11\left(\frac{\Omega _{\rm x}h^2}{0.15}\right)^{1/3}\Big(\frac{m_{\rm x}}{\rm keV} \Big)^{-4/3} \rm Mpc, \end{equation}

(1)

where Ω_x is the total energy density contained in WDM particles relative to the critical density, h is the Hubble constant in units of 100 km s⁻¹ Mpc⁻¹ and m_x is the WDM particle mass.

The resulting modification of the matter power spectrum can be computed by multiplying the CDM power spectrum P_CDM(k) by an additional transfer function (Bode et al. 2001):

\begin{equation} P_{\rm WDM}(k) = P_{\rm CDM}(k)\left[1+(\epsilon k)^{2\mu }\right]^{-5\mu }, \end{equation}

(2)

where μ = 1.12 and

\begin{equation} \epsilon = 0.049 \left(\frac{\Omega _{\rm x}}{0.25}\right)^{0.11}\Big(\frac{m_{\rm x}}{\rm keV}\Big)^{-1.11}\left(\frac{h}{0.7}\right)^{1.22}\,h^{-1}\,{\rm Mpc}. \end{equation}

(3)

Effective pressure

Structure formation in WDM models will be further suppressed by the residual velocity dispersion of the WDM particles, which delays the growth of perturbations. As shown in Barkana et al. (2001) and Mesinger et al. (2005), this ‘effective pressure’ is of comparable importance to the power-spectrum cutoff in determining the WDM mass functions. The pressure can be incorporated in the halo mass function by raising the critical linear extrapolated overdensity threshold at collapse, δ_c(M, z). Using spherically symmetric hydrodynamics simulations and exploiting the analogy between the WDM effective pressure and gas pressure, Barkana et al. (2001) computed the modified δ_c(M, z). They showed that one can define an effective WDM Jeans mass, i.e. the mass scale below which collapse is significantly delayed by the pressure. We will denote this scale as

\begin{equation} M_{\rm WDM} \approx 1.8\times 10^{10}\left(\frac{\Omega _{\rm x}h^2}{0.15}\right)^{1/2}\left(\frac{m_{\rm x}}{1\,{\rm keV}}\right)^{-4} , \end{equation}

(4)

where we make the standard assumption of a fermionic spin-1/2 particle. Note that equation (4) differs from the original one proposed by Barkana et al. (2001) by a factor of 60. With this adjustment, we find that the collapsed fractions, F_coll(z), simulated assuming a sharp mass cutoff at M_WDM are in good agreement with the full random-walk procedure [assuming a more gradual rise in δ_c(M)] of Barkana et al. (2001) (see Fig. 2). The former approach has the advantage of being considerably faster and simpler.

Collapsed fraction of haloes and cosmic star formation

We use the Sheth–Tormen mass function, f_ST (Sheth & Tormen 1999), to estimate the number density of DM haloes, n_ST(M, z), with mass greater than M:

\begin{equation} f_{\rm ST} = A\sqrt{\frac{2a_1}{\pi }} \left[1+\left(\frac{\sigma ^2}{a_1\delta _{{\rm c}}^2}\right)^p\right] \frac{\delta _{{\rm c}}}{\sigma }\exp {\left[-\frac{a_1\delta _{{\rm c}}^2}{2\sigma ^2}\right]}, \end{equation}

(5)

where A = 0.3222, a₁ = 0.707, p = 0.3 and δ_c = 1.686. The mass function f_ST can be related to n_ST(M, z) by

\begin{equation} f_{\rm ST} = \frac{M}{\rho _{\rm m}}\frac{\mathrm{d}n_{\rm ST}(M,z)}{\mathrm{d}\ln {\sigma ^{-1}}}, \end{equation}

(6)

where ρ_m is the total mass density of the background Universe. The variance of the linearly extrapolated density field σ(M, z) is given by

\begin{equation} \sigma ^{2}(M,z) = \frac{b^2(z)}{2\pi ^2}\int _{0}^{\infty } k^2P(k)W^2(k,M)\,\mathrm{d}k, \end{equation}

(7)

where b(z) is the growth factor of linear perturbations normalized to b = 1 at the present epoch and W(k, M) is a real space top-hat filter. In order to calculate the CDM power spectrum P_CDM(k), we use the camb code.³ For WDM models, we compute the power spectrum, P_WDM(k), with equation (2).

The fraction of mass inside collapsed haloes, F_coll(>M_min, z), is then given by

\begin{equation} F_{\rm {coll}}(z)=\frac{1}{\rho _{\rm {m}}}\int _{M_{\rm min}}^{\infty } \mathrm{d}MMn_{\rm {ST}}(M,z), \end{equation}

(8)

and the minimum mass is estimated as

\begin{equation} M_{\rm min} = {\rm Max}[M_{\rm gal}(z), M_{{\rm WDM}}(m_{\rm x})], \end{equation}

(9)

where M_gal corresponds to the minimum halo mass capable of hosting star-forming galaxies. We use |$M_{\rm gal} \equiv \bar{M}_{{\rm min}}$| from equation 11 in Sobacchi & Mesinger (2013a), who present a physically motivated expression for the evolution of M_gal which takes into account a gas cooling criterion as well as radiative feedback from an ionizing UVB during inhomogeneous reionization.⁴ In other words, M_gal is set by astrophysics, whereas M_WDM(m_x) is set by the particle properties of DM. In Fig. 1, we show these two limits as a function of redshift. As we shall see, our conclusions are derived from the regime where M_WDM > M_gal, and hence they are not sensitive to the exact value of M_gal. Another point worth highlighting is that models with m_x ≳ 3.5 keV have M_WDM(m_x) ≥ M_gal up to z ≈ 11. Thus, any constraint beyond ∼3.5 keV will no longer be robust, since galaxy formation can be suppressed even in CDM, below halo masses comparable to M_WDM. Therefore, observations at much higher redshifts are required to obtain tighter constraints.

Figure 1.

The minimum halo masses capable of hosting star-forming galaxies. The solid black line corresponds to the astrophysical limit, M_gal, from Sobacchi & Mesinger (2013a). The horizontal lines correspond to the cosmological cutoffs in WDM models, M_WDM (m_x = 0.5, 1, 1.5, 2, 2.5, 3 and 3.5 keV from top to bottom, respectively).

In Fig. 2, we plot the fraction of the total mass collapsed into haloes of mass >M_min, F_coll(>M_min, z). The shaded region shows the collapsed fraction in the CDM, with a range of low-mass cutoffs corresponding to virial temperatures 300 K <T_vir < 10⁴ K. The other curves correspond to WDM particle masses of m_x = 3.0, 2.5, 2.0, 1.5 and 1.0 keV (top to bottom). This figure is analogous to fig. 2 in Mesinger et al. (2005), serving to motivate equation (4). The fractions F_coll(>M_min, z), computed according to the full random-walk procedure used in Mesinger et al. (2005), are shown as solid black curves, while the approximation of a sharp cutoff at M_WDM (equation 4) corresponds to the red dashed curves.

Figure 2.

Fraction of the total mass collapsed into haloes of mass >M_min as a function of redshift, F_coll(>M_min, z). The shaded region shows the collapsed fraction in CDM, with a range of low-mass cutoffs corresponding to virial temperatures 300 K <T_vir < 10⁴ K. The other curves correspond to WDM particle masses of m_x = 3.0, 2.5, 2.0, 1.5 and 1.0 keV (top to bottom). The figure is an adapted version of fig. 2 from Mesinger et al. (2005) (computed using their cosmology), serving to motivate equation (4). Values of F_coll(>M_min, z) computed according to the full random-walk procedure used in Mesinger et al. (2005) are shown as solid black curves, while the approximation of a sharp cutoff at M_WDM used in this work corresponds to the red dashed curves.

The comoving SFR as a function of redshift is assumed to be proportional to dF_coll/dt:

\begin{equation} {\rm SFR}(z) \propto \frac{{\rm d}F_{\rm {coll}}}{{\rm d}t}. \end{equation}

(10)

The proportionality constant is irrelevant for our analysis, as it is subsumed in our normalization procedure below.

THEORETICAL REDSHIFT DISTRIBUTION OF GRBs

Under the hypothesis that the formation rate of long GRBs (LGRBs; duration longer than 2 s) follows the SFR (e.g. Totani 1997; Bromm & Loeb 2006; Campisi, Li & Jakobsson 2010; de Souza et al. 2011b), the comoving rate of GRBs, Ψ_GRB, can be expressed as

\begin{equation} \Psi _{\rm GRB}(z) =\zeta _{0}(1+z)^{\alpha } {\rm SFR}(z). \end{equation}

(11)

Here, ζ₀ is a constant that includes the absolute conversion from the SFR to the GRB rate. The evolutionary trend described by α may arise from several mechanisms (e.g. Kistler et al. 2009), with a possible explanation provided by the GRB preference for low-metallicity environments.⁵ A metallicity threshold seems to provide a natural explanation for the observed value of α ≈ 0.5–1 at low redshifts (z ≤ 4; e.g. Robertson & Ellis 2012). Such a metallicity threshold increases with redshift, whereas the characteristic halo mass decreases with redshift. In such a scenario, a value of α = 0 would be appropriate for our high-redshift (z > 4) analysis. Nevertheless, we conservatively keep α as a free parameter.

The observed number of GRBs in the range z₁ ≤ z ≤ z₂, N(z₁, z₂), can be expressed by

\begin{equation} N(z_1,z_2) = K\int _{z_1}^{z_2} \frac{\mathrm{d}N_{\rm GRB}}{\mathrm{d}z^{\prime }}\mathrm{d}z^{\prime }, \end{equation}

(12)

with

\begin{equation} \frac{\mathrm{d}N_{\rm GRB}}{\mathrm{d}z} = \Psi _{\rm GRB}(z)\frac{\Delta t}{1+z}\frac{\mathrm{d}V}{\mathrm{d}z}\mathcal {I}(z), \end{equation}

(13)

where the parameter K accounts for the efficiency of finding GRBs and measuring their redshift (e.g. area coverage, the survey flux limit, beaming factor of GRBs, etc.),⁶ dV/dz is the comoving volume element per unit redshift, Δt is the time interval in the observer rest frame⁷ and |$\mathcal {I}(z)$| is the integral over the GRB LF, p(L),

\begin{equation} \mathcal {I}(z)=\int _{\log {L_{{\rm lim}}(z)} }^{\infty }p(L)\,\mathrm{d}\log {L}. \end{equation}

(14)

To remove the dependence on proportionality constants, we construct cumulative distribution functions (CDFs) of GRBs over the redshift range z_i < z < z_max:

\begin{equation} N (< \!\!z|z_{\rm max}) = \frac{N(z_i,z)}{N(z_i,z_{\rm max})}. \end{equation}

(15)

The expected number of observed GRBs in a given redshift interval, for each combination |$\pmb {\theta }\equiv \lbrace m_{\rm x},\alpha \rbrace$|⁠, can be written as

\begin{equation} N(z_1,z_2;\pmb {\theta }) = \zeta _{0;\pmb {\theta }} \int _{z_1}^{z_2}(1+z)^{\alpha } {\rm SFR}(z;m_{\rm x})\frac{\Delta t}{1+z}\frac{\mathrm{d}V}{\mathrm{d}z}\mathcal {I}(z), \end{equation}

(16)

where |$\zeta _{0;\pmb {\theta }}$| is normalized so that each model recovers the observed GRB rate, N_obs(z₁, z₂), at z ≈ 4,

\begin{equation} \zeta _{0;\pmb {\theta }} = \frac{N_{\rm obs}(3,4)}{\int _{3}^{4}(1+z)^{\alpha } {\rm SFR}(z;m_{\rm x})\frac{\Delta t}{1+z}\frac{\mathrm{d}V}{\mathrm{d}z}\mathcal {I}(z)}. \end{equation}

(17)

N_obs(3, 4) is equal to 24 (7) for S1 (S2) samples, respectively. Note that normalizing at a relatively high z (z = 3–4) is indeed a conservative choice. If we had normalized at lower redshifts, say z ∼ 1–2, the absolute number of GRBs predicted by m_x = 0.5 keV and CDM, between z ∼ 3 and 4, would already diverge by a factor of ∼2 for α values in the range 0-2. In addition to being conservative, normalizing at high redshifts allows us to have a relatively short lever arm over which our simple scaling, SFR ∝ dF_coll/dt, is presumed to be accurate. At lower redshifts, mergers and the AGN feedback (missing from our model) are expected to be important in determining the SFR. Hence, to normalize our CDFs, we chose the largest redshift at which our sample is reasonably large (see below). It is important to note however that our main results are based on the z > 4 CDFs, which unlike the predictions for the absolute numbers of bursts do not depend on our choice of normalization.

GRB SAMPLE

Our LGRB data are taken from Robertson & Ellis (2012), corresponding to a compilation from the samples presented in Butler et al. (2007), Perley et al. (2009), Butler, Bloom & Poznanski (2010), Sakamoto et al. (2011), Greiner et al. (2011) and Krühler et al. (2011). They include only GRBs present before the end of the second Swift BAT GRB catalogue, and are comprised of 152 LGRBs with redshift measurements. It is important to consider the completeness of the sample. Several efforts have been made to construct a redshift-complete GRB sample (e.g. Greiner et al. 2011; Salvaterra et al. 2012). However, to do so, many GRBs with measured redshifts are excluded. Such requirements are even more severe for high-z bursts, which makes them of little use for our purposes. To explore the dependence of our results on a possible bias in the GRB redshift distribution, we construct two samples: (S1) we use an LF based on the z < 4 subsample (consisting of 136 GRBs) and (S2) we use a subsample with isotropic-equivalent luminosities bright enough to be observable up to high redshifts (comprised of 38 bursts). The two samples are summarized in Table 1. Since there is a degeneracy between a biased SFR–GRB relation and a redshift-dependent LF, we implicitly assume that any unknown bias will be subsumed in the value of the α parameter.

Table 1.

Set of cases considered in our analysis.

Llim(erg s−1)	α = 0	−1 < α < 2	−∞ < α < ∞
Llim ≥ 0	S1C0	S1C1	S1C2
Llim ≥ 1.34 × 1052	S2C0	S2C1	S2C2

Luminosity function sample (S1)

The number of GRBs detectable by any given instrument depends on the specific flux sensitivity threshold and the intrinsic isotropic LF of the GRBs. In Fig. 3, we show the distribution of log luminosities for z < 4 GRBs, which can be well described by a normal distribution

\begin{equation} p(L) = p_*{\rm e}^{-(L-L^*)^2}/2\sigma ^{2}_L. \end{equation}

(18)

The values L = log L_iso/ergs^{− 1}, L* = log 10^51.16, σ_L = 1.06 and p_* = 1.26 are estimated by maximum-likelihood optimization. The luminosity threshold is given by

\begin{equation} L_{\rm lim} = 4\pi \,d_{\rm L}^{2}\, F_{\rm lim}, \end{equation}

(19)

where d_L is the luminosity distance. Consistent with previous works (e.g. Li 2008), we set a bolometric energy flux limit F_lim = 1.15 × 10⁻⁸ erg cm⁻² s⁻¹ for Swift by using the smallest luminosity of the sample. Due to Malmquist bias, our fitted LF is likely biased towards high luminosities. To minimize this problem and increase the sample completeness, we fit our LF using only z ≤ 4 GRBs, which we show in Fig. 3. We reiterate that the Malmquist bias only serves to make our results even more conservative by predicting a flatter redshift distribution of GRBs.

Figure 3.

Frequency (i.e. fraction in bin) of GRB luminosities for the z < 4 subsample used to construct the LF used for S1 (see the text for details). The red dashed line represents the best-fitting LF.

Luminosity-limited sample (S2)

Another approach, less dependent on the LF parametrization and the Malmquist bias, is to construct a luminosity-limited subsample of the observed bursts bright enough to be seen at the highest redshift of interest. Assuming that the LF does not evolve with redshift, this subset would be proportional to the total number of bursts at any given redshift.

In Fig. 4, we show the redshifts and isotropic luminosities of our entire sample. The dot–dashed blue line corresponds to the effective Swift detection threshold. For our luminosity-limited sample, we only use GRBs with isotropic-equivalent luminosities L_iso ≥ 1.34 × 10⁵² erg s^{− 1}, which comprise all GRBs observable up to z ∼ 9.4. Hereafter, all calculations will correspond to either the complete (LF-derived) sample (S1) or the luminosity-limited sample (S2).

Figure 4.

Isotropic luminosities, L_iso, of 152 Swift GRBs as a function of z from the compilation of Robertson & Ellis (2012). The blue dot–dashed line approximates the effective Swift detection threshold (equation 19). The black dashed horizontal line represents the luminosity limit of L_iso > 1.34 × 10⁵² erg s^{− 1}, used to define our S2 subsample.

OBSERVATIONAL CONSTRAINTS

In this section, we test the WDM models by comparing the predicted absolute detection rates of bursts as well as the CDFs with the observed samples. We consider three different ranges of α: (i) a constant SFR–GRB relation, α = 0 (case 0, C0); (ii) −1 < α < 2 (case 1, C1) and (iii) a flat prior over −∞ < α < ∞ (case 2, C2).⁸ All cases are summarized in Table 1.

Absolute detection rate of bursts

In Tables 2 and 3, we present the absolute number of GRBs at high redshifts in CDM and WDM models with particle masses of 0.5–3.5 keV, as well as the actual number in our sample observed with Swift. All models are normalized to yield the observed number of bursts at 3 < z < 4, as described in equation (17) and the associated discussion.

Table 2.

Absolute number of GRBs per redshift interval predicted by each model for the S1C1 sample.

Model	N(4,6)			N(6,8)			N(8,10)
	α = 0	α = 1	α = 2	α = 0	α = 1	α = 2	α = 0	α = 1	α = 2
mx = 0.5 keV	3.18	3.94	4.88	0.01	0.02	0.04	1.0 × 10−5	2.2 × 10−5	4.6 × 10−5
mx = 1.0 keV	9.34	11.82	15.00	0.42	0.71	1.21	0.01	0.02	0.04
mx = 1.5 keV	14.84	19.10	24.67	1.51	2.58	4.42	0.08	0.17	0.36
mx = 2.0 keV	14.31	18.36	23.65	2.31	4.01	6.96	0.20	0.43	0.93
mx = 2.5 keV	14.44	18.54	23.90	2.09	3.65	6.38	0.39	0.83	1.80
mx = 3.0 keV	14.51	18.62	24.01	2.04	3.54	6.16	0.46	1.00	2.20
mx = 3.5 keV	14.56	18.69	24.09	2.04	3.55	6.18	0.44	0.97	2.13
CDM	14.56	18.69	24.10	2.04	3.55	6.18	0.45	1.00	2.19
Swift	11	11	11	3	3	3	2	2	2

Table 3.

Absolute number of GRBs per redshift bin predicted by each model for the S2C1 sample.

Model	N(4,6)			N(6,8)			N(8,10)
	α = 0	α = 1	α = 2	α = 0	α = 1	α = 2	α = 0	α = 1	α = 2
mx = 0.5 keV	1.33	1.65	2.05	0.01	0.02	0.03	1.5 × 10−5	3.1 × 10−5	6.6 × 10−5
mx = 1.0 keV	4.10	5.23	6.70	0.34	0.59	1.00	0.01	0.03	0.06
mx = 1.5 keV	6.75	8.78	11.46	1.27	2.19	3.76	0.12	0.26	0.55
mx = 2.0 keV	6.47	8.40	10.94	2.02	3.52	6.14	0.31	0.67	1.44
mx = 2.5 keV	6.54	8.49	11.07	1.86	3.27	5.77	0.60	1.30	2.81
mx = 3.0 keV	6.57	8.53	11.12	1.79	3.13	5.49	0.74	1.63	3.58
mx = 3.5 keV	6.59	8.56	11.16	1.80	3.14	5.51	0.72	1.58	3.48
CDM	6.60	8.56	11.16	1.80	3.15	5.51	0.74	1.63	3.58
Swift	6	6	6	3	3	3	2	2	2

As expected, models with small WDM particle masses predict a rapidly decreasing GRB rate towards high redshifts. This exponential suppression can in some cases be partially compensated by an increasing GRB-to-SFR rate (i.e. α > 0). For the S1C1 case, models with m_x ≥ 2.5 keV show good agreement with Swift observations for 0 < α < 1, though values of α ∼ 2 are a better fit to the observations at z > 8. For the case S2C1, all models with m_x ≥ 2.5 keV seem to be consistent with data for α ∼ 1-2. In both cases, the two observed bursts in the interval 8 < z < 10 are already at odds with 1.5 <m_x < 2.5 keV models. Finally, we see that models with m_x ≤ 1 keV predict a dearth of GRBs at z > 6, which is inconsistent with current observations. Extreme models with m_x ∼ 0.5 keV already fail at intermediate redshifts (4 < z < 6), even for values of α as high as 2.

The redshift distribution of z > 4 bursts

Although the absolute rate of bursts is the simplest prediction, it is dependent on the normalization factor between the SFR and GRB rate at 3 < z < 4. Hence, for the remainder of the paper, we focus on comparing the theoretical and observed z > 4 CDFs. The CDFs are not dependent on normalization factors and are therefore more conservative and robust predictions.

In Fig. 5, we plot the CDFs for CDM and WDM (under the assumption of α = 0), as well as the observed Swift distribution. The lighter the WDM particle, the sharper the CDF rise at low z. There is a clear separation between CDM and WDM models with m_x ≲ 1.5 keV. Both the S1 and S2 samples (top and bottom panels, respectively) show the same qualitative trends.

Figure 5.

Cumulative number of GRBs for different values of m_x compared with CDM predictions and Swift observations. The blue dotted line corresponds to m_x = 0.5 keV, green dotted line to m_x = 1.0 keV, red dotted line to m_x = 1.5 keV, purple dotted line to m_x = 2.3 keV, brown dotted line to m_x = 2.5 keV, orange dotted line to m_x = 3.0 keV, cyan dotted line to m_x = 3.5 keV, dark-green two-dashed line to CDM and black to the Swift observations. Top panel: sample S1C0; Bottom panel: sample S2C0.

As we saw above, the high-z suppression of structures in WDM models can be compensated for by allowing the GRB rate/SFR to increase towards higher redshifts. How degenerate are these cosmological versus astrophysical effects? In Fig. 6, we show the CDF for m_x = 0.5 keV for several values of α for the S2 sample. The exponential suppression of DM halo abundances in this model is so strong that an unrealistically high value of α ∼ 15 is required to be roughly consistent with observations. Such a high value is ruled out by low-redshift observations, which imply α ≲ 1 (e.g. Kistler et al. 2009; Robertson & Ellis 2012; Trenti et al. 2012).

Figure 6.

Cumulative number of GRBs for m_x = 0.5 keV as a function of the α parameter. The blue dotted line represents α = 0, green dotted line α = 3, red dotted line α = 6, purple dotted line α = 9, brown dotted line α = 12, orange dotted line α = 15, cyan dotted line α = 18, dark-green two-dashed line CDM and black line the Swift observations.

Constraints from the redshift distribution of z > 4 bursts

To quantify how consistent are these CDFs with the observed distribution from Swift, we make use of two statistics: (i) the one-sample Kolmogorov–Smirnov (K–S) test and (ii) a maximum-likelihood estimation (MLE). Both tests are described in detail in Appendix A.

The K–S test provides a simple estimate of the probability that the observed distribution was drawn from the underlying theoretical one. We compute this probability, for fixed α first, for our models S1C0, S1C1, S2C0 and S2C1. Consistent with the more qualitative analysis from the previous section, models with m_x ≲ 1.0 keV are ruled out at 90 per cent CL assuming −1 ≤ α ≤ 1. For α = 0 (S1C0 and S2C0), the limits are even more restrictive and models with m_x ≲ 1.5 keV are ruled out at 90 per cent CL for both samples.

So far, we have analysed each model individually in order to quantify a lower limit on m_x, given a single value of α. Using a χ² MLE (see Appendix A) allows us to compute posterior probabilities given conservative priors on α. Thus, we are able to construct confidence limits in the two-dimensional (m_x, α) parameter space. The results for cases S1C2 and S2C2 are shown in Fig. 7 at 68, 95 and 99 per cent CL. Both samples show the same qualitative trends, with the data preferring higher values of m_x and CDM. Marginalizing the likelihood over −3 ≤ α ≤ 12, with a flat prior, shows that models with m_x ≤ 1.6-1.8 keV are ruled out at 95 per cent CL for S1C2 and S2C2, respectively.

Figure 7.

Contours over α and m_x, enclosing 68 (green), 95 (orange) and 99 per cent (grey) probability. The asterisks correspond to the best-fitting parameter combinations. Top panel: sample S1C2; bottom panel: sample S2C2. The horizontal dotted lines represent the values α = 0 (Ishida, de Souza & Ferrara 2011; Elliott et al. 2012), α = 0.5 (Robertson & Ellis 2012), α = 1.2 (Kistler et al. 2009) and α = 2 for comparison.

FUTURE CONSTRAINTS

In the previous section, we have quantified the constraints on WDM particle masses using current Swift GRB observations. We obtain constraints of m_x ≳ 1.6–1.8 keV. We now ask how much could these constraints improve with a larger GRB sample, available from future missions? As a reference, we use the Sino-French space-based multiband astronomical variable objects monitor (SVOM)⁹ mission. The SVOM has been designed to optimize the synergy between space and ground instruments. It is forecast to observe ∼ 70-90 GRBs yr^{− 1} and ∼ 2-6 GRB yr^{− 1} at z ≥ 6 (see e.g. Salvaterra et al. 2008).

We first construct a mock GRB data set of 450 bursts with redshifts obtained by sampling the CDM, α = 0 probability density function (PDF) given by equation (13). This sample size represents an optimistic prediction for 5 yr of SVOM observations¹⁰ (see e.g. Salvaterra et al. 2008). We then perform the MLE analysis detailed above on this mock data set at z > 4. The resulting confidence limits are presented in Fig. 8.

Figure 8.

Same as Fig. 7, but assuming a 450-burst mock sample, drawn from the CDM, α = 0 PDF.

This figure shows that ∼5 yr of SVOM observations would be sufficient to rule out m_x ≤ 2.3 keV models (from our fiducial CDM, α = 0 model) at 95 per cent CL, when marginalized over α. This is a modest improvement over our current constraints using Swift observations. As already foreshadowed by Figs 1 and 2, as well as the associated discussion, it is increasingly difficult to push constraints beyond m_x > 2 keV. On the other hand, the α constraint improves dramatically due to having enough high-z bursts to beat the Poisson errors. We caution that the relative narrowness around α = 0 of the contours in Fig. 8 is also partially due to our choice of (CDM, α = 0) as the template for the mock observation.

CONCLUSION

Small-scale structures are strongly suppressed in WDM cosmologies. WDM particle masses of m_x ∼ keV have been invoked in order to interpret observations of local dwarf galaxies and galactic cores. The high-redshift Universe is a powerful testbed for these cosmologies, since the mere presence of collapsed structures can set strong lower limits on m_x. GRBs, being extremely bright and observable to well within the first billion years, are a promising tool for such studies.

Here, we model the collapsed fraction and cosmic SFR in CDM and WDM cosmologies, taking into account the effects of both free-streaming and effective pressure due to the residual velocity dispersion of WDM particles. Assuming that the GRB rate is proportional to the SFR, we interpret 5 yr of Swift observations in order to place constraints on m_x. We conservatively account for astrophysical uncertainty by allowing the GRB rate/SFR to evolve with redshift as ∝(1 + z)^α. In order to fold completeness limits into our analysis, we used a low-z sample to estimate the intrinsic LF, or else restricted our analysis to a luminosity-limited subsample detectable at all redshifts.

For each model (m_x, α), we compute both the absolute detection rates and CDFs, at z > 4. A K–S test between the model and observed CDFs rules out m_x < 1.5 (1.0) keV, assuming α = 0 (<2), at 90 per cent CL. Using a maximum-likelihood estimator, we are able to marginalize over α. Assuming a flat prior in α, we constrain m_x > 1.6–1.8 keV at 95 per cent CL. A future SVOM-like mission would tighten these constraints to m_x ≳ 2.3 keV.

The strong and robust constraints that we derive show that GRBs are a powerful probe of the early Universe. Their utility would be further enhanced with insights into their formation environments and their relation to the cosmic SFR.

We thank Emille Ishida for the careful and fruitful revision of the draft of this work and Andressa Jendreieck for useful comments. RSS thanks the Max Planck Institute for Astrophysics (Garching, Germany) for its hospitality during his visit.

1

Throughout the text, we use ‘conservative’ to imply biasing the GRB distribution towards lower redshifts. This is conservative since it mimics the effects of WDM, thereby resulting in weaker constraints on the particle mass.

2

Throughout the paper, we adopt the cosmological parameters from Wilkinson Microwave Anisotropy Probe 9 (Hinshaw et al. 2012); Ω_m = 0.264, Ω_Λ = 0.736, n_s = 0.97, σ₈ = 0.8 and H₀ = 71 km s⁻¹ Mpc ⁻¹.

3

http://camb.info/

4

Using M_gal(z) from Sobacchi & Mesinger (2013a) for WDM models is not entirely self-consistent, since UVB feedback effects should be smaller in WDM models. However, this effect is smaller than other astrophysical uncertainties. Most importantly, our conclusions are driven by models which at high redshift have M_WDM > M_gal, and are therefore insensitive to the actual choice of M_gal (see Fig. 2).

5

Since host galaxies of long-duration GRBs are often observed to be metal poor, several studies have tried to connect the origin of LGRBs with the metallicity of their progenitors (e.g. Mészáros 2006; Woosley & Bloom 2006; Salvaterra & Chincarini 2007; Salvaterra et al. 2009; Campisi et al. 2011). Such a connection is physically motivated since core-collapse models could not generate an LGRB without the progenitor system having low metallicity (e.g. Hirschi, Meynet & Maeder 2005; Yoon & Langer 2005; Woosley & Bloom 2006). On the other hand, several authors report observations of GRBs in high-metallicity environments (e.g. Levesque et al. 2010; Krühler et al. 2012), suggesting that GRB hosts are not necessarily metal poor. Despite the apparent preference of GRBs towards metal-poor hosts, there is no clear cutoff in metallicity, above which GRB formation should be suppressed.

6

There are several selection effects known to mask the true GRB redshift distribution, e.g. (i) the host galaxy dust extinction; (ii) the redshift desert (a redshift span, 1.4 < z < 2.5, in which it is difficult to measure absorption and emission spectra); (iii) Malmquist bias and (iv) the difference between redshift measurements techniques (e.g. Coward et al. 2012). We therefore expect K to be redshift dependent, and this evolution is subsumed in our parameter α above. Most of the above effects (e.g. obtaining GRB redshifts, dust extinction and Malmquist bias) result in biasing the observed sample towards low redshifts. Since we calibrate the proportionality between the GRB rates and SFRs at z ∼ 3–4, we likely overestimate the efficiency in the redshift determination of z > 4 GRBs. Therefore, we expect that even a non-evolving K (i.e. our results for α = 0) would be a conservative assumption.

7

Which in our case corresponds to 5 yr of observation time by Swift. Again, the value is not relevant for our purposes, since it is subsumed in the normalization.

8

More precisely, C2 was run over the interval −3 < α < 12, which was more than sufficient to capture the likelihood decreasing to → 0 in the tails of the distribution (see Fig. 7).

9

http://www.svom.fr

10

The highest redshift in our mock sample is z_max = 11.6, being the only event at z > 10.

11

The value Δz = 1.5 is chosen to ensure at least one burst per bin.

APPENDIX A: STATISTICS

Kolmogorov–Smirnov one-sample test

A straightforward way to compare the data and WDM models is to perform a one-sample K–S test. The null hypothesis that the observed GRB redshifts are consistent with a model distribution can be evaluated by estimating a p-value, which corresponds to one minus the probability that the null hypothesis can be rejected. The K–S test consists in comparing the statistical parameter

\begin{equation} D = \mathrm{sup}|F(z)-G(z)|, \end{equation}

(A1)

where F(z) and G(z) are the CDF for the theoretical and observed sample and sup is the supremum of a totally or partially ordered set. We estimate the p-value via non-parametric bootstrap, which consists of running Monte Carlo realizations of the observed CDF using a random-selection-with-replacement procedure estimated from the data. This provides a histogram of the statistic D, from which a valid goodness-of-fit probability can be evaluated. The probability distribution function for each model is determined by equation (16).

Maximum-likelihood method

More formally, we can estimate the probability of parameters {α, m_x} given the observed data using a Bayesian technique. Considering that the data are described by the probability density function |$f(x;\pmb {\theta })$|⁠, where x is a variable and |$\pmb {\theta } \equiv \lbrace \alpha ,m_{\rm x}\rbrace$|⁠. We want to estimate |$\pmb {\theta }$|⁠, assuming that the data are independent. So the likelihood will be given by

\begin{equation} {\cal L}(\pmb {\theta }|x_i)\propto \prod _{i=1}^{N}f(x_i|\pmb {\theta }). \end{equation}

(A2)

Given the small number of observed bursts per redshift bin, Δz = 1.5, we use a Poisson error statistics¹¹ (see e.g. Campanelli et al. 2012 for a similar procedure applied to the galaxy cluster number count). Therefore, the likelihood function can be computed as

\begin{equation} {\cal L}(\pmb {\theta }|\kappa _i)\propto \prod _{i=1}^{5}\frac{\Upsilon _i^{\kappa _i}{\rm e}^{-\Upsilon _i}}{\kappa _i!}, \end{equation}

(A3)

where |$\Upsilon _i \equiv N(z_i,z_{i+1};\pmb {\theta })$| and κ_i ≡ N_obs(z_i, z_i+1). Thus, the χ² statistics can be written as

\begin{eqnarray} \chi ^2(\pmb {\theta }) &=& -2\ln {\cal L},\nonumber \\ &=& 2\sum _{i=1}^{5}\Upsilon _i-\kappa _i(1+\ln {\Upsilon _i}-\ln {\kappa _i)}. \end{eqnarray}

(A4)