Coherent radio emission from a twisted magnetosphere after magnetar-quake

TL;DR

This paper proposes a model where twisted magnetospheres in magnetars generate coherent radio emissions through charge fluctuations, explaining observed variability and evolution in radio transient magnetars.

Contribution

It introduces a time-dependent, conal-beam radiative model linking magnetar magnetosphere untwisting to radio emission, aligning with observations of specific magnetars and pulsars.

Findings

Model reproduces observed radio flux and spectrum.

Explains variable radio pulsation behaviors.

Predicts detectable radio emission post X-ray outbursts.

Abstract

Magnetars are a class of highly magnetized, slowly rotating neutron stars, only a small fraction of which exhibit radio emission. We propose that the coherent radio curvature emission is generated by net charge fluctuations from a twist-current-carrying bundle (the j-bundle) in the scenario of magnetar-quake. Two-photon pair production is triggered, which requires a threshold voltage not too much higher than $10^9$ V in the current-carrying bundle, and which can be regarded as the `open field lines' of a magnetar. Continued untwisting of the magnetosphere maintains change fluctuations, and hence coherent radio emission, in the progressively shrinking j-bundle, which lasts for years until the radio beam is too small to be detected. The modeled peak flux of radio emission and the flat spectrum are generally consistent with the observations. We show that this time-dependent, conal-beam, radiative model can interpret the variable radio pulsation behaviors and the evolution of the X-ray hot spot of the radio transient magnetar XTE J1810$-$197 and the high-$B$ pulsar/AXP PSR J1622$-$4950. Radio emission with luminosity of $\lesssim10^{31}\,{\rm erg\,s^{-1}}$ and high-frequency oscillations are expected to be detected for a magnetar after an X-ray outburst. Differences of radio emission between magnetars and ordinary pulsars are discussed.

Equations

$$\Delta \psi =\int d\varphi =\int \frac{B_{\varphi }dl}{Br\sin \theta }.$$

$$\mathbf{j}=\frac{c}{4\pi }\nabla \times \mathbf{B}\simeq \frac{c\mathbf{B}}{4\pi r}{\sin }^2\theta \Delta \psi .$$

$${\rho }_{\mathrm{tw}}=\frac{|\mathbf{j}|}{c}=\beta {\rho }_{\mathrm{GJ}},$$

$${\rho }_{\mathrm{GJ}}=\frac{-\mathbf{\Omega }\cdot \mathbf{B}}{2\pi c}\frac{1}{[1-(r^2{\Omega }^2/c^2){\sin }^2\theta ]}.$$

$$u\equiv R{\sin }^2\theta /r.$$

$$\frac{d{u}_b}{dt}=-\frac{V}{2B{R}^2\Delta {\psi }_0/c+2V^{\prime }t},$$

$$\mu =\frac{0.1194}{1+0.595(173/\Gamma {)}^2}\frac{n_i(Ze{)}^2}{a},$$

$$\mu =0.09(Ze{)}^2{n}_i^{\frac{4}{3}}=9.6\times 10^{27}{\left(\frac{Z}{32}\right)}^{\frac{2}{3}}{\left(\frac{{\rho }_i}{{\rho }_N}\right)}^{\frac{4}{3}}\mathrm{erg}\mathrm{c}{\mathrm{m}}^{-3},$$

$$\mathbf{\xi }={\xi }_x\hat{\mathbf{x}}+{\xi }_y\hat{\mathbf{y}},\nabla \cdot \mathbf{\xi }=0,$$

$$T_{ij}=\mu \left(\frac{\partial {\xi }_i}{\partial x_j}+\frac{\partial {\xi }_j}{\partial x_i}\right),$$

$${\rho }_i\frac{{\partial }^2{\xi }_i}{\partial t^2}=\frac{\partial T_{ij}}{\partial x_j}+\frac{1}{4\pi }{\left[\left(\nabla \times \delta \mathbf{B}+\frac{{\partial }^2\mathbf{\xi }}{c^2\partial t^2}\times \mathbf{B}\right)\times \mathbf{B}\right]}_i,$$

$$-{\rho }_i{\omega }^2\mathbf{\xi }=-\mu \frac{l(l+1)}{R^2}\mathbf{\xi }+\frac{\partial }{\partial r}(\mu \frac{\partial \mathbf{\xi }}{\partial r})+\frac{B^2}{4\pi }(\frac{{\partial }^2\mathbf{\xi }}{\partial r^2}+\frac{{\omega }^2}{c^2}\mathbf{\xi }).$$

$${\Omega }_{\mathrm{osc}}^2=\frac{v_s^2(k_r^2+k_{\perp }^2)+v_{\mathrm{A}}^2{k}_r^2}{1-v_{\mathrm{A}}^2/c^2},$$

$$\mathbf{E}+\frac{\mathbf{v}(r)}{c}\times \mathbf{B}=0,$$

$${\mathbf{v}}_{\mathrm{osc}}=\left(0,\frac{1}{\sin \theta }{\partial }_{\varphi }{Y}_{lm},-{\partial }_{\theta }{Y}_{lm}\right){\Omega }_{\mathrm{osc}}A(r){\mathrm{e}}^{-i{\Omega }_{\mathrm{osc}}t},$$

$$\rho =\frac{\nabla \cdot \mathbf{E}}{4\pi }=\bar{\rho }+\delta \rho ,$$

$$\begin{array}{r}\bar{\rho }={\rho }_{\mathrm{GJ}}+\frac{\mathbf{\Omega }\cdot (\mathbf{r}\times \mathbf{j})}{c^2[1-(r^2{\Omega }^2/c^2){\sin }^2\theta ]},\\ \delta \rho =-\frac{Bl(l+1)Y_{lm}{\Omega }_{\mathrm{osc}}A(r){\mathrm{e}}^{-i{\Omega }_{\mathrm{osc}}t}}{4\pi cr[1-(r^2{\Omega }^2/c^2){\sin }^2\theta ]}.\end{array}$$

$$A(R)\simeq 4.4\times 10^2\left(\frac{E_q}{10^{41}\mathrm{erg}}\right){\left(\frac{R}{10\mathrm{km}}\right)}^{-1}{\left(\frac{M}{1.4M_{\text{sun}}}\right)}^{1/2}\mathrm{cm},$$

$${ϵ}_1{ϵ}_2(1-\cos {\beta }_1)\ge 2,$$

$$\begin{array}{r}{\lambda }_{1\gamma }=9.2\times 10^5{\left(\frac{P}{1\mathrm{s}}\right)}^{1/2}{\left(\frac{r}{10\mathrm{km}}\right)}^{1/2}{\left(\frac{\theta }{0.1}\right)}^{-1}\mathrm{cm},\\ {\lambda }_{\mathrm{sp}}=5.9\times 10^4{\left(\frac{P}{1\mathrm{s}}\right)}^{3/7}{\left(\frac{r}{10\mathrm{km}}\right)}^{3/7}{\left(\frac{\theta }{0.1}\right)}^{-6/7}\mathrm{cm},\end{array}$$

$$\begin{array}{r}N=\frac{p-1}{{\gamma }_1}\frac{\delta \rho }{e}{V}_e\\ =4.1\times 10^{13}\left(\frac{{\Omega }_{\mathrm{osc}}}{1\mathrm{kHz}}\right)\left(\frac{B}{10^{15}\mathrm{G}}\right)\left(\frac{L}{10\mathrm{cm}}\right)\\ \times \left(\frac{\Delta \varphi }{0.1}\right)\left(\frac{\theta }{0.01}\right){\left(\frac{P}{5\mathrm{s}}\right)}^{-1}{\left(\frac{V}{10^9\mathrm{V}}\right)}^{-1},\end{array}$$

$$\begin{array}{r}{F}_{\nu ,\max }=\frac{C(p)e^2}{6\pi c}\frac{N^2{\gamma }_1^4}{D^{2}L}{\left(\frac{\sin \Delta \varphi }{\Delta \varphi }\right)}^2{\left(\frac{{\nu }_{\mathrm{peak}}}{{\nu }_c}\right)}^{2/3}\\ =21.7{\left(\frac{B}{10^{15}\mathrm{G}}\right)}^2{\left(\frac{{\Omega }_{\mathrm{osc}}}{1\mathrm{kHz}}\right)}^2{\left(\frac{V}{10^9\mathrm{V}}\right)}^2{\left(\frac{L}{10\mathrm{cm}}\right)}^2\\ \times {\left(\frac{\Delta \varphi }{0.1}\right)}^2{\left(\frac{\theta }{0.01}\right)}^2{\left(\frac{P}{5\mathrm{s}}\right)}^{-2}{\left(\frac{D}{5\mathrm{kpc}}\right)}^{-2}\mathrm{mJy},\end{array}$$

$${\sin }^2\left(\frac{{\Gamma }_i}{2}\right)=\sin \zeta \sin \alpha {\sin }^2\left(\frac{\Delta \varphi }{4}\right)+{\sin }^2\left(\frac{\zeta -\alpha }{2}\right),$$

$$\tan {\theta }_i=-\frac{3}{2\tan {\Gamma }_i}\pm \sqrt{2+{\left(\frac{3}{2\tan {\Gamma }_i}\right)}^2}.$$

$${\eta }_i\simeq \frac{(1+\sin \zeta )r_i}{2\pi R_{\mathrm{LC}}},$$

$$\begin{array}{rl}{t}_{\mathrm{ev}}& =\frac{B{R}^2\Delta {\psi }_0}{cV}\\ & \simeq 15\left(\frac{\Delta {\psi }_0}{0.1}\right)\left(\frac{B}{10^{15}\mathrm{G}}\right){\left(\frac{R}{10\mathrm{km}}\right)}^2{\left(\frac{V}{10^9\mathrm{V}}\right)}^{-1}\mathrm{yr}.\end{array}$$

$$\begin{array}{r}{L}_X\sim \frac{E_{\mathrm{tw}}}{t_{\mathrm{ev}}}=2.6\times 10^{36}\left(\frac{B}{10^{15}\mathrm{G}}\right)\\ \left(\frac{R}{10\mathrm{km}}\right)\left(\frac{\Delta {\psi }_0}{0.1}\right)\left(\frac{V}{10^9\mathrm{V}}\right)\mathrm{erg}{\mathrm{s}}^{-1}.\end{array}$$

$$R_c=\frac{[r^2+(dr/d\theta ){]}^{3/2}}{|r^2+2(dr/d\theta {)}^2-r{d}^{2}r/d{\theta }^2|}=\frac{r(1+3{\cos }^2\theta {)}^{3/2}}{3\sin \theta (1+{\cos }^2\theta )}.$$

$$R_c=\frac{4r}{3\sin \theta }.$$

$$\alpha =\theta +\arccos \left(\frac{2\cos \theta }{\sqrt{1+3{\cos }^2\theta }}\right).$$

Full text

1. Introduction

Magnetars are highly magnetized and slowly rotating neutron stars (NSs), which are historically identified as two related classes, anomalous X-ray pulsars (AXPs) and soft gamma-ray repeaters (SGRs; see Kaspi & Beloborodov 2017 for a review). They exhibit dramatically variable X-ray and $\gamma$-ray emissions including short bursts, large outbursts, giant flares and quasi-periodic oscillations (QPO), which are believed to be powered by the dissipation of their enormous internal magnetic fields, typically $10^{14}-10^{16}$ G (Duncan & Thompson, 1992). Even their “persistent” emission is far from being steady. Magnetars are often accompanied by glitches which show irregular spin-down evolutions (e.g., Kaspi et al. 2000; Dib et al. 2008; Şaşmaz Muş et al. 2014). These behaviors may be related to the origin of magnetar bursts.

The leading scenario of magnetar bursts invokes quakes in the neutron star crust. Within such a scenario, when the pressure induced by the internal magnetic field exceeds a threshold stress, the magnetic energy releases from the crust into the magnetosphere (Thompson & Duncan, 2001), leading to particle acceleration and a short burst of radiation (Thompson & Duncan, 1995). It is then expected that the magnetar bursts may exhibit characteristics of self-organized criticality, as has been observed in earthquakes (Aschwanden 2011; e.g., Cheng et al. 1996; Duncan 1998; Goğüş et al. 1999).

There have been many attempts to detect radio pulsations from magnetars (e.g., Gaensler et al. 2001; Burgay et al. 2006). However, among the known 23 magnetars, only four have been identified as pulsed radio emitters, with the addition of one more high-B pulsar, PSR J1119$-$6127, which might be a magnetar (Olausen & Kaspi 2014 for a review, Archibald et al. 2016; Goğüş et al. 2016). The radio emission mechanism of magnetars seems to be different from that of ordinary pulsars. For instance, AXP XTE J1810$-$197 was a switched-on radio-transient during 2003 to 2009. The radio emission appeared following an X-ray outburst, and then decayed with the X-ray emission abating (Camilo et al., 2016). 1E 1547.0$-$5408, another AXP, was also identified as a switched-on radio-transient intermittently following its 2009 outburst (Burgay et al., 2009). Their radio emission differs from that of ordinary rotation-powered pulsars by having extremely variable flux densities, flatter spectra, and pulse profiles (e.g.,Camilo et al. 2007b; Lazaridis et al. 2008). Similar distinct characteristics are also found in other radio magnetars (Camilo et al., 2007c; Levin et al., 2010; Shannon & Johnston, 2013). The correlation between the radio and X-ray emission indicates that the radio emission is likely powered by quake-triggered currents in the magnetosphere rather than by the steady spindown power.

The coherent radio emission mechanism for radio pulsars is poorly understood due to their high brightness temperatures. That for magnetars is even more so thanks to their peculiar environment and the different observational properties from the ordinary pulsars. It was proposed that the coherent radio counterparts of the X-ray bursts from the magnetars have fluxes as high as 1 kJy, which is reminiscent of the solar type III radio bursts (Lyutikov, 2002, 2006). However, such a predicted radio flux is much higher than the observed peak fluxes for some radio transients (e.g., Camilo et al. 2016). Alternatively, Lin et al. (2015) argued that stellar oscillations can provide additional voltage in the polar cap region, making a “magnetar” re-active by crossing the pulsar radio emission death line. If the radio emission is originated from the oscillation-induced unipolar induction from the open field line regions, the number of observed radio magnetars would be small because of the narrow open field beam size of slow rotators.

In this paper, a coherent curvature radiative model of magnetar-quake-induced net charge fluctuation in the twist-current-carrying bundle is proposed to explain the radio emission of magnetars. The paper is organized as follows. An introductory description of magnetar-quake-induced, twisted, and oscillating magnetosphere is presented in Section 2. A radiative model is introduced in Section 3. In Section 4, we simulate the radio pulse profile evolution as magnetosphere untwisting and apply the model to XTE 1810$-$197 and PSR J1622$-$4950, respectively. The results are summarized in Section 5 with some discussion. Some detailed calculations are presented in the Appendix.

2. Magnetar quake-induced twisted magnetosphere

2.1. Twisted magnetosphere

The twist of magnetic field lines is strong inside the star and vanishes in the magnetosphere (Thompson et al., 2002). However, during a sudden crust quake, the magnetic field in the outer magnetosphere begins to twist up owing to the magnetic energy release in the crust (Beloborodov, 2009). Basically, the magnetic field lines are anchored to the star crust and the field geometry is determined by the motion of the footpoints. The relaxation of the twist makes a network of fractures in the magnetar crust, which leads to the motion of the footpoints, and the entire magnetosphere twists up. During this process, the ejected current flows along the field lines to the exterior of the star and comes back at other footpoints (Thompson et al., 2000).

The twisted magnetic field is supported by these emerging currents. This configuration can be assumed to be force-free since the energy release is dominated by the magnetic fields in the magnetosphere (Thompson et al., 2002). When the starquake shears outer magnetic fields, a toroidal component $B_{\phi}$ develops and the magnetosphere twists up. In spherical coordinates ($r,\,\theta,\,\phi$), the twist angle can be written as (Beloborodov, 2009)

[TABLE]

Once the whole magnetosphere is twisted, it cannot be untwisted rapidly because a self-induced electric field builds up, which initiates pair production and accelerates the pairs forming current flows. Such twist-maintained electric current is given by (Beloborodov & Thompson, 2007),

[TABLE]

The poloidal components $B_{r}$ and $B_{\theta}$ are not much different from the corresponding components for a normal dipole even when twists are strong, e.g., $\Delta\psi\sim 1$. The main difference is in the toroidal term $B_{\phi}$.

To maintain these currents, unsteady pair production in a short timescale is needed. The minimal charge density that is needed to support the twist-maintained current is

[TABLE]

where $\beta$ is a function of $r$ and $\theta$, and $\rho_{\rm GJ}$ is the Goldreich-Julian (GJ) charge density (Goldreich & Julian, 1969), i.e.,

[TABLE]

Here, the distribution of $\beta$ is plotted in Figure 1. The charge density of the twist-maintained current is a constant for each magnetic field line. In the corotating frame, the neutral condition of the magnetosphere is $(n_{+}-n_{-})e=\rho_{\rm GJ}$, where $n_{+}$ and $n_{-}$ are the density of positrons and electrons respectively, and the current is given by $\bm{j}=en_{+}\bm{v}_{+}-en_{-}\bm{v}_{-}$. In the region without pair production, the net charge is zero, and the electric current is maintained purely by the pairs flowing in the closed filed lines, i.e., $(n_{+}-n_{-})\ll n_{\pm}$ (Beloborodov, 2013). Note that the unstable net charge generation is a necessary condition for coherent emission. Define

[TABLE]

In the region where $u\sim 1$, it is unlikely to create coherent emission.

The energy of the twisted magnetic field is continually dissipated in the magnetosphere. It is worth noting that the magnetosphere can be divided into two parts: the twist-maintained current-carrying bundle with twisted field lines (the j-bundle), and the cavity with untwisted field lines (the u-bundle). The dissipation is mostly ohmic, since the current is maintained by the electric field $E_{\parallel}$ parallel to $\bm{B}$. The dissipation is significant at the transition boundary from the j-bundle to the u-bundle, where the strong currents dissipate as the magnetic field lines twist down rapidly. As a result, the u-bundle expands from the region of $u=1$ to that of $u=0$ (Beloborodov, 2009), as shown in Figure 1. The so-called cavity front, $u_{b}$, is defined as the boundary between the u-bundle and the j-bundle, expanding with an extremely high speed near $u_{b}=1$, then decelerates to $u_{b}\ll 1$. The expansion can be described as (Beloborodov, 2009)

[TABLE]

where $\Delta\psi_{0}$ is the initial twist angle, $V$ is the threshold voltage that can trigger plentiful pairs supplying the electric current, and $V^{\prime}=dV/du$ is the voltage gradient.

2.2. Oscillation of the magnetosphere

In a starquake, the sudden release of energy can create seismic waves. After a previous large quake (main-quake), smaller quakes (aftershocks) occur as the crust around the displaced fault plane adjusts itself to the effects of the main-quake. These quakes (main-quake or aftershock) likely excites significant magnetospheric oscillations. The toroidal modes, preserving the stellar shape, are fundamental modes which are pure shear deformations during stellar oscillations. Other modes also give rise to bulk compression and vertical motion, which have to do work against the much stronger degeneracy pressure or gravity. Therefore, the toroidal modes are most likely excited by starquakes, since the restoring force is due to the Coulomb forces of the crustal ions (Duncan, 1998).

The crust is the outermost $\sim$1 km layer of an NS where ions are locked into a solid lattice. It can store much elastic and magnetic energy. The density scale-height can be neglected for it is typically only a few percent of the crust thickness (Chamel & Haensel, 2008). Within the crust, the ions in the solid crust are arranged in a Coulomb lattice whose shear modulus is (Strohmayer et al., 1991)

[TABLE]

where $\Gamma=(Ze)^{2}/(akT)=173$ (Farouki & Hamaguchi, 1993), $n_{i}$ is the number density of ions, $Z$ is the atomic number and $a$ is the lattice constant. We adopt a bcc crystal lattice for the crust, with the lattice constant $a=(2/n_{i})^{1/3}$. The number density of ions is $n_{i}=\rho_{i}/(Z\mu_{e}m_{u})$, where $\rho_{i}$ is the mass density of the ions, $m_{u}$ is the atomic mass unit, and $\mu_{e}$ is the mean molecular weight per electron, which we adopt an intermediate value 2.5 here. Then, the shear modulus can be written as

[TABLE]

where $\rho_{N}=4\times 10^{11}\,\rm{cm^{-3}}$ is the neutron drip density.

We consider a spherical coordinate system with $r$ as the radial coordinate. The toroidal displacement is defined as

[TABLE]

where the $x$ and $y$ axes are orthogonal in the plane with ${\hat{r}}$ as a normal vector. Within the curst, we assume that the magnetic field $\bm{B}=B\bm{\hat{r}}$ is perpendicular to the curst and constant. A shear stress tensor for the Lagrangian toroidal displacement is given by (Landau & Lifshitz, 1970)

[TABLE]

where $i=x,y$. For an ion in the crustal lattice, one can obtain

[TABLE]

where $\delta\bm{B}=\nabla\times(\bm{\xi}\times\bm{B})$ is the perturbed magnetic field (Piro, 2005). It is assumed that the solution of equation (11) is $AY_{lm}\rm{exp}\it{[i(k_{r}r-\omega t)]}$, where $A$ is the amplitude of the displacement and $Y_{lm}$ stand for spherical harmonics. In the WKB limit, the vertical wave number are given by $\int k_{r}dr=n\pi$ (Piro, 2005). Then, equation (11) can be written as

[TABLE]

One can obtain the solution of equation (12),

[TABLE]

where $k_{\perp}^{2}=l(l+1)/R^{2}$, $v_{s}=(\mu/\rho_{i})^{1/2}$ is the speed of sound, and $v_{\rm A}=B/(4\pi\rho)^{1/2}$ is the Alvén speed. Figure 2 shows a few eigenfrequencies with different modes. These eigenfrequencies are very high so that it is difficult to detect periodic signals in the oscillations.

Because of the anchored field lines, the starquake-induced crust oscillations make the entire magnetosphere oscillating, which drives the fluctuations of the corotating space charges. The force-free condition for a corotation magnetosphere is

[TABLE]

where $\bm{v}(r)=\bm{v}_{\rm spin}+\bm{v}_{\rm osc}$ is the velocity at $r$ in the magnetosphere, $\bm{v}_{\rm spin}$ is the spin velocity, and $\bm{v}_{\rm osc}$ is the oscillation velocity which can be written as (Unno et al., 1989)

[TABLE]

where $A(r)$ is the oscillation amplitude. The net space-charge density is given by

[TABLE]

where

[TABLE]

The surface displacement amplitude for $m=0$ mode is given by (Duncan, 1998)

[TABLE]

where $E_{q}$ is the energy released in the starquake. For the $m\not=0$ modes, the distribution of strain is difficult, while the amplitude is roughly of the same order of magnitude. The amplitude at $r$ is estimated to be $A(r)\simeq A(R)r/R$, since the magnetosphere is corotating with the stellar surface.

3. Radio emission from a magnetar

3.1. Pair production

It is widely believed that pair production is the necessary condition for pulsar radio emission (e.g. Ruderman & Sutherland, 1975; Zhang et al., 2000). For an ordinary radio pulsar with $B\sim 10^{12}$ G, seed electrons are accelerated to an ultra-relativistic energy with $\gamma\sim 10^{7}$, emitting $\gamma-$rays via curvature radiation or inverse Compton scattering, which can be convert to pairs in strong magnetic fields (Ruderman & Sutherland, 1975; Daugherty & Harding, 1996; Zhang & Harding, 2000). The environment of a magnetar is also relevant for pair creation. The standard one-photon ($1\gamma$)-pair production can be suppressed due to magnetic photon splitting when $B$ becomes comparable to or exceeds $B_{\rm QED}\simeq 4.4\times 10^{13}$ G (Baring & Harding, 1998). Photon splitting, even it is a third-order process, can compete with $1\gamma$-pair production in magnetospheres, where photons are below pair creation threshold at their emission points (Harding & Lai, 2006). Strong vacuum dispersion may arise, so that all three CP-conservation photon-splitting modes, i.e., $\perp\to\parallel\parallel$, $\parallel\to\perp\parallel$, $\perp\to\perp\perp$, may operate together, in the special environment of a magnetar (Baring & Harding, 2001). Consequently, photons may split before reaching the $1\gamma$-pair creation threshold.

Nonetheless, two-photon($2\gamma$)-pair production can be an important source for pair production in a magnetar environment (Zhang, 2001). The interaction is between hard $\gamma$-rays and the copious X-ray photons from the magnetar surface. This process occurs mostly in the j-bundle where $E_{\parallel}$ exists. The threshold condition for the $2\gamma$-pair production when $\epsilon_{1}\ll\epsilon_{2}$ is (Gould & Schréder, 1967)

[TABLE]

where $\beta_{1}$ is the angle between the two interacting photons, and $\epsilon_{1}$ and $\epsilon_{2}$ are the soft and hard photon energies in units of $m_{\rm e}c^{2}$, respectively. For a typical blackbody radiation with temperature $kT_{B}$, the required $\gamma$-ray photon energy is $\epsilon_{2}\geq 2\times 10^{3}(kT_{B}/0.5\,{\rm keV})^{-1}$.

To generate such hard $\gamma$-ray photons, we consider the resonant scattering effect between an electron and an ambient X-ray photon (Beloborodov & Thompson, 2007). In the electron rest frame, the resonant scattering happens when the blue-shifted photon frequency matches the cyclotron frequency. The scattered photon acquires energy from Landau levels of electrons, and has a high energy $\gtrsim\gamma^{2}$ keV, which indicates a threshold voltage of $\gtrsim 10^{9}$ V in the j-bundle (Zhang & Qiao, 1996). At the $\gamma$-ray energy $\epsilon_{2}=2\times 10^{3}$, the attenuation length of $2\gamma$-pair production can be estimated as $\lambda_{\rm 2\gamma}=9.0\times 10^{3}\,\rm{cm}$ (Gould & Schréder, 1967). For the same energy, the attenuation length of $1\gamma$-pair production and photon splitting can be approximated as (Baring & Harding, 2001)

[TABLE]

where $P$ is the spin period. Thus, we have $\lambda_{\rm 2\gamma}<\lambda_{\rm sp}<\lambda_{\rm 1\gamma}$ (Zhang, 2001). It means that even if $1\gamma$-pair production is suppressed by photon splitting, the $2\gamma$-pair production process can still proceed to produce copious pairs. The attenuation lengths of each process, with different $\gamma$-ray photon energies, are plotted in Figure 3.

Starquake-induced oscillations can drive a charge deviation from the G-J density. This deviation triggers an electric field $E_{\parallel}$ parallel to $\bm{B}$. The longitudinal voltage along one magnetic field line controls the ohmically released power. To generate $2\gamma$-pair production, this voltage should not be much higher than $10^{9}$ V. The pair plasma is accelerated by this voltage, maintaining the current until it is dissipated via the ohmic effect. The $E_{\parallel}$ can be screened by pairs, but can grow again due to the continued crust oscillations, and the discharging repeats as the pair plasma leaves the discharging region.

3.2. Coherent radio emission

Theoretical models of NS coherent radio emission invoke one of three mechanisms: emission by bunches; a reactive instability, and a kinetic instability (Melrose, 2017). In this paper, we consider the bunch model. To generate coherent emission, it is required that photons are emitted in phase. In view of the observed GHz pulse duration, the curvature radiation time scale $T_{p}\sim 1$ ns is much shorter than that of the observed pulse emission $T_{\rm pul}$, so that there must be more than one bunch sweeping cross the line of sight (Yang & Zhang, 2018). The half-wavelength for the 1 GHz wave is 15 cm. If the bunch scale is smaller than the half-wavelength, the phase of emission radiated by each particle in the bunch would be approximately the same. The $N$ electrons in the bunch are assumed to move along nearly identical orbits, so that they act like a single macro-charge which emits a power $N^{2}$ times the power emitted by a single electron (Melrose, 2017). Consequently, the large fluctuating net charge of $n_{\rm GJ}$ can contribute to the coherent radiation (Yang & Zhang, 2018).

We assume that the oscillation-driven charges in the j-bundle generate coherent radio emission. It is unclear which modes would be excited and be dominant. Here we consider a typical mode with $l=2,m=0$ as an illustration. Basically, modes of $m>1$ are suppressed because $Y_{lm}$ has a term of $\sin^{m}\theta$ which is extremely small for $\theta\ll 1$. From equation (18), the displacement amplitude is $A(R)=2.2\times 10^{3}$ cm, if a quake energy $E_{q}=5\times 10^{41}$ erg is adopted for a typical neutron star. Note that an electric current may flow in closed field lines, which is different from the case of an ordinary radio pulsar. However, in the closed field lines with $r\ll R_{\rm LC}$, the unstable net charges are difficult to create and some significant absorption may exist because of the large value of the charge density there. Hence, the coherent radio emission is suggested to be originated from the regions far from the stellar surface. The curvature radius $R_{c}\sim R_{\rm LC}\simeq 2.4\times 10^{10}(P/5\,{\rm s})\,\rm{cm}$ is adopted at these locations. For $\theta<0.5$, from the equation (A1), the radiation central position can be estimated as $r\simeq 0.75R_{\rm LC}\sin\theta\simeq 1.8\times 10^{10}(P/5\,{\rm s})\sin\theta\,\rm{cm}$ (see Appendix).

It is also assumed that the charges obey a power-law distribution with a spectral index $p$ and energy cut-off at $\gamma_{1}$. The Lorentz factor of particles accelerated by the threshold voltage is $\gamma_{1}=eV/(m_{e}c^{2})\sim 10^{3}(V/10^{9}\rm\,V)$. Therefore, the electron number in the bunch volume $V_{e}\simeq 0.1Lr^{2}\theta\Delta\phi\sin\theta$ is

[TABLE]

where $p=2$ is adopted (see Appendix). The observed peak flux is

[TABLE]

where $D$ is the source distance, and $\nu_{\rm peak}/\nu_{c}={\rm min}[(4R_{c})/(3L\gamma_{1}^{3}),4740/(\gamma_{1}^{3}\theta^{3}),1]$ (see Appendix). This is generally consistent with the observed peak flux of a magnetar.

In this scenario, the radiation spectrum is a broken power law (Yang & Zhang, 2018). With the adopted parameters, the break frequencies in the spectrum can be calculated as $\nu_{l}=1.0\,{\rm GHz},\nu_{\phi}=1.4\,{\rm GHz},\nu_{c}=0.3\,{\rm GHz}$, and the broken power law spectrum is shown in Figure 7 in Appendix. The radio spectra of the magnetars are very flat so that the higher frequency signals can be observed. This flat spectrum is related to a power-law electron distribution with $p\approx 2$. For instance, the galactic center magnetar PSR J1745$-$2900 shows a shallow spectral index of $-0.4\pm 0.1$ at $2.54-225\,$GHz (Torne et al., 2015). As shown in Figure 7, $\nu_{\phi}>225\,$GHz is a necessary condition. The spectrum is flat in the range of $\nu_{l}<\nu<\nu_{c}$. From equation (B4), one can calculate $\Delta\alpha<3\times 10^{-4}$, if the curvature radius is roughly the light cylinder radius. Thus, the emission region is suggested to be at near the light cylinder. Also, in this scenario, the bunch scale is estimated to be 0.04 cm so that the millimeter waves can be coherent. At these bands, pulse scattering and dispersion caused by interstellar medium can be neglected, which is very helpful for pulsar detection.

3.3. Conal beam geometry

Observationally, the radio emission of magnetars has a variety of pulse profiles and often includes multiple emission pulse components. In our model, we adopt a phenomenological conal beam radiative model (see Rankin 1993, for a review). It has been proposed that the conal and core structures can be created via curvature radiation (e.g., Gil & Snakowski 1990; Gangadhara 2004). We assume a circular emission beam for each conal/core component. Let $(r_{i},\theta_{i})$ be the coordinates of the emission point for the $i$th cone. The angle $\Gamma_{i}$ between the pulsar magnetic axis and the magnetic field line tangential direction is calculated at the points of the $i$th cone by (Gil et al., 1984; Thorsett et al., 1991)

[TABLE]

where $\zeta$ is the angle between the line of sight (LOS) and the spin axis, $\alpha$ is the magnetic inclination angle, and $\Delta\phi$ is the corresponding apparent pulse width (or separation) resulting from such a geometric configuration. The relationship between $\Gamma_{i}$ and the emission point is (Gangadhara & Gupta, 2001)

[TABLE]

If one ignores the aberration and retardation effects, photons are emitted tangentially along the field lines and their frequencies are determined by the curvature radius and the Lorentz factor. The location of the emission point is determined by equation (A1). Typically, we assume that there are five emission components which consist of one core and two conal rings.

A shift of the position of conal components with respect to the core component, a.k.a. the so-called aberration and retardation effects, have been observed in many radio pulsars (see Krzeszowski et al. 2009 for a review). These effects address the bending of the radiation beam and the different paths of radiation from the conal emission regions to the observer. The aberration and retardation always play an important roles in low frequency emission, which is believed to originate from the high-altitude regions, i.e., the regions that are far from the stellar surface. The small net phase shift due to aberration and retardation for the $i$th cone is given by (Gangadhara & Gupta, 2001)

[TABLE]

where $\zeta$ is the angle between the line of sight and the spin axis. Both the aberration and retardation effects make the pulse components to appear at earlier longitudinal phases.

4. Case studies: XTE J1810$-$197 and PSR J1622$-$4950

In the following, we apply the physical and geometric radio emission model to two magnetars and interpret their radio emission.

4.1. Observational properties of XTE J1810$-$197

XTE J1810$-$197 is an AXP with spin period $P=5.54$ s and surface magnetic field $B\simeq 3\times 10^{14}$ G (Gotthelf & Halpern, 2007). An outburst occurred some time between 2002 November 17 and 2003 January 23 (Ibrahim et al., 2004), and the X-ray luminosity decays on a timescale of years (e.g., Halpern & Gotthelf 2005). This object switches on as a radio pulsar during the decay of the X-rays, and shows a hard spectrum, strong linear polarization and variable pulse profiles that are very different from ordinary radio pulsars (Halpern et al., 2005). From radio and X-ray observations, the geometry of XTE J1810$-$197 is inferred as $(\alpha,\,\zeta)=(52^{\circ},\,29^{\circ})$ (Bernardini et al., 2011), roughly consistent with the geometry we assumed in Section.3.2.

Camilo et al. (2016) proposed that there are five peaks on the radio pulse profile of XTE J1810$-$197. P1 appeared at the very beginning and disappeared last. This peak was suggested to be the component IV of the outer cone in our model. The component V is missing, perhaps because the line of sight does not sweep through the beam. The phase of P3 does not drift, indicating that it may be the core component I. P2 appeared at the late stage of evolution. Its presence in the early stage was not positively confirmed. The interval phase of P2 and P5 get smaller because the inner cone shrinks as the u-bundle expands. Also, there may be some spectral evolutions which lead to their flux evolution. These peaks may be the inner cone components II and III. However, P4 is very close to P5 and sometimes they are mixed. It is hard to identify P4 because there may be some multi-peak structures caused by noises.

4.2. Radio emission modeling of XTE J1810$-$197

We assume the initial center of the inner cone is at $r_{1}=0.1R_{\rm LC}$ (subscript $i=1$ for the inner cone and $i=2$ for the outer cone). The curvature radius at 1.4 GHz is estimated as $R_{c}\simeq 1.9\times 10^{9}(\gamma/700)^{3}$ cm. For simplicity, we assume that each emitting component has the same curvature radius. Therefore, we have the angular position of the inner cone $\theta_{1}=0.35$. It is worth noting that pair production sharply ends at the surface of $B\approx 10^{13}$ G (Beloborodov, 2013). Multipolar magnetic fields may exist near the stellar surface, leading to possible multi-hollow structures of charge formation regions. Emitting charges come from different formation regions which are thought to be ring-like, i.e., multiple annuli, leading to multi-cones (e.g., Gil & Sendyk 2000), shown in Figure 4. As the u-bundle expands, the shrinkage of annulus $a$ leads to the phases of the outer cone to creep towards the core component. The shape and the location of the annuli $b$ and $c$ are kept constant. The twist-maintained current is intensified, for the j-bundle twists up within the frame of $dV/du>0$ (see Section 4.3). From equation (16), we have $dE_{\parallel}/dt>0$ and $dV/dt>0$. Hence, the emission center for a given frequency moves toward the magnetic axis, while the emission cone shrinks as the u-bundle expands. For simplicity, the angle $\Gamma_{2}$ is adopted as the central angular position of the j-bundle, i.e., $\Gamma_{2}\sim 0.5\arcsin(u_{b})^{0.5}$. The center location of the cone (or core) is determined by equations (23) and (25). As the emission center gets closer to the magnetic axis, the width of each component does not decrease, until the emission cone is eaten by the u-bundle.

Later, when the u-bundle covers the entire outer cone, the outer boundary of the emitting region reaches the tangent point to the line of sight. Note that the area and location of each emission cone is determined by its charge formation region. As a result, components IV and V disappear, and the inner cone starts to shrink. Stages of the evolution of the emission components are shown in Figure 4. The geometry is under the condition of $\alpha=50^{\circ}$ and $\zeta=30^{\circ}$. Phase evolutions are also simulated and plotted in Figure 5. Finally, the entire magnetosphere twists down and the radiation damps.

From the radio profile of XTE J1810$-$197 in late 2008, the position of emission cone is estimated to be $r_{i}\sim 0.24R_{\rm LC}$ (Camilo et al., 2016) with $\theta_{i}\sim 0.65$. Thus, we have $u_{i}>u_{\rm LC}$. The radio emission is created on the closed field lines of the j-bundle, which is different from the case of an ordinary radio pulsar. For an ordinary radio pulsar, charged flow is difficult to form in closed field lines, whereas such global electric currents can be triggered in the closed field lines within the framework of a magnetar-quake. The emission region is suggested to be far from the stellar surface because of its high opacity. Such a high altitude gives rise to significant aberration and retardation. Hence, the j-bundle plays the role of the open field lines within the framework of ordinary radio pulsars, i.e., providing electric currents. It enlarges the beam size that increases the chance that the beam is swept by the LOS.

We have invoked core plus double conal emission components to interpret the radio pulse profiles of XTE 1810$-$197. The outer cone is related to the boundary between the twisted (j-bundle) and untwisted (u-bundle) regions, whereas the inner cone keeps constant before the outer cone is eaten by the u-bundle. Within the pulsar model, the core-double-cone structure was also interpreted within the framework of the inverse Compton scattering model (Qiao & Lin, 1998), in which radiative particles can be generated from one annulus. In this scenario, the shrinkage of the u-bundle leads to that the core and two cones shrink together, so that the core component may diminish because the ring C is separated from LOS. This seems not applicable to the observations of XTE 1810$-$197. Alternatively, the patchy beam model (Lyne & Manchester, 1988) has been applied to interpret the pulse profiles of radio pulsars. This model predicts no frequency dependence for the relative pulse phase between surpluses, which is inconsistent with the systematic variation of pulse components as observed in XTE 1810$-$197. Finally, a fan beam “patchy” model was proposed by Wang et al. (2014) to interpret radio pulsar beams. It predicts that the pulse width increases with the absolute value of the impact angle, in contrast with the trend for the conal beam model. This can be tested by future observations with a larger sample of radio magnetars.

The decay rate of the radio flux is large from the early state of the transient because of the high ohmic dissipation rate. It then decreases because the u-bundle expansion slows down. The spin-down rate decreases during this period of time (Camilo et al., 2007a) because of the growth of the untwisted region and tends to be a constant later as the outflow dissipates. This scenario broadly fits the observed spin-down behavior of XTE J1810$-$197. In addition, some continued large-amplitude day-to-day fluctuations on the flux density of XTE J1810$-$197 are found after 2007. These fluctuations may be caused by some aftershocks. The quake amplitude distribution of magnetars resembles that of earthquakes, e.g. obeying the Gutenberg$-$Richter law.

4.3. X-ray emission of XTE J1810$-$197

Magnetars always show some non-thermal components in the hard X-ray band, e.g., above 20 keV, even though they have Planck-like spectra in the softer band. The X-ray spectra of XTE J1810$-$197 after the outburst in 2003 are plotted in Figure 6 (Albano et al., 2010). The spectra can be well fitted by a three-component blackbody model. It is found that the area of the cold component, which comes from stellar cooling, is getting larger. In the case of twisted magnetosphere, after the starquake, the outflow particles maintain the j-bundle, where the accelerating electric field draws the positrons back from the upper pair formation front. These positrons fall onto the stellar surface, forming a hot spot emitting thermal X-ray photons (e.g., Harding & Muslimov 2001), which shrinks as the untwisted u-bundle expands.

Under the condition of $V(u)\approx$ const, the twist-maintained currents dissipate rapidly from the very beginning and then the decay rate decreases when $u\ll 1$. Therefore, equation (6) determines the evolution timescale of the j-bundle shrinkage. For $V(u)$ = const, one can calculate the timescale as (e.g. Beloborodov, 2009)

[TABLE]

This is generally consistent with the observed timescale of radio luminosity decay (e.g., Camilo et al. 2016). The twist angle is also a constant until it is eaten by the expanding u-bundle front. Additionally, the free energy stored in the twisted magnetosphere is $E_{\rm tw}\simeq B^{2}R^{3}(\Delta\psi)^{2}/24\sim 4\times 10^{44}(B/10^{15}\,{\rm G})^{2}(R/10\,{\rm km})^{3}(\Delta\psi_{0}/0.1)^{2}\,{\rm erg}$ (Beloborodov, 2009). One can calculate an average luminosity

[TABLE]

Consider the condition of $dV/du>0$ in the j-bundle. As discussed in Section.3.1, in principle, a decreased energy of hard $\gamma$-ray photons near the magnetic axis would imply a negative pair number gradient $dn_{\rm pair}/du$. Therefore, the hot spot heating rate decreases as the colatitude increases. One can also obtain $d\Delta\psi/dt>0$ and $dj/dt>0$ at the j-bundle (Beloborodov, 2009). This does not result in an extra net charge fluctuation because the evolution timescale is much longer than the duration of a pulse. A fraction of the released energy during the u-bundle twisting down propagates into the j-bundle, so that the j-bundle magnetic fields twist up. This may be the reason why the observed luminosity one year after the 2003 outburst is below the estimated luminosity $L_{X}$.

4.4. The case of PSR J1622$-$4950

PSR J1622$-$4950 is another radio emitting magnetar with a period of $P=4.3$ s and nearly $100\%$ linear polarization (Levin et al., 2010). It has a flat spectrum, highly variable flux density and pulse profiles during 2009 to 2011, with the X-ray flux decreasing by an order of magnitude (Anderson et al., 2012). Detectable radio emission was observed from 1999 to 2003, and from 2011 November to 2014 March (Scholz et al., 2017). The radio flux decreases (ranging over $\sim 3-80$ mJy) and finally disappeared as the entire magnetosphere is untwisted. The peak flux during the evolution is a few to several tens mJy, similar to that observed in XTE 1810$-$197, matching the calculation of equation (22).

However, the pulse profile of this object consists of two main bright peaks, so only one emission cone is needed without the need of introducing a core component. From equation (23) and (24), if one assumes $\alpha=20^{\circ}$ and $\zeta=10^{\circ}$ (e.g., Levin et al. 2012), the position of emission cone is estimated as $r_{i}\sim 0.1R_{\rm LC}$ with $\theta_{i}=0.71$. Therefore, one has $u_{i}\sim 5u_{\rm LC}$. The two components are emitted from the closed field lines of the j-bundle. These peaks tend to get closer (Scholz et al., 2017) since the emission cone shrinks while the u-bundle expands. The timescale of the pulse profile variation is similar to that of the flux density and the torque, which meets the estimation from equation (26). The variable torque of the magnetar, as a common property for magnetars following outbursts, was also observed in this object (Camilo et al., 2018).

The X-ray evolution also shares similar properties with the case of XTE J1810$-$197. However, the X-ray observations are not sufficient to indicate that the area of the hot spot is getting smaller. Also, no X-ray outburst was detected before its radio decay. From the exponential decay of the X-ray flux during 2007-2011, Anderson et al. (2012) argued that an undetected X-ray outburst occurred not long before mid-2007. Possibly, another X-ray outburst would have also occurred not long before 1999.

5. Summary and Discussion

A model of coherent curvature radiation for magnetars is proposed within the scenario of a magnetar-quake-induced, twisted, and oscillating magnetosphere. We show that radio emission originates from the twisted current-carrying bundle (the j-bundle), which is similar to the open field line region of normal pulsars (cf. Beloborodov 2009, who suggested that the radio emission is from the closed field line region). Continued oscillations are excited in this region due to the aftershocks and untwisting of the magnetosphere, with pairs generated via two-photon processes. The j-bundle shrinks with time as the untwisted region (the u-bundle) increases, so that the radio emission beam progressively becomes narrower. The radio emission profile evolves with time as a consequence of the shrinkage of the j-bundle, and disappears when the radio emission beam is small enough to escape the LOS. This model predicts a peak flux of radio emission and a flat spectrum that are generally consistent with the observations of magnetar radio emission. The shrinkage of the j-bundle is also consistent with the evolution of the X-ray hot spot of the magnetars during the radio active phase. We apply a time-dependent conal-beam radiative model to successfully explain the variable radio pulsation behaviors of XTE J1810$-$197 and PSR J1622$-$4950.

Our study shows that magnetars most likely have a different radio emission mechanism from ordinary radio pulsars. Even though the coherent mechanism is similar (bunched curvature radiation), the mechanisms to excite bunches are different. Whereas radio pulsars likely trigger bunches through unsteady pair production from the polar cap region defined by the open field lines, magnetars trigger bunches through quake-driven oscillations and continued untwisting of the magnetosphere. The twisted magnetosphere serves as an effective open field line region, which shrinks as a function of time. The shrinkage of this effective open field line region is the ultimate reason for the transient nature of magnetar radio emission.

The difference of the emission mechanism between magnetars and normal pulsars is also reflected on their pair production mechanisms. PSR J1119$-$6127, for instance, is a highly magnetized radio pulsar. Archibald et al. (2016) found that the persistent X-ray flux increased by a factor of 160 with a large glitch following the X-ray bursts. Unlike radio transient magnetars, the radio emission quenches following an X-ray burst, and reappears roughly two weeks later (Burgay et al., 2016a, b). Even the magnetic field is very high for this object, $1\gamma$-pair production may not be suppressed. The X-ray burst may have formed a fireball, making pair plasma density exceeding the G-J density by orders of magnitude. The leakage of these pairs to the polar cap region may have screened the parallel electric field and quenched the radio emission (Archibald et al., 2017). This is different from the magnetar case, because the triggering mechanism for pair production and radio emission is very different for magnetars.

We are grateful to Jiguang Lu and Yuan Shi at National Astronomical Observatories, Chinese Academy of Sciences, Hao Tong and Hongguang Wang at Guangzhou University, all of the members in the pulsar group at Peking University, and an anonymous referee for helpful comments and discussions about the model. W.Y.W. and X.L.C. acknowledge the support of MoST 2016YFE0100300, Natural Science Foundation of China (NSFC) 11473044, 11633004, 11653003, CAS QYZDJ-SSW-SLH017. R.X.X. acknowledges the support of National Key R&D Program of China (No. 2017YFA0402602), NSFC 11673002 and U1531243, and the Strategic Priority Research Program of CAS (No. XDB23010200).

Appendix A Dipolar geometry

We make some brief calculations of the pulsar geometry. For a dipole magnetic field, the curvature radius is

[TABLE]

It can be reduced to that of a simple relationship when $\theta\ll 1$, i.e.,

[TABLE]

The angle between the magnetic axis and the magnetic field is

[TABLE]

For $\theta\ll 1$, its derivation is

[TABLE]

The bunch volume can be estimated as (Yang & Zhang, 2018)

[TABLE]

where $\Delta\alpha$ and $\Delta\phi$ are the bunch opening angles.

Appendix B Coherent radio emisssion

For a single electron, the energy radiated per unit frequency interval per unit solid angle is given be (Jackson, 1998)

[TABLE]

where two polarized components of amplitude are

[TABLE]

where $\xi=\omega R_{c}(1/\gamma^{2}+\theta^{2})^{3/2}/(3c)$ and $K(\xi)$ is the modified Bessel function. For electrons with a power-law of $N_{e}=N(\gamma/\gamma_{1})^{-p}$ from $\gamma_{1}<\gamma<\gamma_{2}$, a coherent sum of amplitudes is

[TABLE]

The normalization of electron distribution is $N=(p-1)\gamma_{1}^{-1}\delta nV$. In this picture, the spectrum is a broken power-law (Yang & Zhang, 2018). The break frequencies are

[TABLE]

For the given parameters in Section.3.1, the spectrum are plotted in Figure 7. The peak flux is

[TABLE]

where $C(p)=4^{p/3}[\Gamma(2/3)\Gamma((p-1)/3)]^{2}$ and $T$ is the observation time.