The First Day in the Life of a Magnetar: Evolution of the Inclination Angle, Magnetic Dipole Moment, and Braking Index of Millisecond Magnetars during Gamma-Ray Burst Afterglows

The spin-down (Spitkovsky 2006) and alignment (Philippov et al. 2014) torques in the presence of a corotating plasma (Goldreich & Julian 1969) are

$$\begin{eqnarray}&&I\displaystyle \frac{d{\rm{\Omega }}}{{dt}}=-\displaystyle \frac{{\mu }^{2}{{\rm{\Omega }}}^{3}}{{c}^{3}}(1+{\sin }^{2}\alpha )\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&I\displaystyle \frac{d\alpha }{{dt}}=-\displaystyle \frac{{\mu }^{2}{{\rm{\Omega }}}^{2}}{{c}^{3}}\sin \alpha \cos \alpha .\end{eqnarray} \tag{ 2 }$$

The alignment in the presence of the corotating plasma is slower compared to the vacuum case (Philippov et al. 2014). In the plasma-filled magnetosphere model the braking index is given as (Arzamasskiy et al. 2015)

$$\begin{eqnarray}&&n=3+2{\left(\displaystyle \frac{\sin \alpha \cos \alpha }{1+{\sin }^{2}\alpha }\right)}^{2}.\end{eqnarray} \tag{ 3 }$$

The magnetic dipole moment of nascent magnetars could be changing in the timeframe of our analysis either because the magnetic field (see, e.g., Duncan & Thompson 1992) or the radius (Burrows & Lattimer 1986) of a millisecond nascent magnetar could change. The evolution of the magnetic field within the star is governed by a diffusive type partial differential equation and depends on complicated microscopic physics of the crust (see, e.g., Urpin & Muslimov 1992). To simplify matters we assume that this initial transient stage can be described by the ordinary differential equation

$$\begin{eqnarray}&&\dot{\mu }=-K{(\mu -{\mu }_{\infty })}^{q+1},\end{eqnarray} \tag{ 4 }$$

where K and q are constants and ${\mu }_{\infty }$ is the value the magnetic moment tends to relax. The long-term decay of the magnetic field (Duncan & Thompson 1992; Colpi et al. 2000) with a timescale of 10²–10⁴ yr is beyond the scope of this paper and cannot be constrained with the short-term (≲10 days) data we use in this work. It means that ${\mu }_{\infty }$ which we assume is constant in the timeframe we are interested in actually changes slowly. This nonlinear differential equation has the solution

$$\begin{eqnarray}&&\mu ={\mu }_{\infty }+({\mu }_{0}-{\mu }_{\infty }){\left(1+q\displaystyle \frac{t}{{t}_{{\rm{m}}}}\right)}^{-1/q}\end{eqnarray} \tag{ 5 }$$

where μ₀ is the initial magnetic moment, and

$$\begin{eqnarray}&&{t}_{{\rm{m}}}={({\mu }_{0}-{\mu }_{\infty })}^{-q}/K\end{eqnarray} \tag{ 6 }$$

is the field decay timescale. Note that at the linearity limit ($q\to 0$) the solution becomes

$$\begin{eqnarray}&&\mu ={\mu }_{\infty }+({\mu }_{0}-{\mu }_{\infty })\exp (-t/{t}_{{\rm{m}}}).\end{eqnarray} \tag{ 7 }$$

In case the magnetic dipole moment evolves simultaneously with the inclination angle, the braking index is given as (Ekşi 2017)

$$\begin{eqnarray}&&n=3+2{\left(\displaystyle \frac{\sin \alpha \cos \alpha }{1+{\sin }^{2}\alpha }\right)}^{2}+4\displaystyle \frac{{\tau }_{{\rm{c}}}}{{\tau }_{\mu }}\end{eqnarray} \tag{ 8 }$$

where ${\tau }_{{\rm{c}}}\equiv -{\rm{\Omega }}/2\dot{{\rm{\Omega }}}$ and ${\tau }_{\mu }\equiv -\mu /\dot{\mu }$. Using Equation (7) we obtain

