Magnetar Central Engine and Possible Gravitational Wave Emission of Nearby Short GRB 160821B

GRB 160821B triggered the Burst Alert Telescope (BAT) at 22:29:13 UT on 2016 August 21 (Siegel et al. 2016). We developed an IDL script to automatically download the Swift BAT data. We use the standard HEASOFT tools (version 6.12) to process the data. We run bateconvert from the HEASOFT software release to obtain the energy scale for the BAT events. The light curves and spectra are extracted by running batbinevt (Sakamoto et al. 2007). Then, we calculate the cumulative distribution of the source counts using the arrival time of a fraction between 5% and 95% of the total counts to define T₉₀. The time bin size is fixed to 64 ms in this case, due to the short duration. The background is extracted using two time intervals, one before and one after the burst. We model the background as Poisson noise, which is the standard background model for prompt emission. We invoked Xspec to fit the spectra. For technical details, please refer to Sakamoto et al. (2007). XRT began observing the field 57 s after the BAT trigger (Siegel et al. 2016). We made use of public data from the Swift archive.⁷ The Ultra-Violet Optical Telescope (UVOT; Roming et al. 2005) observed the field at T₀ + 76 s, but no optical afterglow was consistent with the XRT position (Evans et al. 2016). There was also no detection in the initial UVOT exposures (Xu et al. 2016). Upper limits of 3σ on the data are preliminarily obtained by using the UVOT photometric system for the first finding chart (FC) exposure (Breeveld & Siegel 2016). On the other hand, r- and z-band afterglow images were obtained using the William Herschel Telescope on La Palma (Levan et al. 2016). In the spectrum of the candidate host galaxy, several prominent emission lines were found (Hβ, [O iii], and Hα), at a redshift of z = 0.16. The physical offset of the afterglow from the candidate host galaxy is approximately 15 kpc (Levan et al. 2016).

The Fermi Gamma-ray Burst Monitor (GBM) triggered and located GRB 160821B at 22:29:13.33 UT on 21 August 2016 (Stanbro & Meegan 2016). GBM has 12 sodium iodide (NaI) and two bismuth germanate (BGO) scintillation detectors, covering the energy range 8 keV–40 MeV (Meegan et al. 2009). We downloaded GBM data of this GRB from the public science support center at the official Fermi Web site.⁸ Each of the GBM detectors collected data of three different types: CTIME, CSPEC, and TTE. Then, we extracted the light curves and performed spectral analysis based on the package gtBurst. By invoking the heasoft command fselect and the ScienceTools command gtbin, we extracted light curves with a time bin of 64 ms in a user-specified energy band from the GBM. We clicked "Tasks $\to$ Make spectra for XSPEC" in gtBurst to extract the source spectrum of the GBM data. The background spectra are extracted from the time intervals before and after the prompt emission phase and modeled with a polynomial function, and the source spectrum is extracted by applying the background model to the prompt emission phase. This GRB occurred right at the edge of the Large Area Telescope field of view (Atwood et al. 2009), about 61° from boresight, so we do not expect to detect the high-energy signal if it exists.

