Flares from merged magnetars:
their prospects as a new population of gamma-ray counterparts of binary neutron star mergers

Magnetars are a kind of neutron stars (NSs) which have extremely strong magnetic fields. The magnetar’s magnetic field can be as strong as $\sim 10^{13-15}$ G (Ferrario & Wickramasinghe, 2008; Rea & Esposito, 2011; Mereghetti et al., 2015; Turolla et al., 2015; Kaspi & Beloborodov, 2017), while that of an ordinary NS is $\sim 10^{10-12}$ G (but see “low-magnetic field” magnetars (Rea et al., 2010; Turolla & Esposito, 2013)). The typical radiation activities are believed to be powered by the huge energy reservoir in the magnetic fields of magnetars, rather than their rotational energy or gravitational energy as those in spin powered or accretion powered NSs. Such magnetar radiation activities were observed as anomalous X-ray pulsars (AXPs) and Soft Gamma-ray Repeaters (SGRs). AXPs appear to be isolated pulsars with X-ray emission, whose spin down luminosity are thought to be insufficient to power their observed luminosity (Fahlman & Gregory, 1981; Gavriil et al., 2002; Kaspi & Gavriil, 2004); SGRs are thought to be magnetars which give off bursts in gamma-ray at irregular time intervals (Golenetskii et al., 1984; Norris et al., 1991; Hartmann, 1995; Thompson & Duncan, 1995; Mereghetti, 2008). Besides, there are rare cases, where much more energetic flares are emitted from magnetars, which are referred to as “Giant Flares” (Palmer et al., 2005; Hurley et al., 2005; Minaev & Pozanenko, 2020; Zhang et al., 2020; Roberts et al., 2021; Svinkin et al., 2021).

The latter two flaring activities are believed to originate from the release of magnetic energy during occasional magnetic field reconnection in magnetars. There are various theories to explain the underlying triggers of such recombination. Following the dichotomy of Sharma et al. (2023), the first class of mechanisms attributes the trigger to crustal destructive or defective events (the “star quake” paradigm, see models e.g., Blaes et al. 1989; Thompson & Duncan 1996; Levin 2006; Beloborodov 2020; Bransgrove et al. 2020; Yuan et al. 2020); while the second attributes to the reconfiguration in the twisted magnetosphere due to magnetohydrodynamical (MHD) instabilities¹¹1Specific MHD instabilities to be expected in this scenario are tearing mode instability (Lyutikov, 2003; Komissarov et al., 2007), plasmoid instability (Ripperda et al., 2019), sausage/varicose mode and kink instability (Lasky et al., 2011) etc.. (the “solar flare” paradigm, e.g., Lyutikov 2003; Komissarov et al. 2007; Ripperda et al. 2019; Mahlmann et al. 2019).

Unlike those isolated evolved magnetars, there is a population of magnetars that were born in the remnants of binary neutron star (BNS) collisions. We refer to these magnetars as merged-magnetars, which are the focus of this manuscript. It is widely believed that, BNS collision, if the remnant is not massive enough to cause a prompt collapse into a black hole (BH), will result a massive magnetar in millisecond time scale (Duncan & Thompson, 1992; Usov, 1992; Thompson, 1994; Yi & Blackman, 1998; Blackman & Yi, 1998; Kluźniak & Ruderman, 1998; Nakamura, 1998; Spruit, 1999; Wheeler et al., 2000; Ruderman et al., 2000; Kiuchi et al., 2014). Such a magnetar inherits most of the orbital angular momentum of the progenitor NSs, therefore initially possessing a short spin period of milliseconds, much shorter than the typical seconds-long spin period of magnetars. Because of their faster spin, larger mass and younger age, it is intuitive to suspect that they are even less stable than isolated evolved magnetars, and therefore, could be more likely to emit gamma-ray flares.

