MAGNETICALLY DRIVEN WINDS FROM DIFFERENTIALLY ROTATING NEUTRON STARS AND X-RAY AFTERGLOWS OF SHORT GAMMA-RAY BURSTS

Binary neutron star (BNS) mergers represent a leading scenario to generate the physical conditions necessary for the launch of short gamma-ray bursts (SGRBs; see, e.g., Paczynski 1986; Eichler et al. 1989; Narayan et al. 1992; Barthelmy et al. 2005; Fox et al. 2005; Gehrels et al. 2005; Rezzolla et al. 2011; Berger et al. 2013; Tanvir et al. 2013). Additionally, the inspirals and mergers of BNSs are promising sources for the detection of gravitational waves (GWs) with advanced interferometers such as LIGO and Virgo (Harry et al. 2010; Accadia et al. 2011) that have a predicted rate of 40 yr⁻¹ (Abadie et al. 2010). Coincident electromagnetic (EM) and GW observations are needed to confirm the association of SGRBs with BNS mergers and to unravel the physical processes involved. Despite observational evidence in some cases (e.g., Burrows et al. 2006; Soderberg et al. 2006; Berger 2013) and a requirement coming from modeling merger rates (Metzger & Berger 2012), it is not known whether the prompt gamma-ray emission of an SGRB is always collimated in a relativistic jet. If it was, a large fraction of SGRB events would be missed, as their jets would be beamed away from us. Therefore, understanding the potential EM signatures from BNS mergers and, in particular, the non-collimated emission (e.g., Yu et al. 2013), is essential.

The Swift satellite (Gehrels et al. 2004) has recently revealed a large subclass of SGRBs that show phases of roughly constant luminosity in their X-ray afterglow lightcurves, referred to as "extended emission" and "X-ray plateaus" (e.g., Rowlinson et al. 2013; Gompertz et al. 2014). The associated X-ray emission can last from 10 to 10⁴ s and the fluence can be comparable to or larger than that of the prompt emission. These features are taken as a signature of a long-lived central engine as proposed in Zhang & Mészáros (2001), where such emission is assumed to be powered by the magnetic spin-down of a uniformly rotating object formed in a BNS merger. In general, depending on the total mass and mass ratio of the binary, the binary-merger product (BMP) can either be a stable neutron star, a supramassive neutron star (SMNS, i.e., a star with mass above the maximum mass for nonrotating configurations $M_{_{\rm TOV}}$, but below the maximum mass for uniformly rotating configurations M_max, with $M_{\rm max}\simeq (1.15\hbox{--}1.20)\,M_{_{\rm TOV}}$; Lasota et al. 1996), a hypermassive neutron star (HMNS; i.e., a star more massive than an SMNS), or a black hole.³

In this Letter, we propose a different scenario to explain the X-ray emission soon after the merger, in which the central engine is a long-lived, differentially rotating BMP and, in particular, an HMNS. Recent observations of neutron stars with masses as large as ≃ 2 M_☉ (Demorest et al. 2010; Antoniadis et al. 2013) indicate a rather stiff equation of state, while the mass distribution in BNSs is peaked around 1.3–1.4 M_☉ (Belczynski et al. 2008). This combined evidence suggests that the BMP is almost certainly an SMNS or an HMNS.

On the basis of these considerations and using three-dimensional (3D) general-relativistic ideal magnetohydrodynamic (MHD) simulations, we investigate the EM emission of an initially axisymmetric HMNS endowed with different initial magnetic fields spanning the range of reasonable configurations. In all cases, we find a very luminous, quasi-stationary EM emission from the HMNS, which is associated with a baryon-loaded outflow, driven by magnetic winding in the stellar interior. The emission, which is compatible with the observed X-ray afterglows, is almost isotropic, though a collimated, mildly relativistic flow can be produced if the magnetic field is mainly oriented along the spin axis. The isotropic emission is present even in random-field, initial configurations, making it a robust feature of a BNS merger and an important EM counterpart to the GW signatures.