As shown in Figure 1, the BAT light curve shows a single short peak with duration T₉₀ = 0.48 ± 0.07 s, and there is no evidence of extended emission detected in the BAT energy range up to 100 s after the BAT trigger (T₀). The time-integrated BAT spectrum can be fit by a single power law with photon index Γ_γ = 1.98 ± 0.11. The BAT band (15−150 keV) peak flux is (1.7 ± 0.2) photons cm⁻² s⁻¹, and the total fluence is (1.1 ± 0.1) × 10⁻⁷ erg cm⁻² (Palmer et al. 2016). The GBM light curve of GRB 160821B is also shown in Figure 1 with a 64 ms time bin. The profile of the light curve is similar to that of BAT data, having a bright peak with duration T₉₀ ∼ 1.2 s in the 8−1000 keV range. The GBM spectra can be fit by a power-law function due to lack high-energy photons.⁹ Here, we jointly fit the time-averaged spectra of Fermi/GBM+Swift/BAT with a power-law model, and obtained ${{\rm{\Gamma }}}_{\mathrm{BAT}+\mathrm{GBM}}\,=\,1.88\pm 0.12$ with χ² = 0.84 (Figure 2). The total fluence in the 8–1000 keV range is (2.52 ± 0.19) × 10⁻⁶ erg cm⁻². The isotropic energy ${E}_{\gamma ,\mathrm{iso}}\,=\,(2.1\pm 0.2)\times {10}^{50}\,\mathrm{erg}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Norris et al. (2000) discovered an anti-correlation between the GRB peak luminosity and the delay time (t_lag) in different energy bands, meaning softer photons usually arrive later than harder photons. This spectral lag is always significant in long-duration GRBs (Norris et al. 2000; Gehrels et al. 2006; Liang et al. 2006), but is typically negligible in short-duration GRBs (Norris & Bonnell 2006; Zhang et al. 2009). We extracted 4 ms binned light curves in the following three BAT energy bands: 15–25 keV, 25–50 keV, and 50–100 keV.¹⁰ Then, we used the cross-correlation function method (Norris et al. 2000; Ukwatta et al. 2010) to calculate the lags between the 25–50 keV and 50–100 keV light curves in order to address the questions of whether the short GRB 160821B is consistent with typical type I and other short-hard GRBs within the spectral lag distribution. Figure 4 shows the peak luminosity as a function of spectral lag for typical type II, type I, and other short-hard GRBs and GRB 160821B. Here, type II GRBs correspond to a confirmed supernova (SN) association, or have a high specific star formation rate (SSFR) and do not have a large offset from the galaxy center. On the contrary, type I GRBs occur in elliptical or early-type host galaxy without an SN signature or that has a relatively low local SSFR and large offset from the host galaxy center. The other short-hard GRBs do not satisfy either of the two criteria of the type I sample and do not have their host galaxy identified, but these GRBs have a short duration and hard spectra (Zhang et al. 2009). The type II GRB samples are from Norris & Bonnell (2006), and give the best power-law model fit, with a 2σ region of the fitting. Type I and other short-hard GRBs are collected from Zhang et al. (2009). For GRB 160821B, we found t_lag = (10 ± 6) ms, which is consistent with the result of Palmer et al. (2016). It deviates from the 2σ region of the type II GRB fitting, but both the peak luminosity and lag value of this case are comparable with those of type I GRBs.

The initial X-ray light curve is best fit by a broken power law, which reads

$$\begin{eqnarray}&&F\,=\,{F}_{0}{\left[{\left(\displaystyle \frac{t}{{t}_{{\rm{b}}}}\right)}^{\omega {\alpha }_{1}}+{\left(\displaystyle \frac{t}{{t}_{{\rm{b}}}}\right)}^{\omega {\alpha }_{2}}\right]}^{-1/\omega },\end{eqnarray} \tag{ 1 }$$

