Constraining the Long-lived Magnetar Remnants in Short Gamma-Ray Bursts from Late-time Radio Observations

The most promising models for short gamma-ray bursts (SGRBs) are mergers of two compact objects, such as double neutron stars (NS–NS) or a neutron star–black hole systems. The detection of the gravitational-wave event GW170817 from an NS–NS merger (Abbott et al. 2017a), and its associated SGRB 170817A (Abbott et al. 2017b; Goldstein et al. 2017; Zhang et al. 2018) unambiguously confirmed that at least a fraction of SGRBs originate from binary neutron star mergers. However, whether or not a long-lived NS remnant could be formed during this merger event remains an open question (e.g., Ai et al. 2018, 2019).

The recent discovery of millisecond pulsar MSP J0740+6620, with a mass ${2.14}_{-0.09}^{+0.10}{M}_{\odot }$ (Cromartie et al. 2019), posed a strong constraint on the equation of state (EOS) of high-density matter. This mass could be used as a lower limit on the maximum NS mass and could rule out soft EOSs of NS that cannot produce such a high-mass NS. For a relatively small total mass of a binary neutron star (BNS) system, a long-lived remnant could be formed (Dai et al. 2006; Giacomazzo & Perna 2013; Zhang 2013). A rapidly spinning magnetar has been suggested as the central engine of gamma-ray bursts (GRBs; Duncan & Thompson 1992; Usov 1992; Dai & Lu 1998; Zhang & Mészáros 2001; Dai et al. 2006). In the case of SGRBs, a long-lived magnetar can help to interpret several interesting X-ray activities following the GRBs, such as X-ray plateaus (Dai & Lu 1998; Zhang & Mészáros 2001; Rowlinson et al. 2013), extended emission (Metzger et al. 2008), and X-ray flares (Dai et al. 2006). Ciolfi et al. (2019) recently performed general relativistic magnetohydrodynamic simulations of a BNS merger system up to ∼100 ms after the merge and followed the evolution of the rotational and magnetic energy of a long-lived magnetar in great detail.

Numerical simulations of NS–NS mergers indicated the typical masses of the merger ejecta M_ej ∼ 10⁻³M_⊙ to a few 10⁻²M_⊙, and the velocities of the ejecta v_ej ∼ 0.1 − 0.3c (e.g., Rezzolla et al. 2010; Hotokezaka et al. 2013; Rosswog et al. 2013; Siegel & Metzger 2017). The interaction between the subrelativistic merger ejecta with the surrounding medium would give rise to synchrotron radio emission on longer timescales of ∼a few years (Nakar & Piran 2011; Gao et al. 2013a; Hotokezaka & Piran 2015). If the merger remnant is a rapidly rotating magnetar, richer electromagnetic signals are expected. These include GRB-less X-ray transients (Zhang 2013; Sun et al. 2017), magnetar-boosted kilonova-like events known as "merger-novae" (Yu et al. 2013; Metzger & Piro 2014; Gao et al. 2015a, 2017), and the brighter forward and reverse shock emission from the interaction between the engine-powered ejecta and the surrounding medium (Gao et al. 2013a; Wang et al. 2015; Liu et al. 2016). Recently, a GRB-less X-ray transient CDF-S XT2 was reported by Xue et al. (2019), which can be interpreted as originating from the internal magnetic dissipation process in an ultrarelativistic wind of a newborn magnetar (Sun et al. 2019; Xiao et al. 2019).

A magnetar would deposit a significant fraction of its rotational energy into the merger ejecta to increase its kinetic energy. Radio observations on the timescale of ∼years after the bursts provide a probe of the total kinetic energy of ejecta. Several groups tried to search for late-time radio emission following SGRBs and use the nondetection upper limits to constrain the existence of a magnetar central engine (Metzger & Bower 2014; Fong et al. 2016; Horesh et al. 2016; Klose et al. 2019). An upper limit of a few times 10⁵¹ erg of the kinetic energy was claimed for some SGRBs, which was used to argue against a magnetar engine (Fong et al. 2016; Horesh et al. 2016). However, there is a high level of degeneracy between the kinetic energy and other model parameters. For example, in these calculations, large values of the shock microscopic parameters (e.g., _B = 0.1) have been adopted. In addition, some simplifications of the model have been adopted (e.g., in Metzger & Bower 2014 and Fong et al. 2016), which led to tighter constraints on the magnetar model.

In order to more precisely calculate the radio emission flux following SGRBs in the timescale of ∼1−10 yr after the bursts, we developed a more sophisticated model by invoking several important physical processes not fully incorporated in previous models, e.g., generic hydrodynamics, relativistic effects, and the deep Newtonian phase. We collect the late-time radio observational data of 15 SGRBs from the literature and constrain the allowed parameter space for a long-lasting NS as the BNS merger remnant using the observations. In Section 2, we describe our model in detail. In Section 3, we show the applications of our model to the observations. Our conclusions and discussion are presented in Section 4.