As a typical HMNS, resulting from a BNS merger, we consider an axisymmetric initial model constructed using the RNS code (Stergioulas & Friedman 1995), assuming a polytropic equation of state p = Kρ^Γ, where p is the pressure, ρ the rest-mass density, Γ = 2, and K = 2.124 × 10⁵ cm⁵ g⁻¹ s⁻². The corresponding maximum gravitational mass for a uniformly rotating (nonrotating) model is ≈2.27 (1.97) M_☉. Our initial HMNS has a mass of M = 2.43 M_☉, an equatorial radius of R_e = 11.2 km, and it is differentially rotating according to a "j-constant" law (Komatsu et al. 1989), with the differential rotation parameter A/R_e = 1.112 and central period P_c = 0.47 ms.

To study the influence of the magnetic field geometry on the EM emission, we endow this model in hydrodynamic equilibrium with three different initial magnetic field geometries, which are shown in the left panels of Figures 1–3. The first model (hereafter dip-60) employs a dipolar configuration very similar to the one in Shibata et al. (2011) and Kiuchi et al. (2012), specified by the azimuthal component of the vector potential $A_\phi =A_{0,d}\varpi ^2/(r^2+\varpi _{0,d}^2/2)^{3/2}$, where ϖ is the cylindrical radius, r² ≡ ϖ² + z², A_{0, d} tunes the overall field strength, and ϖ_{0, d} ≃ 5.3 R_e ≃ 60 km indicates the radial location of the magnetic neutral point on the equatorial plane. The second model (henceforth dip-6) is endowed with the same field geometry, but with ϖ_{0, d} ≃ 0.53 R_e ≃ 6 km. In contrast to model dip-60, where the star is threaded by a magnetic field that is artificially uniform on lengthscales ≳ R_e (cf. inset in Figure 1), for model dip-6, the neutral point and maximum magnetic field are located in the interior of the star, corresponding to a more realistic distribution of currents. Finally, the third model (henceforth rand) represents a "random" magnetic field, which we compute from a vector potential built as the linear superposition of ∼6 × 10⁴ modes with random amplitudes and phases

$$\begin{eqnarray} \boldsymbol{A}_{ijk} &=& \frac{A_{0,r}\sqrt{\gamma }}{\left(r^2+\varpi _{0,r}^2\right)^{3/2}} \sum _{\ell m n=0}^{n_k} \boldsymbol{a}_{\ell m n}\cos\, (\boldsymbol{x}_{ijk} \cdot \boldsymbol{k}_{\ell m n} + 2\pi \boldsymbol{b}_{\ell m n})\nonumber \\ &&+\boldsymbol{c}_{\ell m n}\sin (\boldsymbol{x}_{ijk} \cdot \boldsymbol{k}_{\ell m n} + 2\pi \boldsymbol{d}_{\ell m n}). \end{eqnarray} \tag{ 1 }$$