In this work, we investigate the possibility that flares from merged-magnetars be identified, especially as electromagnetic wave counterparts (EMC) of gravitational wave (GW) signals of BNS mergers. The manuscript is arranged as follows: In section 2, we give semi-quantitative arguments that the flare activities should dominantly occur in the merged magnetar’s early phase. Then in the next section, the detectable event rate of merged-magnetar flares is estimated from a population of BNS mergers. In section 4, we consider their potential to be identified as EMC of GW, followed by a section where the absorption from ejected matter is taken into consideration. We discuss several relevant aspects and summarize the main findings in the last two sections.

In the “star quake” paradigm, where gamma-ray flares are triggered by crustal defective events, when the magnetar is fast spinning down due to the gigantic magnetic breaking torque, the centrifugal force decreases rapidly; see Figure 1 (a). The crust of the NS will re-adjust to its new equilibrium configuration, where the centrifugal force is counteracting against the self-gravity. In this re-configuring process, crustal defective events can be expect to be much frequent than when the NS enters the later slowly spinning-down stage (such spin-down induced star quake was first discuss by Baym & Pines 1971); in the second paradigm, where gamma-ray flares are attributed to the instabilities in the magnetosphere, the boundary of the magnetar’s magnetosphere (near the light cylinder, whose radius $R_{\rm LC}=cP$, where $c$ is the velocity of light, and $P$ is the spin period) expand fast as its spin period rapidly increase; see Figure 1 (b). The magnetic field lines near the boundary will have to re-adjust accordingly, and thus are expected to be more likely to give-off gamma-ray flares²²2unlike in those conventional “solar flare” paradigm, where the magnetic twist leaks from the NS interior into the magnetosphere, through the slow crustal deformation, here the magnetic twist fast expands along with the magnetosphere. Therefore we refer to it as the “magnetosphere instability” paradigm for distinction. In fact, from our following semi-quantitative estimation we show that, the rate of gamma-ray flares in the fast spinning-down millisecond period phase is much higher than its second period phase, so that the total flare energy released in the former phase (with the time scale of days) is comparable or higher than the latter phase (with the time scale of 10⁵ years).

In the crust crack scenarios, let us assume that each crustal destructive event which is energetic enough to trigger a magnetar gamma-ray flare has a characteristic energy $E_{\rm crack}$. Denote the change in the centrifugal force as $\Delta F_{\rm cent}$ in a small interval of time $\Delta t$. The corresponding linear deformation of the NS is $\Delta l$, which we suppose is proportional to $\Delta F_{\rm cent}$ as

where $\kappa$ is the elastic factor of the crust. The work done by the gravity is thus:

where $F_{\rm G}$ denote the gravitational force acting on the crust, which remains almost constant as the deformation is small. If we equalize the work done by gravity with that released in crust cracks, we have the following equation:

where $\Delta N_{\rm crack}$ is the number of crust cracks in $\Delta t$, where their ratio is the magnetar gamma-ray flare rate

which, according to equation (2), has the following proportional relationship:

Note that if the spinning-down torque is dominated by magnetic dipole braking³³3GW braking will only dominate over the magnetic braking when the spin period is less than 0.1 s (Usov, 1992; Blackman & Yi, 1998). In this case, GW braking will fast spin-down the magnetar to a regime where the magnetic braking takes over (Zhang & Mészáros, 2001)., then we have:

where $B_{\rm s}$ the surface magnetic field strength, we assume to be an constant. As a result, equation (3) becomes:

Now, since in the first phase, the spin period is in the order of $\sim 10$ ms, and in the second phase it is $\sim$s, the $R_{\rm B}$ in the first phase can be eight orders of magnitude larger than that in the second phase. On the other hand, the time-span of the first phase is about $10^{-8}-10^{-7}$ of that of the second phase. Therefore, the magnetar gamma-ray flare energy releasing in the first phase is the same or one order of magnitude larger than that in the second phase, as we claimed in above.

If we are in the magnetosphere instability scenarios, in a short time interval $\Delta t$, the boundary of the magnetosphere (where close and open magnetic field lines transit) expand in distance: $\Delta R_{\rm LC}$. The volume which has been swept is:

The volume times the magnetic field energy density is the energy got involved. This energy is likely to be released by processes such as magnetic field re-connection near the light cylinder. We have:

where $B_{r=R_{\rm LC}}$ is the magnetic field strength at the light cylinder. For a dipole magnetic field,

Consequently, equation (7) can be reformed into:

Therefore, the number of bursts with characteristic energy $E_{\rm{burst}}$ within $\Delta t$, i.e., the characteristic burst rate $R_{\rm{B}}$ is:

Using the similar argument as in the crust crack scenarios, we can see the ratio of the $R_{\rm B}$ between the first and second phase can be ten orders of magnitude, and thus the ratio of the corresponding total energy releasing can be $10^{3}$ in the magnetosphere instability scenarios.

Define the magnetar flare energy differential number density distribution from a single merged-magnetar as: $n(E)$, where $E$ is the energy release during the flare. The total number of flares above some certain energy limit ($E_{\rm limit}$) during the life time of the magnetar is:

and the total energy released is:

which should be less than the total energy stored in the magnetosphere.

Now the rate of bursts from all merged-magnetars in the local Universe within a sphere shell of radius from $D$ to $D+dD$, in the energy range from $E$ to $E+dE$ is:

where $R_{\rm m}$ is the merger rate density of double neutron stars whose remnants are NSs instead of prompt collapsed BHs. The above equation can be further formulated to:

where $F$ is the fluence. Since $E=4\pi D^{2}F$, we have from the above equation that:

As a result, the rate of such bursts with in a limiting distance $D_{\rm u}$ and above a limiting fluence is:

where $F_{\rm limit}$ is the fluence limit of the gamma-ray detector.

The volumetric integral in equation (12) should be limited to in local Universe, where the merger rate density can be viewed as a constant, and cosmic expansion has a negligible effect. When considering the joint observation of such flares with GW detection of the BNS mergers, the integral over the luminosity distance in equation (13) should be truncated at the BNS horizon of the GW detector.

The key quantities are $n(\tau,E)$ and $R_{\rm m}$, where we reformulate the latter:

, where $\mathcal{R}_{\rm m}$ is the merger rate density of all BNS population, and $\eta$ is the fraction of those have long-lived magnetar remnants. $\mathcal{R}_{\rm m}$ can be constrained by previous GW observation at $39-1900$ Gpc^-3s^-1 (The LIGO Scientific Collaboration et al., 2021).

We assume a power-law form of $n(E)$:

Studies (Cheng et al., 2020) found the index in broad consistency with that expected from a Self-Organized Criticality (SOC, see e.g., Bak et al. 1987; Lu & Hamilton 1991; Olami et al. 1992; Aschwanden 2011) process ($\beta=5/3$), and the factor

Now the equation (15) can be simplified to:

The total energy to be released should be limited by the magnetic energy stored in the magnetosphere:

the approximant in the above equation is valid only when $1<\beta<2$. Cheng et al. (2020) found $\beta\sim 1.66$, which meets the above mentioned conditions. If we assume that the flares in ordinary SGR and GFs follows the same energy distribution law, the energy of those flares can range more than five orders of magnitudes, with the $E_{\rm u}$ corresponding to the most energetic giant flare being $E\sim 10^{46}$ ergs. Therefore, from equation (19) we find that: $\overline{E}\sim 4\times 10^{43}E^{\beta-1}_{l,42}\,\rm ergs$, where $E_{l,42}$ is the lower energy end of the $n(E)$ in unit of $10^{42}$ ergs.

The magnetic energy stored in the magnetosphere is (Zhang et al., 2022):

where $B$ is the surface magnetic field strength of the magnetar scaled with $10^{15}$ G. If we equalize the both sides of inequality (18), we can have a rough estimation of $N_{\rm B}$ as:

Taking the $N_{\rm B}$ from above, and taking $E_{l}=10^{42}$ergs, with the expression of $f_{0}$ into equation (17), we have:

If we insert the numbers into the above equation with $\beta=5/3$, we obtain that:

where $\mathcal{R}_{\rm m}$ is in unit of yr^-1/Gpc³, $D_{u,100}$ is the distance limit in unit of 100 Mpc, $F_{\rm{limit},-8}$ is the fluence limit in unit of $10^{-8}$ ergs/cm². The detection horizon $D_{u}$ is limited by the fluence cut of a gamma-ray detector as:

when taking the $E_{u}=10^{46}$ ergs, which corresponds to a conservative estimation of the total magnetic energy stored in the magnetosphere (Zhang et al., 2022). It should be noted that, since we equalized inequality 18, the estimated occurrence rate $R_{\rm B}$ is an upper limit.

A magnetar which was born with a millisecond spin period will experience two evolutionary phases. In its first phase of millisecond period of spinning, the magnetar’s spin period rapidly slows down to seconds by the strong magnetic braking torque. In the later phase, the spin period is settled to seconds scale and evolves less rapidly. We can define a transition time between the first and the second phases as:

As argued in previous sections, the burst rate before $\tau_{\rm trans}$ overwhelms that after it. Therefore, the rate in equation (23) mostly describes those bursts occuring before $\tau_{\rm trans}$. Equivalently, those bursts to be detected is likely to following a GW chirp from BNS merger with a time lag $\tau_{\rm lag}<\tau_{\rm trans}$. On the other hand, $\tau_{\rm lag}$ should be larger than $\tau_{\rm limit}$, which is the time limit less than which, the ejecta from the BNS merger is still optically thick, thus the flares from the magnetars will be largely absorbed and the temporal structure within the flares is smeared. $\tau_{\rm limit}$ is also in time scale of days (Li & Paczyński 1998, and see discussion in following section).

Since the duration of a typical magnetar GF is $\sim 0.1$-1 s, the flares detected in this phase can exhibit significant spin modulation, which can serve as an unambiguous evidence of the existence of a merged magnetar. In this case, the $D_{\rm u}$ in equation (23) is the minimum between the gamma-ray detection horizon and the GW horizon:

The flare rate as a function of fluence limit is plotted in Figure 2, see figure caption for the detailed description of the plot. When plotting Figure 2, we calculate the rate using equation (23) with Monte Carlo samplings of $\eta$, $R_{\rm m}$ and $D_{\rm GW}$: $\eta$ is uniformly randomly sampled in log-space from 0.01 to 0.1; $R$ is sampled from a log-Gaussian distribution with 1-$\sigma$ upper and lower limits correspond to 39 and 1900 yr^-1/Gpc³; $D_{\rm GW}$ is sampled from a Gaussian random with mean 300 Mpc and a standard deviation of 40 Mpc, which corresponds to the BNS detection horizon of a GW detector network with LIGO-Virgo-KAGRA (LVK) in O4 period⁴⁴4as simulated here: https://emfollow.docs.ligo.org/userguide/capabilities.html.

During the collision of BNS, abundant material will be ejected from both the tidal tail and the disk (Bovard et al., 2017; Just et al., 2015). Actually BNS are the confirmed site for rapid neutron capture nucleosynthesis ($r$-process) (Abbott et al., 2017a, b; Cowperthwaite et al., 2017; Kasen et al., 2017), which is responsible for about half of the elements heavier than iron measured in our solar system (Burbidge et al., 1957; Sneden et al., 2008). Thus, it is expected that the BNS will be surrounded by dense $r$-process material at early time, with a total ejected mass ranging from $\sim 0.005-0.1M_{\odot}$ (Bovard et al., 2017; Radice et al., 2018; Côté et al., 2018; Just et al., 2015; Fernández et al., 2015), and is optically thick to the flare gamma-ray radiation from the center remnant due to the absorption/interaction processes of photons when propagating through the material, including Compton scattering, photoelectric absorption, pair production, etc.. In the mean time, as the $r$-process material is ejected from BNS merger with high speed, the ejecta will become optically thin at $\sim$ days after the merging event (Li & Paczyński, 1998; Korobkin et al., 2020; Wang et al., 2020b). The kilonova models of GW170817 observation suggested that such ejecta has a speed ranging from $0.1c$ to $0.3c$ on average (e.g., Kasen et al., 2017; Rosswog et al., 2018; Wollaeger et al., 2018; Watson et al., 2019), as expected from previous theoretical work (e.g., Li & Paczyński, 1998; Tanaka & Hotokezaka, 2013).