If the EOS of NSs is stiff enough, at least a fraction of the BNS mergers will leave behind a supramassive or even a stable NS that spins rapidly with a strong magnetic field (Dai et al. 2006; Zhang 2013; Gao et al. 2016; Piro et al. 2017; Margalit & Metzger 2019). Such a magnetar would deposit a significant fraction of its rotational energy into the merger ejecta. The kinetic energy of the merger ejecta would significantly increase. The interaction between the merger ejecta and the ambient medium produces radio emission via synchrotron radiation of relativistic electrons. Due to the additional energy injection from the long-lasting magnetar remnant, the radio brightness would be significantly enhanced (Gao et al. 2013a; Metzger & Bower 2014; Fong et al. 2016; Horesh et al. 2016).

The rotational energy of an NS formed by a BNS merger is

$$\begin{eqnarray}&&{E}_{\mathrm{rot}}=\displaystyle \frac{1}{2}I{{\rm{\Omega }}}^{2}\simeq 2\times {10}^{52}{I}_{45}{\left(\displaystyle \frac{{P}_{0}}{1\mathrm{ms}}\right)}^{-2}\mathrm{ergs},\end{eqnarray} \tag{ 1 }$$

where I is the moment of inertia of the proto-NS, and for a massive NS formed from a BNS merger, one has I₄₅ ∼ 1.5. All quantities are in c.g.s. units, and the convention Q_n = Q/10ⁿ has been adopted throughout the paper. Because the merging BNS has a high orbital angular momentum, the post-merger proto-NS would be rotating extremely rapidly, with an initial rotation period close to the centrifugal breakup limit, e.g., P₀ ∼ 1 ms. The rotation energy E_rot in Equation (1) presents a characteristic energy for the magnetar model to be tested with the radio data. Since the millisecond pulsar energy injection is essentially isotropic, the injected energy can be regarded as the isotropic equivalent energy in the ejecta–medium interaction model discussed in the rest of the paper.

Due to the dissipation of the newborn magnetar wind, a fraction of E_rot would be radiated to power early bright X-ray and optical emissions (Zhang 2013; Sun et al. 2017). It is possible that some fractions of the energy are radiated by secular gravitational waves (Fan et al. 2013; Gao et al. 2016) or fall into the black hole for a supramassive NS that collapses before fully spinning down (Gao et al. 2016). In any case, a good fraction of the rotation energy would be transferred into the merger ejecta, as ${E}_{\mathrm{inj}}=\xi {E}_{\mathrm{rot}}$, where ξ < 1 is the fraction of rotation energy that is injected into the shock. Whether or not the ejecta can be accelerated to a relativistic speed depends on E_inj and the ejecta mass M_ej. With ${E}_{\mathrm{inj}}\sim {M}_{\mathrm{ej}}{c}^{2}$, one can define a characteristic ejecta mass (Gao et al. 2013a)

$$\begin{eqnarray}&&{M}_{\mathrm{ej},{\rm{c}}}\sim 1.1\times {10}^{-2}{M}_{\odot }\xi \left(\displaystyle \frac{{E}_{\mathrm{rot}}}{2\times {10}^{52}\,\mathrm{ergs}}\right).\end{eqnarray} \tag{ 2 }$$

The ejecta lighter than ${M}_{\mathrm{ej},{\rm{c}}}$ can be accelerated to a relativistic speed. For such a case, some relativistic effects should be taken into account.

In order to calculate the blast-wave dynamics in both relativistic and nonrelativistic (Newtonian) phases, we use the generic dynamical model proposed by Huang et al. (1999).³ Consider the energy injection from the magnetar and deceleration of the ejecta due to interaction with the circumburst medium. The bulk Lorentz factor of the shock Γ evolves with the ejecta radius R as (Liu & Chen 2014)

$$\begin{eqnarray}&&\displaystyle \frac{d{\rm{\Gamma }}}{{dR}}=\displaystyle \frac{4\pi {R}^{2}{{nm}}_{p}}{{M}_{\mathrm{ej}}+2{\rm{\Gamma }}{M}_{\mathrm{sw}}}\left[\displaystyle \frac{{L}_{\mathrm{inj}}(t)}{{c}^{2}}\displaystyle \frac{{dt}}{{{dM}}_{\mathrm{sw}}}-({{\rm{\Gamma }}}^{2}-1)\right],\end{eqnarray} \tag{ 3 }$$