Here, i, j, k label the points on the computational grid, $\boldsymbol{x}_{ijk}$ denotes the associated position vectors, and ℓ, m, n refer to a grid in wave vector $\boldsymbol{k}$ space, defined by $\lbrace \boldsymbol{k}_{\ell m n}\rbrace \equiv [0,2\pi /\lambda _\mathrm{min},2\pi /(\lambda _\mathrm{min}+\Delta \lambda),\ldots,2\pi /(\lambda _\mathrm{min}+(n_k-1)\Delta \lambda)]^3$, where n_k = 30, λ_min ≈ 3 km is the smallest wavelength employed, and Δλ ≈ 2.2 km. Furthermore, $\boldsymbol{a}_{\ell m n},\ldots,\boldsymbol{d}_{\ell m n}$ are random numbers between 0 and 1. The constant A_{0, r} sets the maximum field strength, the factor $\sqrt{\gamma }$, with γ as the determinant of the spatial metric, helps to concentrate stronger magnetic fields in regions of larger space–time curvature, while the denominator scales the magnetic field as ∼r⁻³ for r > ϖ_{0, r} ≈ 2 R_e. We built the random-field model to resemble the actual magnetic field configuration resulting from a BNS merger (cf. Rezzolla et al. 2011). In all the above cases, A_{0, d} and A_{0, r} are adjusted to yield maximum initial field strengths of B₀ = 2 × 10¹⁴ G, with the magnetic-to-fluid pressure ratio being ≪10⁻⁵ inside the star (cf. Figure 4). The resulting initial magnetic energy $E_{_\mathrm{M}}$ differs with the magnetic field geometry, being $E_{_\mathrm{M}}\simeq (530,\,1.5,\,5.1)\times 10^{44}\,\mathrm{erg}$ for models dip-60, dip-6 and rand, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The time evolution of these initial data is performed using the publicly available Einstein Toolkit with the McLachlan spacetime-evolution code (Löffler et al. 2012), combined with the fully general-relativistic ideal-MHD code WhiskyMHD (Giacomazzo & Rezzolla 2007). The ideal-fluid equation of state p = ρ(Γ − 1) is used for the evolution, where is the specific internal energy and Γ = 2. The star is initially surrounded by a low-density atmosphere with ρ_atm ≃ 6 × 10⁻⁹ρ_c, where ρ_c ≃ 1.0 × 10¹⁵ g cm⁻³ is the central rest-mass density. We carry out the MHD computations within the modified Lorenz gauge approach (Farris et al. 2012; Giacomazzo & Perna 2013), with damping parameter ξ = 2/M. The computational grid consists of a hierarchy of seven nested boxes, extending to ≈105 R_e ≈ 1180 km in the x- and y-directions and to ≈72 R_e ≈ 817 km in the z-direction. The finest refinement level corresponds to a spatial domain of [0, 18.4] × [0, 18.4] × [0, 12.8] km and covers the HMNS at all times. The highest spatial resolution is h ≃ 140 m, but lower-resolution simulations have been performed to check convergence. To make these long-term simulations computationally affordable in 3D and at high resolution, a π/2 rotation symmetry around the z-axis and a reflection symmetry across the z = 0 plane have been employed.

Figures 1–3 show snapshots of the norm of the magnetic field and rest-mass density contours in the (x, z) plane for the three models discussed above at representative times during the evolution. Due to the differential rotation in the HMNS, strong toroidal fields are generated via magnetic winding (linearly growing at the beginning). Within a few rotational periods, the build-up of magnetic pressure (mostly associated with steep toroidal-field radial gradients) is sufficient to overcome the gravitational binding in the vicinity of the stellar surface, powering an outflow of highly magnetized low-density matter with mildly relativistic velocities (cf. also Figure 4).

For model dip-60, since the initial magnetic field is oriented along the rotation axis on lengthscales ≳ R_e (cf. Figure 1, left panel), the outflow is efficiently channeled along this axis and thus significantly collimated, though this is more the result of the ordered initial magnetic field configuration than of genuine self-confining processes. During the evolution, the non-collimated part of the outflow at lower latitudes tends to push the collimated part toward the axis, effectively shrinking the inner opening angle of this jet-like structure (cf. Figure 1). This triggers irregularities in the collimated outflow that manifest themselves as small bubbles of material with low density and magnetization being released along the vertical axis, possibly pointing to a kink instability (Kiuchi et al. 2012).

The initial magnetic field geometry of dip-6 has a much smaller curvature scale and is less effective at funneling the outflow; the resulting collimation is thus much less pronounced than for dip-60. Furthermore, while model dip-60 shows a clear distinction between the collimated and the non-collimated part of the outflow, this is no longer the case for model dip-6 (cf. Figures 1 and 2). Finally, the model with the most realistic magnetic field configuration, i.e., model rand, is characterized by an almost isotropic outflow (cf. Figure 3), though we cannot exclude the fact that some form of collimation could emerge at later times. This represents an important result of our simulations: the amount of collimation in the magnetic-driven wind from the BMP will depend sensitively on the magnetic field geometry and could be absent if the field is randomly distributed.