$$\begin{eqnarray}&&{\tau }_{\mu }=\displaystyle \frac{{t}_{{\rm{m}}}}{{\mu }_{\infty }/\mu -1}\end{eqnarray} \tag{ 9 }$$

and by using Equation (1)

$$\begin{eqnarray}&&{\tau }_{{\rm{c}}}=\displaystyle \frac{{{Ic}}^{3}}{2{\mu }^{2}{{\rm{\Omega }}}^{2}(1+{\sin }^{2}\alpha )}.\end{eqnarray} \tag{ 10 }$$

According to the model employed in this work, the braking index of a magnetar tends to relax to n = 3, the canonical theoretical value for spinning-down magnetic dipoles, as the second and the third terms in Equation (8) vanish in time. By this we do not imply that mature magnetars will have an index of n = 3 as there could be other processes, such as winds (Gao et al. 2016), ambipolar diffusion, and Hall drift (Goldreich & Reisenegger 1992) leading to deviations from the simple magnetic dipole spin-down law in those cases. The processes that may dominate the spin-down following the afterglow stage is beyond the scope of this work. We assume that the isotropic equivalent X-ray luminosity of the afterglow scales with the spin-down luminosity ${L}_{\mathrm{sd}}=-I{\rm{\Omega }}\dot{{\rm{\Omega }}}$ so that

$$\begin{eqnarray}&&{L}_{{\rm{X}}}=\eta \displaystyle \frac{{\mu }^{2}{{\rm{\Omega }}}^{4}}{{c}^{3}}(1+{\sin }^{2}\alpha ),\end{eqnarray} \tag{ 11 }$$

where η is a constant. This factor not only represents how efficiently the spin-down luminosity is converted into X-rays in the observed band, but also takes care of the beaming and hence can be greater than unity.

Our parameter sets are {P₀, $\sin {\alpha }_{0}$, μ₀} for the constant magnetic dipole moment case and {P₀, $\sin {\alpha }_{0}$, μ₀, ${\mu }_{\infty }$, t_m} for the changing magnetic dipole moment case. In order to estimate the nascent magnetar parameters and their credible regions we followed the Bayesian approach (see, e.g., Sivia & Skilling 2006) with the posterior probability distribution defined as

$$\begin{eqnarray}&&\mathrm{posterior}\propto \mathrm{likelihood}\times \mathrm{prior}.\end{eqnarray} \tag{ 12 }$$

The Swift/XRT light-curve production phase has a set of criteria in order to achieve Gaussian statistics (Evans et al. 2007). Therefore, it is reasonable to assume error distribution is Gaussian. Here, we used a Gaussian ln-likelihood function¹ such that

$$\begin{eqnarray}&&\mathrm{ln}\left(\mathrm{likelihood}\right)=-\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{N}\left[\displaystyle \frac{{\left({y}_{i}-f({x}_{i},\theta )\right)}^{2}}{{\sigma }_{y,i}^{2}}+\mathrm{ln}\left(2\pi {\sigma }_{y,i}^{2}\right)\right]\end{eqnarray} \tag{ 13 }$$

where N is the number of data points, and y_i and σ_y,i represent the ith data point and its uncertainty, respectively. f(θ) is the model with a set of θ parameters corresponding to the x_i abscissa value. We used uniform prior probability in specified parameter ranges (see Table 1) defined as