Our calculation in the previous section did not include the absorption of the surrounding BNS ejecta to the flare radiation. When such effect is included, together with a finite work energy range of the detector, it effectively replaces the limiting fluence in equations (23) and (24) with a new limiting fluence $\tilde{F}_{\rm limit}$, which is related with the original $F_{\rm limit}$ as:

Here we define an absorption factor $\xi$ to describe the effect of the surrounding ejecta in absorbing the high-energy photons from the BNS flare, which is a function of time after BNS collision ($\tau$) and is defined as:

which is the ratio between the flux after absorption (observed) at time $\tau$ and the total emitted flux from the flare, and $E_{\rm low}$ and $E_{\rm high}$ denote the energy range where a specific detector works. $f_{E}=dE_{\gamma}/dEdAdt$ is the differential energy flux emitted from the flare. We assume a spectrum shape of $f_{E}$ as a power law of index -0.2 with an exponentially-cutoff at $0.48$ MeV, i.e., $f_{E}=f_{0}E^{-0.2}\exp(-E/0.48\rm{MeV})$, and $f_{0}$ is the normalization factor with $\int^{\infty}_{0}f_{E}dE=F_{\rm total}=L_{\rm total}/4\pi D^{2}$, where $D$ is the distance of the BNS, and $L_{\rm total}$ is the luminosity of the magnetar flare. This spectrum shape is taken from that of the GF from magnetar SGR 1806-20 (Palmer et al., 2005).

For approximation, we assume a uniform spherical $r$-process ejecta distribution as in Wang et al. (2020a, b) to calculate the observed flare emission (after the ejecta absorption). Then, the emitted gamma rays after propagation through the ejecta (due to various photon interactions in the ejecta) is

where $\rho_{\rm ej}$ is the ejecta density, $\kappa(E)$ is the opacity of the ejecta to a photon with energy $E$, path-length $l$ is the distance of the photons travelling through the ejecta, in this case, $l\sim v\tau$, with $v$ to be the expanding/ejected velocity of the ejecta. Only non-interacted photons are included in the observed gamma-ray signal here; scattered photons are ignored as their effects are minimal at late times when the ejecta is nearly optically thin.

We obtain the $r$-process nuclei abundance distribution in the BNS merger ejecta using the nuclear reaction network code Portable Routines for Integrated nucleoSynthesis Modeling, or PRISM (Mumpower et al., 2018) as in Wang et al. (2020a, b). We adopt a BNS merger dynamical ejecta with robust $r$-process productions (Rosswog et al., 2013) for the baseline calculation. The opacity values for the total BNS collision ejecta are calculated using the ejecta’s composition with a mixture of the opacity values of five characteristic isotopes (Fe, Xe, Eu, Pt, and U), the detailed procedure is described in Wang et al. (2020b). The opacity values of individual $r$-process nuclei are adopted from the XCOM website⁵⁵5https://www.nist.gov/pml/xcom-photon-cross-sections-database, with photon interactions including coherent (Rayleigh) scattering, incoherent (Compton) scattering, photoelectric absorption, and pair production. For photons with energy above $\sim 10$ MeV, the main interaction with the ejecta material is pair production; while at energy range $\sim 0.5-\sim 10$ MeV, the dominant process is incoherent (Compton) scattering; for lower energy photons (below $\sim 0.5$ MeV), photoelectric process mostly takes place, manifested as edges shown at energy below $\sim 0.5$ MeV in Figure 3, especially for the blue line that reach the energy below 0.1 MeV. Such feature arises from spikes in the BNS ejecta opacity/cross section at the same photon energies. The resulting spectra after absorption by the ejecta at 5 different times (0.2 days, 0.3 days, 0.6 days, 1 day and 2 days) after BNS collision are shown in figure 3. Compared to the emitted spectrum shown as a black line, we conclude that the detection window of such bursts should be in the energy range from $\sim$1 MeV to 10 MeV, and in the time window between 0.5 and 2 days after BNS collision.