In all of the configurations considered, the magnetized baryon-loaded outflow has rest-mass densities ∼10⁸–10⁹ g cm⁻³ and is ejected from the star with velocities v/c ≲ 0.1, in the isotropic part, and v/c ≲ 0.3, in the collimated part.

Defining the isotropic luminosity as

$$\begin{equation} L_{_\mathrm{EM}} \equiv -\oint _{r=R_\mathrm{d}} \hspace{-11.111pt}d\Omega \,\sqrt{-g}\, (T^{^{\mathrm{EM}}})^{r}_{\phantom{r}t}\,, \end{equation} \tag{ 2 }$$

where dΩ is the solid-angle element, g is the determinant of the spacetime metric, and $T^{^{\mathrm{EM}}}_{\mu \nu }$ is the EM part of the energy-momentum tensor, all of the different initial magnetic field geometries yield high, stationary EM luminosities associated with the outflow of $L_{_\mathrm{EM}}\sim 10^{48}\hbox{--}10^{50}\,\mathrm{erg}\,\mathrm{s}^{-1}$ (cf. Figure 5, upper left panel). Model dip-60 has a luminosity that is two orders of magnitude larger than those of models dip-6 and rand, which have surprisingly similar luminosities. These differences can be understood in terms of the different magnetic energies associated with the three models (see below).

Zoom In Zoom Out Reset image size

We have also verified that when computed at sufficiently large distances, the luminosity depends only weakly on the radius R_d, chosen to compute the integral in Equation (2) (see Figure 5, lower left panel). The remaining discrepancies among different R_d choices are due to the bulk of the non-collimated part of the outflow having not yet crossed the outer detection spheres at t = 60 ms. The effective bulk speed of the outflow at larger distances (as inferred from the time delay in the onset of $L_{_\mathrm{EM}}$ at different R_d) is smaller than the ejection speed of ≲ 0.1 c, due to a nonzero atmosphere rest-mass density floor. This effect is illustrated in the lower right panel of Figure 5, where a lower floor yields a higher bulk speed (cf. black and red lines). Note, however, that the stationary luminosity is not affected by the atmosphere level even when the latter is changed by one order of magnitude. Furthermore, the lower right panel of Figure 5 reports the luminosity for three different resolutions, showing that the outflow speed is affected but the stationary levels are in reasonably good agreement.

Another important conclusion to be drawn from our simulations is that, keeping the background fluid model fixed, the EM luminosity essentially depends on the initial magnetic energy of the system. In this way, the difference of two orders of magnitude between the luminosities of model dip-60 and those of models dip-6 and rand (cf. upper left panel of Figure 5) is simply due to the corresponding difference in $E_{_\mathrm{M}}$. This is demonstrated in the upper right panel of Figure 5, where we show the luminosity of a model dip-60 with an initial magnetic field strength that is rescaled to match the initial magnetic energy of model dip-6. Once the initial magnetic energies coincide, the resulting stationary luminosities are the same and therefore independent of the magnetic field geometry and degree of collimation.

Our results lead to a simple expression for the EM luminosity from the BMP,

$$\begin{equation} L_{_\mathrm{EM}} \simeq 10^{48} \chi \left(\frac{B_0}{10^{14}\,\mathrm{G}}\right)^{2} \left(\frac{R_e}{10^{6}\,\mathrm{cm}}\right)^{3} \left(\frac{P}{10^{-4}\,\mathrm{s}}\right)^{-1} \mathrm{erg}\,\mathrm{s}^{-1}, \end{equation} \tag{ 3 }$$

where P denotes the (central) spin period and B₀ denotes the initial maximum magnetic field. Our simulations reveal the need for a fudge factor, χ, that is ∼1 for the more realistic magnetic field configurations rand and dip-6, and ∼100 for dip-60, showing that initial models with the same B₀ but different field geometries lead to luminosities that differ by orders of magnitude (cf. upper left panel of Figure 5). The role of the equatorial radius, R_e, is to define the volume in which the conversion of rotational energy into EM energy occurs, while the period, P, sets the timescale of such conversion.