$$\begin{eqnarray}\mathrm{ln}\left(\mathrm{prior}\right)=\left\{\begin{array}{ll}0, & \mathrm{if}\ {\theta }_{l}\,\lt \,\theta \,\lt \,{\theta }_{u}.\\ -\infty , & \mathrm{otherwise}.\end{array}\right.\end{eqnarray} \tag{ 14 }$$

Table 1.  Parameter Boundaries for Uniform Prior Distributions

Parameters	Lower Limit	Upper Limit
P0 (ms)	0.70	102
sin α0	10−3	0.99
μ0 (1033 G cm3)	10−3	10
${\mu }_{\infty }$ (1033 G cm3)	10−3	10
tm (day)	10−4	102

Download table as:  ASCIITypeset image

We do not have any strong prior information for the value of the parameters. Therefore, we choose to use uniform distribution. The upper and lower limits of the parameters (θ_l and θ_u values in Equation (14)) are chosen such that they cover a large enough parameter region for a nascent millisecond magnetar. The lower value period is chosen above the breakup limit (see, e.g., Cook et al. 1994). Lower values of magnetic dipole moments correspond to typical magnetic field values for normal pulsars (10¹² G). Upper values correspond to field values an order of magnitude greater than the highest inferred dipole field values of mature magnetars. The magnetic dipole moment decay timescale range is chosen to cover variability from ∼10 s to 100 days. We think this range is large enough for the short-term rapid evolution of magnetic field that we consider in our model. We increased the upper values of parameters (except the inclination angle) an order of magnitude (two orders of magnitude for the decay timescale) and found that our results are not affected.

The luminosity model, f(θ), we employed is given in Equation (11). We solved the coupled Equations, (1) and (2), described in Section 2.1 using scipy.integrate.odeint to find Ω(t) and α(t) and calculated the light curve with Equation (11). The efficiency factor η in Equation (11) is not a parameter that could be determined independent of μ₀ and P₀. Because of the underlying symmetry of the equations, ${P}_{0}\to {P}_{0}\sqrt{\eta }$ and ${\mu }_{0}\to {\mu }_{0}\sqrt{\eta }$ will result with the same light curve. The X-ray efficiency factors obtained from observations for mature pulsars range between 10⁻⁵–10⁻¹ (see, e.g., Kargaltsev et al. 2012). On the other hand, the jet correction factor, if we consider the relativistic beaming of the radiation, can reach up to ∼500 (see, e.g., Frail et al. 2001). Therefore, we do our analysis for η = 1. However, using the underlying symmetry above one can easily estimate the P₀ and μ₀ values that would give the same results for a different value of η. In Figure 3, we present the possible values of P₀ and μ₀ for different values of η and mark the values that correspond to η = 0.1 and η = 10 as well as η = 1. Although we fixed the moment of inertia to its canonical value, i.e., I = 10⁴⁵ g cm², in all calculations, P₀ and μ₀ for different I values will transform under ${P}_{0}\to {P}_{0}\sqrt{{I}_{45}}$ and ${\mu }_{0}\to \,{\mu }_{0}{I}_{45}$ where I₄₅ ≡ I/(10⁴⁵ g cm²). We present possible values of P₀ and μ₀ as I₄₅ changes for η = 1 in Figure 3 and specifically mark I₄₅ = 0.5 and I₄₅ = 3.0 values since these values represent approximate lower and upper values given by calculations based on different equations of state (see, e.g., Haensel et al. 2007).

In order to calculate the posterior probability distributions of our parameter sets, i.e., {P₀, sinα₀, μ₀} and {P₀, sinα₀, μ₀, ${\mu }_{\infty }$, t_m}, we used an affine-invariant Markov Chain Monte-Carlo Ensemble sampler, emcee (Foreman-Mackey et al. 2013, 2018). We used 500 walkers for each parameter. In the first run of the emcee we used 50,000 steps for the burn-in phase. In the main phase we adjusted the number of steps according to the integrated autocorrelation time (τ_f, Goodman & Weare 2010; Foreman-Mackey et al. 2013) which is calculated using the acor (Foreman-Mackey 2012) package. At every 100 steps we calculated τ_f and its difference from the previous step. When the step number is smaller than 100τ_f and the differences in the new and old τ_f values are smaller than 0.01, we assume that the chain has converged.² In order to be more conservative we continued sampling until either the step number had reached up to 500 times the maximum τ_f or 100,000 steps. Then we "thinned" the sample by half of the mean τ_f. In the second run of the emcee we initialized the chains around a small Gaussian ball of most probable values of the parameters determined in the first run of the emcee. The parameter values and their uncertainties are obtained from the median and standard deviation of the posterior probability distributions. Finally, we calculated the evolution of the braking index from Equations (3) and (8) for constant and changing magnetic dipole moment cases, respectively. The 1D and 2D posterior probability distributions of the parameters are plotted with the getdist (Lewis et al. 2018) package (i.e., the triangle plot).

GRBs display a bimodal distribution of duration (Kouveliotou et al. 1993): the short- (t < 2 s) and long-duration bursts (t ≳ 2 s). Short GRBs have long been considered, on theoretical grounds, to arise from the merger of two neutron stars (Blinnikov et al. 1984; Eichler et al. 1989), an idea which was spectacularly confirmed by the detection of gravitational waves (Abbott et al. 2017a) released during inspiral of two neutron stars associated with GRB 170817 (Abbott et al. 2017b). The merger of binary neutron stars may result with the formation of a rapidly rotating (P ∼ 1 ms) magnetar (Duncan & Thompson 1992) which is centrifugally supported and may collapse to a black hole upon slow-down. The origin of the long-duration GRBs are understood as the deaths of massive stars in core-collapse supernovae of Type 1c (Woosley 1993). Some long-duration GRBs, nevertheless, are believed to be associated with magnetars (see, e.g., Gompertz & Fruchter 2017). Therefore, we do not use the prompt emission duration as a selection criteria for candidate searches since both types could host magnetars.

We obtained the 0.3–10 keV unabsorbed flux values of the GRBs from the Swift/XRT GRB light curve repository³ (Evans et al. 2007, 2009) including flux values obtained from the photon counting mode and also from the windowed timing mode if available. We computed the luminosity values from flux light curves using

$$\begin{eqnarray}&&L=4\pi {d}_{{\rm{L}}}^{2}(z){F}_{{\rm{X}}}k(z).\end{eqnarray} \tag{ 15 }$$

Here, d_L(z) is the luminosity distance which depends on the redshift of the source, F_X is the unabsorbed flux values, and k is the cosmological k-correction due to cosmological expansion (Bloom et al. 2001) which is given as

$$\begin{eqnarray}&&k{(z)=(1+z)}^{({\rm{\Gamma }}-2)}.\end{eqnarray} \tag{ 16 }$$

Here Γ is the photon index of the power law. We obtained the photon index values from the Swift/XRT GRB light-curve repository. We calculated the luminosity distance using a flat ΛCDM cosmological model with cosmological parameters H₀ = 71 km s⁻¹ Mpc⁻¹ and Ω_M = 0.27 with astropy.cosmology subpackage (Price-Whelan et al. 2018). We obtained most of the redshift values of the sources from the same GRB light-curve repository. For the sources where there is no redshift information in that repository, we used published values from the literature. We list these values and their references for the candidates presented in this work in Table 2.

Exploring all the available GRBs we focused on the ones which contain a "plateau" phase. The presence of this "plateau" phase is what rekindled the millisecond magnetar model (Nousek et al. 2006; Zhang et al. 2006), so we assumed those events would be the most probable candidates for hosting a magnetar as the central engine. Yet, our model has limits, e.g., it does not involve the physics of fallback accretion that is expected to form in the core collapse (see, e.g., Colgate 1971) and debris disks expected to form in the aftermath of the merger of two neutron stars (Rosswog 2007), so we do not expect it to fit all the light curves. We thus eliminated many afterglows with complicated light curves that cannot be addressed with the physics involved here. In order to decide whether the light curve can be fitted with the model, we employed an initial nonlinear least squares fit using the lmfit package (Newville et al. 2014). Thus, the subset of the GRBs modeled here likely do not have dynamically significant disks that could apply a torque comparable with the magnetic dipole torques. Here we present the strongest candidates, GRBs 060510A, 070420, 070521, 140629A, 070508, 091018, and 161219B, for which the magnetar central engine would be favored. We note that although we did not use the duration of GRBs as a selection criterion, we find that all strong candidates belong to the long-GRB class (Barbier et al. 2006; Barthelmy et al. 2007; Guidorzi et al. 2007; Palmer et al. 2007, 2016; Parsons et al. 2007; Stamatikos et al. 2007; Markwardt et al. 2009; Cummings et al. 2014; D'Ai et al. 2016).