where n is the number density of the surrounding medium, m_p is the proton mass, c is the speed of light, and L_inj(t) is the injected luminosity from the magnetar. We characterize the injection luminosity as ${L}_{\mathrm{inj}}(t)=\xi {L}_{\mathrm{sd}}(t)$. Assuming that the main channel of protomagnetar energy loss is via dipole radiation, the spin-down luminosity can be written as ${L}_{\mathrm{sd}}={L}_{\mathrm{sd},0}{(1+t/{T}_{\mathrm{sd}})}^{-2}$. The characteristic spin-down luminosity ${L}_{\mathrm{sd},0}$ and timescale T_sd critically depend on the magnetic field strength of the magnetar (given a particular E_rot that is defined by the initial period P₀).

The evolution of the mass of the swept-up medium M_sw and the radius of the ejecta R are given by (Huang et al. 1999)

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{\mathrm{sw}}}{{dR}}=4\pi {R}^{2}{{nm}}_{p}\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{dR}}{{dt}}=\displaystyle \frac{\beta c}{1-\beta },\end{eqnarray} \tag{ 5 }$$

where β is the velocity of the ejecta divided by the speed of light c. Initially, the kinetic energy of the ejecta would increase because of energy injection. When the ejecta collects a mass comparable to its own, the shock begins deceleration at the characteristic timescale

$$\begin{eqnarray}&&{t}_{\mathrm{dec}}\sim 4.9\times {10}^{5}{\xi }^{-7/3}{E}_{\mathrm{rot},52}^{-7/3}{M}_{\mathrm{ej},-3}^{8/3}{n}^{1/3}{\rm{s}}.\end{eqnarray} \tag{ 6 }$$

The radio lightcurve usually peaks at the timescale t ∼ t_dec (Nakar & Piran 2011; Metzger & Bower 2014), so the radio observations at this timescale offer an important probe of the total kinetic energy in the shock.

Nakar & Piran (2011) calculated synchrotron radio emission lightcurve in the black hole scenario without energy injection from the central engine. They adopted the kinetic energy of the ejecta ${E}_{{\rm{k}}}\sim {10}^{49}\mbox{--}{10}^{50}$ erg. The velocity of the ejecta is nonrelativistic. The shock is in a free-coasting phase early on and enters the subsequent Sedov–Taylor self-similar evolution later. Metzger & Bower (2014) and Fong et al. (2016) used a similar method to calculate the dynamical evolution for the case of a magnetar. They used the rotation energy of the magnetar E_rot instead of the kinetic energy of the ejecta E_k in the calculations. Using Equations (3)–(5), we can calculate the evolution of the bulk Lorentz factor of the ejecta with the generic hydrodynamics model. The numerical results are shown in the left panel of Figure 1. The blue solid line is the evolution of the ejecta bulk Lorentz factor based on the generic hydrodynamics model, and the green dotted line is the dynamical evolution adopted by Fong et al. (2016). We find that the two models obtain the same maximum bulk Lorentz factor of the ejecta ${{\rm{\Gamma }}}_{\max }\approx \xi {E}_{\mathrm{rot}}/{M}_{\mathrm{ej}}{c}^{2}$, but the deceleration timescale in our model is shorter than that of Fong et al. (2016).

Zoom In Zoom Out Reset image size

In the synchrotron blast-wave model (Sari et al. 1998), the observed spectra reflect the distribution of the shock-accelerated electrons Lorentz factor γ_e. It is usually assumed that the electron energy spectrum is a power law with slope p, i.e.,

$$\begin{eqnarray}&&{dN}\propto {\gamma }_{e}^{-p}d{\gamma }_{e},\qquad {\gamma }_{e}\geqslant {\gamma }_{m},\end{eqnarray} \tag{ 7 }$$

for a mildly relativistic shock p ≈ 2.1 − 2.5 (Nakar & Piran 2011). The minimum electron Lorentz factor can be obtained based on the total energy of the accelerated electron, i.e.,

$$\begin{eqnarray}&&{\gamma }_{m}-1=\displaystyle \frac{p-2}{p-1}\displaystyle \frac{{m}_{p}}{{m}_{e}}{\epsilon }_{e}({\rm{\Gamma }}-1),\end{eqnarray} \tag{ 8 }$$

where _e is the fraction of the total internal energy of the shocked medium carried by electrons.