We note that in addition to the burst flare radiation, the BNS collision ejecta itself also emit gamma-ray photons through the decays of the radioactive $r$-process nuclei. The total gamma radiation rate from the $r$-process ejecta is estimated to be $\epsilon_{0}(\tau)\sim 2\times 10^{10}{\rm erg\ g}^{-1}s^{-1}(\tau/{\rm day})^{-1/3}$ (Metzger & Berger, 2012; Korobkin et al., 2020), and the $r$-process gamma-ray energy at $\sim 1$day is then $\sim 10^{41}$ erg/s for a 0.01 $M_{\odot}$ BNS merger ejecta. Thus, such signal would be small compared to the flare emissions, and the BNS ejecta spectrum shapes with nuclear decay lines (Korobkin et al., 2020; Wang et al., 2020b) are also different from the flare signal discussed here.

Then we conduct the integral in equation (28) to obtain $\xi$ as function of $\tau$ in Figure 4. The uncertainty bands are due to variations in BNS ejecta properties, including velocity, ejecta mass, and the components. Here we varied the ejecta mass between 0.005 to 0.03 $M_{\odot}$, and the velocity between 0.1c to 0.3c. To test the sensitivity of the signal to the ejecta component, we adopted the parameterized BNS outflow conditions (Just et al., 2015; Radice et al., 2018) with a range of initial electron fractions as in Wang et al. (2020b), so that the ejecta component varies from the weak $r$-process (no third peak and heavier actinides elements) to robust $r$-process (with actinides). From Figure 4, we can see that at the detection window discussed above, the corresponding absorption factor is $\xi\sim 0.5-1$ with an order of magnitude uncertainty. The burst detection rate after absorption of BNS ejecta considered is re-plotted in Figure 5. When plotting the figure, we calculate the rate with a $\xi$ randomly drawn from its corresponding range evaluated above with uniform distribution.

During the variation test, we find that the flare signal is more sensitive to the ejecta velocity and mass. Therefore, on the other hand, the detection rate obtained in the real observation could enable us to put a constraint on the BNS ejecta property.

From our above argument and calculation, we conclude that flares from merged massive magnetars can be expected as a population of gamma-ray transients, which associate GW chirp events of BNS mergers. Such a gamma-ray counterpart of GW may look like a short gamma-ray bursts (sGRB) according to its duration, but it can be found with several distinct features:

• •

  it tends to be weak in flux, and the time-lag between the burst and GW chirp is $\sim$1-2 days, rather than $\sim$ s as in sGRB;

• •

  its spectrum has a lower energy cut at $\sim$100 keV, due to absorption of the ejecta.

Besides, it may show spin modulation with a significant spin down, although the potential of significantly observing such short scale temporal structures is very challenging in reality.

Due to the absorption by the ejecta from the BNS collision, such flares are to be optimally observed in the energy range from 0.1 to 10 MeV. The estimated detection rate is increasing towards a fluence flux limit, with a power law with an index -1.5, while the rate of such bursts as GW association will also be limited by the detecting reach of GW detector networks, when the fluence limit of the HE detector is below some turn-over sensitivity. Below this turn-over fluence limit, the rate follows another power law with index -2/3. When observing with a detector of energy range 0.1-1 MeV, the turn-over flux limit is at $\sim 2\times 10^{-9}$ erg/cm², while for a detector of 1-10 MeV is at $\sim 10^{-8}$ erg/cm². Based on our evaluation from a population of BNS, a GRB monitor with energy range of $\sim 0.1-1$ MeV and a fluence limit of $\sim 2\times 10^{-8}$ erg/cm² could detect such flares as gamma-ray counterparts of GW events at a rate from 0.01 to 1 per year. To raise the detection rate of such event to a few to a few tens per year, we expect a future MeV detector working in a range from $\sim 1-10$ MeV with a fluence limit $\sim 10^{-9}$ erg/cm² for a 1s exposure time.