where ω describes the sharpness of the break and is taken to be 3 in this analysis (Liang et al. 2007). There is an initial decay slope α₁ = 0.21 ± 0.14, followed by a steeper decay of α₂ = 4.52 ± 0.45 with a break time t_b = 180 ± 46 s after the BAT trigger. No significant X-ray flare was detected during the observation time. The X-ray spectrum in the 0.3–10 keV energy band is best fit by an absorbed power law with ${{\rm{\Gamma }}}_{X}\,=\,{1.95}_{-0.08}^{+0.21}$ and column density N_H = (7.5 ± 2.1) ×10²⁰ cm⁻². The X-ray light curve, along with the Γ_X evolution, is shown in Figure 3. About 1000 s after the BAT trigger, another component emerged, which is likely normal decay and a post-jet break. We used a broken power law to fit this component and found α₃ ∼ 0.45, α₄ ∼ 3.5 with a break time around 35,000 s. We follow the method discussed in Zhang et al. (2007) to calculate ${E}_{{\rm{K}},\mathrm{iso}}$, which is almost constant during the normal decay phase in the X-ray afterglow. We assume that it is in the ν > max(ν_m, ν_c) region, where the afterglow flux expression does not depend on the medium density. In our calculations, the microphysics parameters of the shock are assigned standard values derived from observations, i.e., ε_e = 0.01 and ε_B = 0.001, and thus ${E}_{{\rm{K}},\mathrm{iso}}\sim 8\times {10}^{52}\,\mathrm{erg}$. If the later break is assumed to be a post-jet break,¹¹ one can estimate the jet opening angle θ_j ∼ 0.063 rad ∼ 3.6 with medium density n = 0.1 cm⁻³.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The XRT light curve of the short GRB 160821B, which has a short plateau emission following an abrupt decay, is unusual, but not odd. This temporal behavior is similar to that of the short GRB 090515 (Rowlinson et al. 2010). In Figure 5(a), we collect all of the short GRB light curves without extended emission to compare with the X-ray emission of GRB 160821B. We found that most short GRBs appear to have a steeper decay around several hundred seconds. In particular, the X-ray emission behavior of GRB 160821B is similar to that of GRB 090515, which is the first short GRB claimed to have a magnetar central engine origin (Rowlinson et al. 2010). Also, the plateau flux of GRB 160821B is the highest among other short GRBs. In Figure 5(b), we show the fluence in the BAT band (15–150 keV) and the flux in the XRT band (0.3–10 keV) at T₀ + 100 s for all of the short GRBs in the Swift sample. GRB 160821B is shown as a filled circle. As expected, the higher fluence GRBs tend to have higher flux X-ray afterglows.

Zoom In Zoom Out Reset image size

According to Zhang & Mészáros (2001), the characteristic spin-down luminosity L₀ and timescale τ are written as

$$\begin{eqnarray}&&{L}_{0}\,=\,1.0\times {10}^{49}\,\mathrm{erg}\,{{\rm{s}}}^{-1}({B}_{{\rm{p}},15}^{2}{P}_{0,-3}^{-4}{R}_{6}^{6}),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\tau \,=\,2.05\times {10}^{3}\,{\rm{s}}\,({I}_{45}{B}_{{\rm{p}},15}^{-2}{P}_{0,-3}^{2}{R}_{6}^{-6}),\end{eqnarray} \tag{ 3 }$$

where I is the moment of inertia of a typical NS with mass M_NS = 1.4 M_⊙, P₀ is the initial spin period, B_p is the magnetic field strength, R is the radius of the NS, and the convention Q = 10^x Q_x is adopted in cgs units for all other parameters throughout the paper. The spin-down timescale can be generally identified as the lower limit of the observed break time, i.e.,

$$\begin{eqnarray}&&\tau \gt {t}_{b}/(1+z),\end{eqnarray} \tag{ 4 }$$

where t_b is the break time after the internal plateau, obtained using a broken power-law function fitting. The redshift z = 0.16 is adopted. The bolometric luminosity at the break time t_b is

$$\begin{eqnarray}&&{L}_{b}\,=\,4\pi {D}_{L}^{2}{F}_{b}\cdot k,\end{eqnarray} \tag{ 5 }$$

where F_b = (1.6 ± 0.82) × 10⁻⁹ erg cm⁻² s⁻¹ is the X-ray flux at t_b, D_L² is luminosity distance, and k is the k-correction factor. The characteristic spin-down luminosity is essentially the plateau luminosity, which may be estimated as

$$\begin{eqnarray}&&{L}_{0}\simeq {L}_{{\rm{b}}}\,=\,(1.8\pm 0.6)\times {10}^{48}\,\mathrm{erg}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 6 }$$

Based on Equations (2) and (3), one can derive the magnetar parameters ${B}_{{\rm{p}}}$ and P₀:

$$\begin{eqnarray}&&{B}_{{\rm{p}},15}\,=\,2.05({I}_{45}{R}_{6}^{-3}{L}_{0,49}^{-1/2}{\tau }_{3}^{-1})\,{\rm{G}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{P}_{0,-3}\,=\,1.42({I}_{45}^{1/2}{L}_{0,49}^{-1/2}{\tau }_{3}^{-1/2})\,{\rm{s}}.\end{eqnarray} \tag{ 8 }$$

Using the lower limit of τ we derive upper limits for P₀ and B_p, which are, respectively, P₀ < 8.5 × 10⁻³ s and B_p < 3.12 ×10¹⁶ G. Figure 6(a) shows the B_p−P₀ diagram for GRB 160821B, and compares it other short GRBs.

Zoom In Zoom Out Reset image size

Another relevant timescale is the collapse time of a supra-massive magnetar, t_col. The post-internal plateau decay slope α₂ is steeper than 2, which is the standard spin-down luminosity evolution with time (Zhang & Mészáros 2001). The break time is therefore defined by the collapse time t_col, and one can write

$$\begin{eqnarray}&&\tau \simeq {t}_{{\rm{b}}}.\end{eqnarray} \tag{ 9 }$$

The maximum gravitational mass (M_max) depends on the spin period, which increases with time. Using the same method described in Lasky et al. (2014) and Lü et al. (2015), we can write down t_col as a function of M_p^:

$$\begin{eqnarray}{t}_{\mathrm{col}} & = & \displaystyle \frac{3{c}^{3}I}{4{\pi }^{2}{B}_{{\rm{p}}}^{2}{R}^{6}}\left[{\left(\displaystyle \frac{{M}_{{\rm{p}}}-{M}_{\mathrm{TOV}}}{\hat{\alpha }{M}_{\mathrm{TOV}}}\right)}^{2/\hat{\beta }}-{P}_{0}^{2}\right]\\ & = & \displaystyle \frac{\tau }{{P}_{0}^{2}}\left[{\left(\displaystyle \frac{{M}_{{\rm{p}}}-{M}_{\mathrm{TOV}}}{\hat{\alpha }{M}_{\mathrm{TOV}}}\right)}^{2/\hat{\beta }}-{P}_{0}^{2}\right].\end{eqnarray} \tag{ 10 }$$

where $\hat{\alpha }$, $\hat{\beta }$, and M_TOV are dependent on the EOS. Therefore, we can use t_col to constrain the NS EOS. Here, we only consider five EOSs (SLy, APR, GM1, AB-N, and AB-L) for the given protomagnetar mass distribution derived from the total mass distribution of Galactic NS−NS binary systems (Figure 6(b)). (1) SLy: effective nuclear interaction of neutron-rich matter with M_TOV = 2.05 M_⊙ and R = 9.97 km. (2) APR: assume that the inner material includes both dense nucleon admixture of quark matter with M_TOV = 2.20 M_⊙ and R = 10.00 km. (3) GM1: relates scalar and vector couplings of the hyperons for saturated nuclear matter with M_TOV = 2.37 M_⊙ and R = 12.05 km. (4) AB-N: neutrons' nuclear attraction due to pion exchange tensor with M_TOV = 2.67 M_⊙ and R = 12.90 km. (5) AB-L: neutrons' nuclear attraction due to scalar exchange with M_TOV = 2.71 M_⊙ and R = 13.70 km (Lasky et al. 2014). Our results show that the GM1 model gives an M_p band falling within the 2σ region of the protomagnetar mass distribution, such that the GM1 EOS is the best candidate for a non-rotating NS with maximum mass M_TOV = 2.37 M_⊙.

One of the most important necessary conditions of a magnetar central engine candidate for GRBs is that the sum of the prompt emission energy (${E}_{\gamma ,\mathrm{iso}}$), internal plateau energy (${E}_{{\rm{X}},\mathrm{iso}}$), and kinetic energy (${E}_{{\rm{K}},\mathrm{iso}}$) after jet correction should be less than the total rotation energy (energy budget of magnetar) if we assume the magnetar wind is isotropic. The total rotation energy of the millisecond magnetar is

$$\begin{eqnarray}&&{E}_{\mathrm{rot}}\,=\,\displaystyle \frac{1}{2}I{{\rm{\Omega }}}_{0}^{2}\simeq 2\times {10}^{52}\,\mathrm{erg}\,{M}_{1.4}{R}_{6}^{2}{P}_{0,-3}^{-2},\end{eqnarray} \tag{ 11 }$$

where Ω₀ = 2π/P₀ is the initial angular frequency of the neutron star. The total energy of the prompt emission is ${E}_{\gamma ,\mathrm{iso}}\,=\,(2.1\pm 0.2)\times {10}^{50}\,\mathrm{erg}$ within the energy range 8–1000 keV. The X-ray internal plateau energy can be roughly estimated using the break time and break luminosity (Lü & Zhang 2014), i.e.,

$$\begin{eqnarray}{E}_{{\rm{X}},\mathrm{iso}} & \simeq & {L}_{b}\cdot \displaystyle \frac{{t}_{b}}{1+z}\\ & \simeq & (2.79\pm 0.9)\times {10}^{49}\,\mathrm{erg}.\end{eqnarray} \tag{ 12 }$$

To estimate the kinetic energy ${E}_{{\rm{K}},\mathrm{iso}}$, which is used in the standard forward afterglow model, one has ${E}_{{\rm{K}},\mathrm{iso}}\sim 8\,\times {10}^{52}\,\mathrm{erg}$ (see Section 2.3). Therefore ${E}_{\mathrm{rot}}\gg \tfrac{1}{2}{\theta }_{{\rm{j}}}^{2}({E}_{\gamma ,\mathrm{iso}}\,+{E}_{{\rm{X}},\mathrm{iso}}+{E}_{{\rm{X}},\mathrm{iso}})$, which satisfies the magnetar central engine energy budget requirement.

The coalescence of double neutron stars is believed to be one of the most likely sources for powering GW radiation with associated EM signals. These events have promising detectability prospects with current and future GW detectors like Advanced LIGO/Virgo (Fan et al. 2013; Gao et al. 2013; Yu et al. 2013; Zhang 2013). If indeed a magnetar drove GRB 160821B, why was the total rotation energy of the magnetar much larger than the sum of the prompt emission energy, internal energy, and kinetic energy? Several possible reasons may be used to interpret the gap in the energy compared to the energy budget. One is that the efficiency is as low as ∼0.01, and such low efficiency may disfavor the magnetic energy dissipation process (Fan et al. 2013). Another possibility is the missing energy must have been carried away by non-electromagnetic GW radiation (Fan et al. 2013; Ho 2016; Lasky & Glampedakis 2016) or carried to the black hole before spin down.

Following Fan et al. (2013) and Lasky & Glampedakis (2016), a magnetar loses rotational energy through two channels: magnetic dipole torques (L_m) GW radiation (L_w):

$$\begin{eqnarray}-\,{{dE}}_{\mathrm{rot}}/{dt} & = & {L}_{{\rm{m}}}+{L}_{{\rm{w}}}\\ & = & \displaystyle \frac{{B}_{{\rm{p}}}^{2}{R}^{6}{{\rm{\Omega }}}^{4}}{6{c}^{3}}+\displaystyle \frac{32{{GI}}^{2}{\epsilon }^{2}{{\rm{\Omega }}}^{6}}{5{c}^{5}},\end{eqnarray} \tag{ 13 }$$

where = 2(I_xx − I_yy)/(I_xx + I_yy) is the ellipticity in terms of the principal moments of inertia, assuming the magnetar has a triaxial shape. Following the method of Lasky & Glampedakis (2016), GW radiation can be more efficient than magnetic dipole radiation because of its stronger dependence on the neutron star spin rate Ω, i.e., Ω⁶ and Ω⁴ respectively. The upper limit on the ellipticity () can be expressed simply as a dependence on the observed plateau luminosity and break time (Lasky & Glampedakis 2016),

$$\begin{eqnarray}{\epsilon }_{\mathrm{obs}} & \leqslant & {\left(\displaystyle \frac{15{c}^{5}\eta I}{512{{GL}}_{0}^{2}{t}_{{\rm{b}}}^{3}}\right)}^{1/2}\\ & = & 0.33\eta {\left(\displaystyle \frac{I}{{10}^{45}{\rm{g}}{\mathrm{cm}}^{2}}\right)}^{1/2}{\left(\displaystyle \frac{{L}_{0}}{{10}^{49}\mathrm{erg}{{\rm{s}}}^{-1}}\right)}^{-1}{\left(\displaystyle \frac{{t}_{{\rm{b}}}}{100{\rm{s}}}\right)}^{-3/2}.\end{eqnarray} \tag{ 14 }$$

Using the typical NS mass and radius, η = 0.1, L₀ ∼ 5.38 ×10⁴⁷ erg s⁻¹, and t_b ∼ 180 s, one has _obs < 0.07.

If most of the rotation energy is released via GW radiation with a frequency f, the GW strain for a rotating neutron star at distance D_L can be expressed as

$$\begin{eqnarray}&&h{(t)=\displaystyle \frac{4{GI}\epsilon }{{D}_{{\rm{L}}}{c}^{4}}{\rm{\Omega }}(t)}^{2}.\end{eqnarray} \tag{ 15 }$$

The noise power spectral density of the detector, S_h(f), and the stationary phase approximation implies $\tilde{h}{(f)}^{2}\,=\,h{(t)}^{2}| {dt}/{df}|$, where $\tilde{h}(f)$ is the Fourier transform of h(t). Following the method of Lasky & Glampedakis (2016), $\tilde{h}(f)$ can be expressed as

$$\begin{eqnarray}\tilde{h}(f) & = & \displaystyle \frac{1}{{D}_{{\rm{L}}}}\sqrt{\displaystyle \frac{5{GI}}{2{c}^{3}f}}\\ & \approx & 2.6\times {10}^{-25}{\left(\displaystyle \frac{I}{{10}^{45}{\rm{g}}{\mathrm{cm}}^{2}}\displaystyle \frac{1\mathrm{kHz}}{f}\right)}^{1/2}{\left(\displaystyle \frac{{D}_{{\rm{L}}}}{100\mathrm{Mpc}}\right)}^{-1}.\end{eqnarray} \tag{ 16 }$$

So $\tilde{h}(f)$ is independent of the neutron star ellipticity, but depends on the angular frequency evolution with time. The characteristic amplitude ${h}_{{\rm{c}}}\,=\,{fh}(t)\sqrt{{dt}/{df}}\,=\,f\tilde{h}(f)$ (Corsi & Mészáros 2009; Hild et al. 2011) is

$$\begin{eqnarray}{h}_{{\rm{c}}} & = & \displaystyle \frac{f}{{D}_{{\rm{L}}}}\sqrt{\displaystyle \frac{5{GI}}{2{c}^{3}f}}\\ & \approx & 8.22\times {10}^{-24}{\left(\displaystyle \frac{I}{{10}^{45}{\rm{g}}{\mathrm{cm}}^{2}}\displaystyle \frac{f}{1\mathrm{kHz}}\right)}^{1/2}{\left(\displaystyle \frac{{D}_{{\rm{L}}}}{100\mathrm{Mpc}}\right)}^{\mathrm{-}1}.\end{eqnarray} \tag{ 17 }$$

For GRB 160821B, its redshift z = 0.16 corresponds to D_L = 765 Mpc. Using this and f = 1000 Hz, one can estimate the maximum value of the strain h_c, which is less than 1.1 × 10⁻²⁴. In Figure 7, we plot the GW strain sensitivity for Advanced LIGO and the Einstein Telescope (ET), from Figure 3 of Lasky & Glampedakis (2016). It is clear that the strain of GRB 160821B is below the initial LIGO or Advanced-LIGO noise curve. However, it is comparable to the proposed detectability limit of ET and such a signal may be detected by ET in the future.

Zoom In Zoom Out Reset image size

On the other hand, keeping the total energy constant and moving the event to a lower redshift allows one to estimate the minimum detectability distance of such an event. The GW strain amplitude will be stronger if the event occurs at a lower redshift. One can estimate the cosmological distances where the GW signal can be detected by the current Advanced LIGO. We simulate this source at different distances and calculate the GW strain amplitude. Then we compare that value with the current sensitivity of Advanced LIGO. We find that this GW signal could be detected if shifted to about 100 Mpc, which corresponds to redshift z ∼ 0.023 (Figure 7).