Equation (3) still holds if expressed in terms of quantities referring to the stage when the outflow has become roughly stationary, by substituting χ (B₀/10¹⁴ G)² with $(\bar{B}/10^{15}\,\mathrm{G})^2$, where $\bar{B}$ is the magnetic field strength reached in the outer layers of the star (R_e and P are essentially unchanged from t = 0). The fudge factor is no longer necessary, as this new expression does not depend on the magnetic field geometry.

Two remarks should be made. First, scalings similar to Equation (3) have already been reported in the literature (e.g., Meier 1999; Shibata et al. 2011), where, however, no discussion was made on the dependence of the luminosity on the initial magnetic field geometry and thus on the importance of the factor χ. Ignoring this dependence can easily lead to over/underestimates of the luminosity by orders of magnitude. Second, the scaling with radius and period in Equation (3) are very different from those resulting from the magnetic dipole spin-down emission considered in the model of Zhang & Mészáros (2001), where the luminosity is given by $L_{_\mathrm{EM}}\propto B^2\,R_e^6\,P^{-4}$. These differences in scaling may be used to discriminate between the scenario proposed here and that by Zhang & Mészáros (2001).

Defining the conversion efficiency of EM luminosity into observed X-ray luminosity as $\eta \equiv \,L_{_\mathrm{EM}}^{\mathrm{obs}}/L_{_\mathrm{EM}}$, we can use the results of our simulations for characteristic radii and periods (R_e ∼ 10⁶ cm and P ∼ 10⁻⁴ s) to deduce that magnetic field strengths in the range $\bar{B}\sim \eta ^{-1/2}\times 10^{14}\hbox{--}10^{16}\,\mathrm{G}$ are needed to produce luminosities $L_{_\mathrm{EM}}^{\mathrm{obs}}\sim 10^{46}\hbox{--}10^{51}\,\mathrm{erg}\,\mathrm{s}^{-1}$, which characterize the extended emission and X-ray plateaus of SGRBs (e.g., Rowlinson et al. 2013; Gompertz et al. 2014). Assuming an efficiency of, e.g., η ∼ 0.01–0.1 yields magnetic field strengths $\bar{B}\sim 10^{14}\hbox{--}10^{17}\,\mathrm{G}$. These strengths are not reached in the progenitor neutron stars, but they can be built from much weaker initial magnetic fields in a number of ways: the compression of stellar cores (Giacomazzo et al. 2011), the Kelvin–Helmholtz instability developing during the merger (Price & Rosswog 2006; Baiotti et al. 2008; Anderson et al. 2008; Giacomazzo et al. 2011; Zrake & MacFadyen 2013), magnetic winding, and the magnetorotational instability (MRI; Siegel et al. 2013). In particular, the MRI might enhance the EM emission, though its saturation level is unknown and its role in determining luminosity levels remains unclear. Unfortunately, resolving the MRI in long-term global simulations is out of reach. In our simulations, only magnetic winding is at play, yet magnetic fields are amplified by at least one order of magnitude.