It is common to relate the magnetic field energy density (B²/8π) and the internal energy of the post-shocked medium (U') with a shock microphysics parameter _B. Based on the relativistic shock jump conditions, the internal energy density of the shocked medium can be written in the form of $U^{\prime} =(4{\rm{\Gamma }}+3)({\rm{\Gamma }}-1){{nm}}_{p}{c}^{2}$. Under this assumption, the magnetic field strength in the shock can be estimated by (Sari et al. 1998)

$$\begin{eqnarray}&&B=\sqrt{8\pi {\epsilon }_{B}(4{\rm{\Gamma }}+3)({\rm{\Gamma }}-1){{nm}}_{p}{c}^{2}}.\end{eqnarray} \tag{ 9 }$$

The late-time radio spectrum produced by the shock is irrelevant to the cooling frequency ν_c. The spectrum is determined by two characteristic frequencies, one is the typical synchrotron of electrons ν_m with the minimum electron Lorentz factor γ_m, i.e.,

$$\begin{eqnarray}&&{\nu }_{m}\simeq \displaystyle \frac{3}{4\pi }{\rm{\Gamma }}{\gamma }_{m}^{2}\displaystyle \frac{{eB}}{{m}_{e}c}.\end{eqnarray} \tag{ 10 }$$

The factor of Γ is introduced to transfer the shock comoving frame to the frame of the observer. The other one is the synchrotron self-absorption frequency ν_a, which can be estimated by requiring that the optical depth be equal to unity. In the case we are interested in, ν_a can be expressed as (Zhang 2018)

$$\begin{eqnarray}{\nu }_{a}=\left\{\begin{array}{ll}{\left(\displaystyle \frac{{{ \mathcal C }}_{1}{enR}}{B{\gamma }_{m}^{5}}\right)}^{3/5}{\nu }_{m}, & {\nu }_{a}\lt {\nu }_{m},\\ {\left(\displaystyle \frac{{{ \mathcal C }}_{2}{enR}}{B{\gamma }_{m}^{5}}\right)}^{2/(p+4)}{\nu }_{m}, & {\nu }_{a}\geqslant {\nu }_{m}\end{array}\right.,\end{eqnarray} \tag{ 11 }$$

where the coefficients depend on the electron power-law index p, i.e.,

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal C }}_{1}=\displaystyle \frac{16\pi \sqrt[3]{4}}{3{\rm{\Gamma }}\left(\tfrac{1}{3}\right)}\displaystyle \frac{p+2}{3p+2},\quad \mathrm{and}\\ {{ \mathcal C }}_{2}=\displaystyle \frac{{2}^{3+\tfrac{p}{2}}}{3\sqrt{3}}{\rm{\Gamma }}\left(\displaystyle \frac{1}{6}+\displaystyle \frac{p}{4}\right){\rm{\Gamma }}\left(\displaystyle \frac{11}{6}+\displaystyle \frac{p}{4}\right).\end{array}\end{eqnarray} \tag{ 12 }$$

The peak specific synchrotron emission power of a single electron in the observer frame can be expressed as (Sari et al. 1998)

$$\begin{eqnarray}&&{P}_{\nu ,\max }=\displaystyle \frac{{m}_{e}{c}^{2}{\sigma }_{T}}{3e}{\rm{\Gamma }}B,\end{eqnarray} \tag{ 13 }$$

which is independent of the electron Lorentz factor γ_e. The total number of swept-up electrons in the post-shock gas is N_e = 4π R³n/3. The observed peak flux at the luminosity distance D_L can be written as

$$\begin{eqnarray}&&{F}_{\nu ,\max }=(1+z)\displaystyle \frac{{N}_{e}{P}_{\nu ,\max }}{4\pi {D}_{L}^{2}}.\end{eqnarray} \tag{ 14 }$$

The radio-band synchrotron spectrum from the shock is governed by the relative orderings between ν_a and ν_m. There are two possible types of radio spectra (see Gao et al. 2013b; Piran et al. 2013, their Figure 4): for ν_a < ν_m, the observed flux at an observational frequency ν_obs is given by