The timescale for the persistence of differential rotation and the survival time for the persistence of differential rotation are still very much uncertain (Lasky et al. 2014) and it is hard to estimate the timescale $\tau _{_{\mathrm{EM}}}$ over which these luminosities can be sustained. Simulations show that the timescale for the extraction of angular momentum via GWs from a bar-deformed BMP is ≲ 1 s (Baiotti et al. 2008; Bernuzzi et al. 2013; but see also Fan et al. 2013). Hence, if the BMP survives for more than ∼1 s, it must have been driven toward an almost axisymmetric equilibrium by GW emission. A back of the envelope estimate of magnetic braking leads to a timescale of ∼1–10 s for the magnetic fields considered here (Shapiro 2000). This is confirmed by the evolution of the angular velocity in our simulations, which shows that $\Omega /\dot{\Omega }\sim 10$ s in the stellar interior. Both of these timescales should be meant as lower limits and differential rotation will be removed on timescales that could be 10 times larger. However, even 100 s are not sufficient to explain X-ray luminosities lasting 10³–10⁴ s. A viable scenario is one in which the BMP is an HMNS at birth, but evolves into an SMNS while differential rotation is removed via magnetic braking and rest mass is lost via the wind (our simulations show mass-loss rates of ∼10⁻³–10⁻² M_☉ s⁻¹). In this case, the BMP can survive on much longer timescales as it needs less angular momentum to prevent gravitational collapse. Once differential rotation has been removed (or is very small), the spin-down via dipolar emission through the global magnetic field produced over these timescales would power the EM emission (Zhang & Mészáros 2001; Gao & Fan 2006). Note that a similar evolution also applies if the differentially rotating BMP is an SMNS at birth.

Considering a reference luminosity of $L_{_\mathrm{EM}}\sim \,10^{48}\,\mathrm{erg}\,\mathrm{s}^{-1}$, the timescale needed to exhaust the reservoir of rotational energy T ∼ 5 × 10⁵² erg of our BMP is readily given by $\tau _{_{\mathrm{EM}}}\lesssim \,T/L_{_\mathrm{EM}}\sim 5\times 10^4\,\mathrm{s}$. Through fitting the data reported in Table 3 of Rowlinson et al. (2013), we find a power-law correlation between the observed plateau luminosities and durations: $L_{_\mathrm{EM}}^{\mathrm{obs}}\,[\mathrm{erg}\,\mathrm{s}^{-1}]\sim 10^{52}\,(\tau _{_{\mathrm{EM}}}\,[\mathrm{s}])^{-a}$, where a = 1.36 ± 0.11. Given a luminosity of $L_{_\mathrm{EM}}=L_{_\mathrm{EM}}^{\mathrm{obs}}/\eta \sim 10^{48}\,\mathrm{erg}\,\mathrm{s}^{-1}$, the fit gives a duration of the plateau emission $\tau _{_{\mathrm{EM}}}\sim (10^4/\eta)^{1/a}\,\mathrm{s}\sim 5\times 10^3\,{\rm s}$ for η ∼ 0.1. It is reassuring that even these crude estimates provide emission timescales that are well below the upper limit of $T/L_{_\mathrm{EM}}\sim 5\times 10^4\,\mathrm{s}$, thus showing the compatibility with the observations.

We also note that the estimates made above are based on the assumption that the observed emission is essentially isotropic. If there is collimation within a solid angle Ω_coll, the duration $\tau _{_{\mathrm{EM}}}$ estimated from the fit would be reduced by a factor of ∼(Ω_coll/4π)^1/a.

Using ideal MHD simulations we have investigated the EM emission of an initially axisymmetric HMNS endowed with different initial magnetic field configurations, spanning the range of geometries expected from BNS mergers. Despite the different initial configurations, we have found that, in all cases, differential rotation in the HMNS generates a strong toroidal magnetic field and a consequent baryon-loaded outflow with bulk velocities ≲ 0.1 c. This emission is almost isotropic, though an additional collimated, mildly relativistic flow is produced if the initial magnetic field has a dominant dipole component along the spin axis.

Since the emission we observe emerges as a robust feature of a BNS merger and given the consistency of the luminosity levels and duration with the observations, we conclude that the proposed physical mechanism represents a viable explanation for the X-ray afterglows of SGRBs.

We thank F. Pannarale and F. Galeazzi for discussions, W. Kastaun for help with the visualizations, and the referee for valuable comments. R.C. is supported by the Humboldt Foundation. Additional support comes from the DFG Grant SFB/Transregio 7, from "NewCompStar," COST Action MP1304, and from HIC for FAIR. The simulations have been performed on SuperMUC at LRZ Garching and on Datura at the AEI.