$$\begin{eqnarray}{F}_{\nu }={F}_{\nu ,\max }\left\{\begin{array}{ll}{\left(\displaystyle \frac{{\nu }_{a}}{{\nu }_{m}}\right)}^{\tfrac{1}{3}}{\left(\displaystyle \frac{{\nu }_{\mathrm{obs}}}{{\nu }_{a}}\right)}^{2}, & {\nu }_{\mathrm{obs}}\leqslant {\nu }_{a}\lt {\nu }_{m},\\ {\left(\displaystyle \frac{{\nu }_{\mathrm{obs}}}{{\nu }_{m}}\right)}^{\tfrac{1}{3}}, & {\nu }_{a}\lt {\nu }_{\mathrm{obs}}\leqslant {\nu }_{m},\\ {\left(\displaystyle \frac{{\nu }_{\mathrm{obs}}}{{\nu }_{m}}\right)}^{-\tfrac{p-1}{2}}, & {\nu }_{a}\lt {\nu }_{m}\lt {\nu }_{\mathrm{obs},}\end{array}\right.\end{eqnarray} \tag{ 15 }$$

and for ν_a > ν_m, the observed flux F_ν is

$$\begin{eqnarray}{F}_{\nu }={F}_{\nu ,\max }\left\{\begin{array}{ll}{\left(\displaystyle \frac{{\nu }_{m}}{{\nu }_{a}}\right)}^{\tfrac{(p+4)}{2}}{\left(\displaystyle \frac{{\nu }_{\mathrm{obs}}}{{\nu }_{m}}\right)}^{2}, & {\nu }_{\mathrm{obs}}\leqslant {\nu }_{m}\lt {\nu }_{a},\\ {\left(\displaystyle \frac{{\nu }_{a}}{{\nu }_{m}}\right)}^{-\tfrac{(p-1)}{2}}{\left(\displaystyle \frac{{\nu }_{\mathrm{obs}}}{{\nu }_{a}}\right)}^{\tfrac{5}{2}}, & {\nu }_{m}\lt {\nu }_{\mathrm{obs}}\leqslant {\nu }_{a},\\ {\left(\displaystyle \frac{{\nu }_{\mathrm{obs}}}{{\nu }_{m}}\right)}^{-\tfrac{p-1}{2}}, & {\nu }_{m}\lt {\nu }_{a}\lt {\nu }_{\mathrm{obs}.}\end{array}\right.\end{eqnarray} \tag{ 16 }$$

As the shock wave sweeps across the ambient circumburst medium, the shock slows down to a nonrelativistic speed (i.e., Γ − 1≪ 1). The dynamics can be then described by the nonrelativistic Sedov–Taylor self-similar solution, β ∝ t^−3/5. If the minimum electron Lorentz factor still satisfies γ_m ≫ 1, the synchrotron flux in the radio band would decay as ${F}_{\nu }\propto {t}^{-3(5p-7)/10}$ (Frail et al. 2000). Once the majority of the shock-accelerated electrons are no longer highly relativistic, the blast wave would enter the so called "deep Newtonian phase" as studied by Huang & Cheng (2003). In this situation, according to the theory of Fermi acceleration in nonrelativistic shock, the electron spectrum is likely to be a power-law distribution in the momentum space rather than in the energy space (Sironi & Giannios 2013).

The deep Newtonian phase would begin at the time t_DN when γ_m − 1 ∼ 1 (Sironi & Giannios 2013), corresponding to the velocity of the shock $\beta \sim 0.22{\epsilon }_{e,-1}^{-1/2}$. This is at

$$\begin{eqnarray}&&{t}_{\mathrm{DN}}\sim 370\xi {E}_{52}^{1/3}{n}^{-1/3}\mathrm{days}.\end{eqnarray} \tag{ 17 }$$

When t > t_DN, most of the electron energy is contributed by the electrons with γ_e ∼ 2 and the electron spectrum follows a power-law distribution in the momentum space. In the deep Newtonian regime, the radio flux decays as ${F}_{\nu }\propto {t}^{-3(1+p)/10}$ (Granot et al. 2006; Sironi & Giannios 2013). This temporal index is shallower than the one derived by ignoring this effect.

The comparisons of our model with the four previous relevant works (i.e., Nakar & Piran 2011; Metzger & Bower 2014; Fong et al. 2016, and Horesh et al. 2016) are given in Table 1. In these previous papers, some of the important physical processes discussed here were not taken into account.

The numerical results are shown in Figure 1. The radio lightcurve peak time at 6 GHz calculated by our model (the blue solid line) is about one order of magnitude earlier than that of Fong et al. (2016). On the timescale of ∼1 − 10 yr after the bursts, which are the time windows for the observations, the theoretical luminosity calculated by our model is about one order of magnitude lower than the lightcurves predicted in Fong et al. (2016), and several orders of magnitude higher than the black hole case (Nakar & Piran 2011). With the detailed treatment of the "deep Newtonian" phase, at late times the decline rate predicted in our model is shallower than those presented in Fong et al. (2016) and Nakar & Piran (2011).

We collect 15 SGRBs with radio observations on timescales of ∼years from the literature (Metzger & Bower 2014; Fong et al. 2016; Horesh et al. 2016; Hajela et al. 2019). No radio source was detected in either case, and upper limits of the radio flux F_ν at the level of (8.4−510) μJy on the timescale of 189–3500 days after the bursts were obtained, which correspond to the luminosity upper limits $\nu {L}_{\nu }\sim (1.18\times {10}^{35}-1.83\times {10}^{40})\ \mathrm{erg}\,{{\rm{s}}}^{-1}$ (Table 2).

All the events in our sample have prior observations showing X-ray excess emission that could be a sign of the existence of a magnetar central engine (see column 7 of Table 2). Nine events have extended emission, and seven show an X-ray plateau. The X-ray afterglow of GRB 100117A shows both an X-ray plateau and flares. In particular, Lü et al. (2015) fitted the X-ray lightcurves of GRB 050724A and GRB 090510 with an "internal plateau" model. GRB 170817A displays an extended emission and a low-significance temporal feature in the X-ray afterglow, which is consistent with the reactivation of the central NS (Piro et al. 2019). Three events (i.e., GRB 050724, GRB 051221A, and GRB 130603B) have two radio observations on different frequencies and different times.

Table 2.  Late-time Radio Observations of SGRBs in Our Sample

GRB	z	νobs	Tobsa	Fνb	$\nu {L}_{\nu }$	X-Ray behavior	Referencec
		(GHz)	(days)	(μJy)	(1038 erg s−1)
050709	0.16	1.4	924	350	3.4	Extended emission	1
050724d	0.257	1.4	913	240	6.7	Extended emission	1, 2
		6.0	3500	22.1	2.7
051221Ad	0.547	1.4	759	210	34.7	Extended emission/Plateau	1, 2
		6.0	3350	19.5	14
051227	0.8	1.4	753	240	101.5	Extended emission	1
060313	0.75	1.4	677	510	183.4	Extended emission	1
060505	0.089	1.4	624	330	0.89	Extended emission	1
070714B	0.923	1.4	189	190	114.7	Extended emission	1
070724A	0.457	6.0	2768	19.1	9.1	Plateau	2
080905A	0.122	6.0	2363	22.2	0.52	Plateau	2
090510	0.903	6.0	2127	26.5	66	Extended emission	2
090515	0.403	6.0	2117	22.7	8.0	Plateau	2
100117A	0.915	6.0	1867	32	83	Plateau & Flares	2
101219A	0.718	6.0	1528	17.5	25	Plateau	2
130603Bd,e	0.356	3.0	619	60	8.6	Excess emission/Plateau	2, 3
		6.0	639	20.6	5.4
170817Af	0.00978	6.0	724	8.4	1.18 × 10−3	Extended emission	4

Notes.

^aT_obs is the observational time after the GRB in observer frame, the rest frame time after burst T_rest = T_obs/(1+z). ^bThe upper limit flux F_ν inferred by nondetection of late-time radio emission. ^cReferences for radio observations: (1) Metzger & Bower (2014), (2) Fong et al. (2016), (3) Horesh et al. (2016), and (4) Hajela et al. (2019). ^dGRB 050724, GRB 051221A, and GRB 130603B have twice the radio observations on different frequencies and different times. ^eGRB 130603B is a possible kilonova candidate. ^fPiro et al. (2019) reported a low-significance X-ray variability in GRB 170817A at 155 days after the merger.

Download table as:  ASCIITypeset image

The free parameters in our model include the injected energy from the magnetar E_inj, the mass of the merger ejecta M_ej, the initial spin-down luminosity ${L}_{\mathrm{sd},0}$, the density of the surrounding medium n, the power-law index of electron distribution p, and the shock microphysics parameters _B and _e. How various parameters might affect the properties of radio emission are shown in Figure 2. We find that there is a high level of degeneracy between the model parameters.

Zoom In Zoom Out Reset image size

In Figure 3, we show the constraint on the long-lived magnetar from the upper limit of GRB 080905A. Assuming the injected energy from the magnetar E_inj = 10⁵² erg, the ejecta mass M_ej = 0.01M_⊙, and _e = 0.1, we present the parameter space in the n − _B plane with the color indicating the contours of the observed flux F_ν. The lower left part of each panel is the allowed parameter space from the nondetection of radio emission from GRB 080905A. Since the theoretical luminosity predicted in our model is one order of magnitude lower than that of Fong et al. (2016) at the observational time of this event, T_rest = 5.769 yr, our model predicts a larger allowed parameter space for the magnetar model to survive.

Zoom In Zoom Out Reset image size

Due to the high degeneracy between model parameters, adopting different values of E_inj, M_ej, and _e would change the allowable parameter space for the same observational upper limit. The constraints on the parameter space with the nondetection of radio emission from GRB 060505 are shown in Figure 4. In the upper left panel, we take ${E}_{\mathrm{inj}}={10}^{52}$ erg, M_ej = 0.01M_⊙, and _e = 0.1 as the fiducial values of parameters. In each plot we vary one parameter while keeping the other parameters to the fiducial values. Afterglow modeling of GRB 060505 by Xu et al. (2009) suggested that the surrounding medium density is n ∼ 1 cm⁻³. From the constraint of the upper left panel, one can estimate the maximum magnetic field fraction ${\epsilon }_{B,\max }\approx 3.7\times {10}^{-4}$ in the magnetar scenario that is allowed to satisfy the radio upper limit of GRB 060505. In the upper right panel, we increase the injected energy by a factor of 10 to E_inj = 10⁵³ erg, reaching a tighter constraint on the allowed parameter space in the n − _B plane. In the lower left panel, we take a lower value of the ejecta mass M_ej = 10⁻³M_⊙. Compared with the fiducial parameters, this case has a slightly smaller allowed parameter space. In the lower right panel, we adopt a lower value of ${\epsilon }_{e}=0.01$.⁴ As shown in the panel (f) of Figure 2, lowering _e by one order of magnitude would lower the radio flux by about one order of magnitude. Therefore, in this case the allowed parameter space in the n − _B plane is greatly enlarged.

Zoom In Zoom Out Reset image size

Recently, Hajela et al. (2019) presented the VLA observations of GW170817 at ∼2 yr after the merger; they obtained the upper limit flux at 6 GHz as F_ν = 8.4 μJy. By modeling the thermal UV–optical–NIR kilonova (AT 2017gfo) associated with GW170817, Villar et al. (2017) constrained the total ejecta mass M_ej ∼ 0.08M_⊙ within the radioactive-power-dominated scenario. Such a high value of the ejecta mass is higher than the typical dynamical ejecta obtained by numerical-relativity simulations for BNS mergers (Shibata et al. 2017). The spin-down of the long-lived remnant NS offers additional energy to power the kilonova. Therefore, the required mass of the merger ejecta could be somewhat smaller than that required by the single radioactive power model (Ai et al. 2018; Li et al. 2018; Yu et al. 2018). By invoking the energy from the long-lived remnant NS, a relatively normal ejecta mass of M_ej = 0.03 ± 0.002M_⊙ could account for the kilonova. Hajela et al. (2019) modeled the broadband afterglow of GRB 170817A, indicating the circumburst medium density $n={2.5}_{-1.9}^{+4.1}\times {10}^{-3}$ cm⁻³. In Figure 5, we show the constraint on the parameter space for GRB 170817A, indicating that there is still a reasonably large parameter space to allow the existence of a long-lived NS with E_rot ∼ 10⁵² erg to satisfy the radio observation constraint.

Zoom In Zoom Out Reset image size

Broadband modeling of the SGRB afterglows could provide the measurements of the circumburst density and the shock microphysics parameters, which can be used as independent constraints on the allowed parameter space for the magnetar model. SGRBs prefer to occur in relatively low density environments with a median surrounding circumburst density of n ∼ 4 × 10⁻³ cm⁻³ (Fong et al. 2015). There is a narrow distribution of the _e values from the literature. About 62% of the GRBs in the sample adopted by Santana et al. (2014) have _e ∼ 0.1–0.3. It seems likely that _e does not change by much from burst to burst. However, there is a much wider range in the distribution of _B values. Santana et al. (2014) did a systematic study on the magnetic fields in GRB external shock based on a large X-ray and optical afterglow sample and found that the distribution of _B has a range of ∼10⁻⁸–10⁻³ with a median value of ∼few ×10⁻⁵. Gao et al. (2015b) found that the value of the magnetic equipartition parameter in the external shock ranges from 10⁻⁶ to 10⁻². We collect the inferred densities n for each burst from the literature; the values are listed in Table 3. The constraints on the parameter space for the other 12 SGRBs in our sample are shown in Figure 6. There is no detailed modeling of the afterglows of GRB 0051227 and GRB 090515 to constrain the densities. Assuming that a magnetar injects 10⁵² erg of the rotational energy into the surrounding medium and an ejecta mass M_ej = 0.01M_⊙, the maximum allowed values of _B are listed in Table 3. We find that the constraints on the maximum _B for GRB 060313, 070714B, 070724A, 090510, and 101219A reach the upper limit we set a prior. For all 15 SGRBs with nondetection of the radio emission, the constraints on the upper limit of _B in the magnetar scenario are consistent with the expectations from the modeling of GRB afterglows.

Zoom In Zoom Out Reset image size

Table 3.  Constraints on the Magnetic Equipartition Parameter _B

GRB	na	${\epsilon }_{B,\max }$ b	Referencec
	(cm−3)
050709	10−4 − 0.1	3.6 × 10−2	Panaitescu (2006)
050724	0.4 − 1.47	5.7 × 10−2	Fong et al. (2015)
051221A	2.4 × 10−3 − 0.5	1.5 × 10−3	Soderberg et al. (2006)
051227d	⋯	⋯	⋯
060313e	${3.3}_{-0.5}^{+1.0}\times {10}^{-3}$	0.1	Fong et al. (2015)
060505	1.0	3.7 × 10−4	Xu et al. (2009)
070714Be	${5.6}_{-1.1}^{+2.4}\times {10}^{-2}$	0.1	Fong et al. (2015)
070724A	${1.9}_{-1.6}^{+12}\times {10}^{-5}$	0.1	Fong et al. (2015)
080905A	${1.3}_{-1.2}^{+33}\times {10}^{-4}$	7.6 × 10−3	Fong et al. (2015)
090510e	${1.2}_{-1.0}^{+5.5}\times {10}^{-5}$	0.1	Fong et al. (2015)
090515d	⋯	⋯	⋯
100117Ae	${4.0}_{-1.0}^{+3.0}\times {10}^{-2}$	0.1	Fong et al. (2015)
101219Ae	${4.6}_{-4.3}^{+59}\times {10}^{-5}$	0.1	Fong et al. (2015)
130603B	${9.0}_{-3.0}^{+4.0}\times {10}^{-2}$	2.3 × 10−3	Fong et al. (2015)
170817Af	${2.5}_{-1.9}^{+4.1}\times {10}^{-3}$	8.6 × 10−4	Hajela et al. (2019)

Notes.

^aThe circumburst density n based on GRB afterglow modeling from the literature. ^bMaximum allowed _B by the observation assuming E_inj = 10⁵² erg, ${M}_{\mathrm{ej}}={10}^{-2}{M}_{\odot }$, ${L}_{\mathrm{sd},0}={10}^{48}$ erg s⁻¹, _e = 0.1, and p = 2.3, and we adopt the maximum value of n. ^cReferences for the circumburst densities. ^dNo afterglow modeling was available for GRB 051227 and GRB 090515. ^eThe constraints on the maximum allowed _B for GRB 060313, 070714B, 070724A, 090510, and 101219A reach the upper limit we set prior. ^fFor GRB 170817A, we adopt M_ej = 0.03M_⊙, while the rest of the parameters are kept fixed.

Download table as:  ASCIITypeset image

A long-lived magnetar remnant has been wildly invoked to explain the observational properties of the X-ray afterglows of SGRBs. Late-time radio observations of SGRBs provide a potential way to place a constraint on the existence of a long-lived magnetar remnant. We developed a sophisticated model to calculate the radio emission from the interaction between the merger ejecta and the circumburst medium in the magnetar scenario. Our model invokes several important physical processes, e.g., generic hydrodynamics, relativistic effects, and the deep Newtonian phase. The theoretical lightcurves predicted by our model in the timescale of ∼1–10 yr is about one order of magnitude lower than those predicted in previous oversimplified models (Metzger & Bower 2014; Fong et al. 2016), which used the nonrelativistic calculations following Nakar & Piran (2011) but with a higher kinetic energy (∼10⁵² erg) of the ejecta. Our generic dynamical model applies to both the relativistic and the nonrelativistic phases. Our calculations also extend to the deep Newtonian phase when the minimum Lorentz factor of the electrons, γ_m, drops below to unity, which results in a shallower decline rate of the lightcurve.

We collected 15 SGRB late-time radio observational data from the literature (Metzger & Bower 2014; Fong et al. 2016; Horesh et al. 2016; Hajela et al. 2019). All the events show an X-ray emission signature (e.g., X-ray plateau, extended emission, or X-ray flares) that may be interpreted as being powered by a magnetar central engine. No radio source was detected from any GRB in our sample. We derive the constraints on the maximally allowed _B in the magnetar scenario. Our results show that all the nondetections can be accommodated within the magnetar engine model with a reasonably large allowed parameter space, which also overlaps with that inferred from SGRB afterglow modeling. Considering the possibility of low values of shock microphysics parameters as inferred from GRB multiband afterglow observations and the simplified modeling by previous authors, the radio upper limits reported in previous works (Metzger & Bower 2014; Fong et al. 2016; Horesh et al. 2016) may not necessarily pose severe constraints on the existence of a long-lived magnetar remnant in these short GRBs.

More extreme parameters (e.g., E_rot ∼ 10⁵³ erg for a ∼2.2M_⊙ NS with a spin period close to 1 ms) are ruled out for some bursts. However, it is unlikely that a newborn supramassive NS can eject a kinetic energy of such an order. The newborn NS may possess a large ellipticity, which would release energy through secular gravitational waves (Gao et al. 2016; Ai et al. 2018). Strong gravitational-wave emission is expected in the post-merger phase, due to deformations of the core caused by the high magnetization (e.g., Dall'Osso et al. 2015). A long-lived remnant with a typical injected energy E_rot ∼ 10⁵² erg may be more likely. Such magnetars are generally allowed for all 15 SGRBs studied in our sample.

Future radio telescopes such as Square Kilometre Array and the Next Generation Very Large Array with their sub-μJy level sensitivity will be able to improve the current limits of the afterglows. The detection of late-time radio emission from the interaction the merger ejecta with the circumburst medium would confirm the existence of a long-lived magnetar remnant. Nondetections, on the other hand, would substantially tighten the parameter space allowed by the magnetar model and rule out the existence of such an engine in some cases.

We thank the referee for valuable suggestions. We thank Yun-Wei Yu for helpful discussions. This work is supported by the National Natural Science Foundation of China (NSFC) under grant Nos. 11722324, 11690024, 11603003, 11633001; the Strategic Priority Research Program of the Chinese Academy of Sciences, grant No. XDB23040100; and the Fundamental Research Funds for the Central Universities. L.D.L. is supported by the National Postdoctoral Program for Innovative Talents (grant No. BX20190044), China Postdoctoral Science Foundation (grant No. 2019M660515), and "LiYun" postdoctoral fellow of Beijing Normal University.