Axisymmetric force-free magnetosphere in the exterior of a neutron star

TL;DR

This paper models axisymmetric force-free magnetospheres of neutron stars, incorporating general relativistic effects, revealing how magnetic energy and helicity accumulate and lead to flux rope detachment and potential giant flares.

Contribution

It provides new solutions for twisted magnetospheres considering relativistic effects and explores the conditions leading to magnetic eruptions.

Findings

Relativistic effects enhance magnetic flux rope formation.

Magnetic energy and helicity have upper bounds in the magnetosphere.

Catastrophic events like giant flares may occur when helicity exceeds limits.

Abstract

Magnetar magnetospheres gradually become twisted due to shearing by footpoint motion. The axisymmetric solutions for the force-free field with a power-law current model are calculated by taking into account general relativistic effects. Here we show how the magnetic energy and helicity are accumulated along a sequence of equilibria. In a strongly twisted case, a magnetic flux rope, in which a large amount of toroidal field is confined, detaches in the vicinity of the star. It is found this kind of magnetic field line structure is easily produced due to strong relativistic effects compared with the results in flat spacetime. Although the structure is not inherent to the relativity, the confinement is possible even for a smaller power-law index. There is an upper bound on the energy and helicity stored in the magnetosphere. When the helicity is further accumulated, a catastrophic event such as a giant flare may occur.

Tables

Table 1. Table 1: Model and result. A family of the solutions is characterized by power-law index n 𝑛 n and relativistic factor M / R 𝑀 𝑅 M/R . The magnetic energy for the potential field is E 0 subscript 𝐸 0 E_{0} , and maximum energy increment Δ ​ E Δ 𝐸 \Delta E from E 0 subscript 𝐸 0 E_{0} , relative helicity H R subscript 𝐻 R H_{\rm R} and twist angle Δ ​ ϕ max Δ subscript italic-ϕ max \Delta\phi_{\rm max} are also listed.

Model	$n$	$M/R$	$E_0/\left(B_0^2{R}^3\right)$	${\left(\Delta E/E_0\right)}_{\max }$	${\left(H_R/\left(4\pi B_0^2{R}^4\right)\right)}_{\max }$	$\Delta {\varphi }_{\max }$
F5	5	0	0.33	0.34	0.74	1.05
F7	7	0	0.33	0.58	0.77	1.49
G5	5	0.25	0.71	1.32	3.25	1.38
G7	7	0.25	0.71	2.78	3.50	0.88

Equations

$$d{s}^2=-{\alpha }^{2}d{t}^2+{\alpha }^{-2}d{r}^2+r^{2}d{\theta }^2+{\varpi }^{2}d{\varphi }^2,$$

$${\alpha }^2=1-\frac{2M}{r},\varpi =r\sin \theta .$$

$$\vec{\nabla }\cdot \vec{B}=0,$$

$$\vec{\nabla }\times (\alpha \vec{B})=4\pi \alpha \vec{j}.$$

$$\vec{B}=\vec{\nabla }\times \left(\frac{G}{\varpi }{\vec{e}}_{\hat{\varphi }}\right)+\frac{S}{\alpha \varpi }{\vec{e}}_{\hat{\varphi }}=\frac{\vec{\nabla }G\times {\vec{e}}_{\hat{\varphi }}}{\varpi }+\frac{S}{\alpha \varpi }{\vec{e}}_{\hat{\varphi }},$$

$$[B_{\hat{r}},B_{\hat{\theta }},B_{\hat{\varphi }}]=\left[\frac{G{,}_{\theta }}{r\varpi },-\frac{\alpha G{,}_r}{\varpi },\frac{S}{\alpha \varpi }\right].$$

$$\vec{A}=\vec{\nabla }\times \left(\frac{F}{\varpi }{\vec{e}}_{\hat{\varphi }}\right)+\frac{G}{\varpi }{\vec{e}}_{\hat{\varphi }}=\frac{\vec{\nabla }F\times {\vec{e}}_{\hat{\varphi }}}{\varpi }+\frac{G}{\varpi }{\vec{e}}_{\hat{\varphi }}.$$

$$\alpha \frac{\partial }{\partial r}\left(\alpha \frac{\partial F}{\partial r}\right)+\frac{\sin \theta }{r^2}\frac{\partial }{\partial \theta }\left(\frac{1}{\sin \theta }\frac{\partial F}{\partial \theta }\right)=-\frac{1}{\alpha }S.$$

$${\alpha }^2\frac{\partial }{\partial r}\left({\alpha }^2\frac{\partial G}{\partial r}\right)+\frac{{\alpha }^2\sin \theta }{r^2}\frac{\partial }{\partial \theta }\left(\frac{1}{\sin \theta }\frac{\partial G}{\partial \theta }\right)=-S{S}^{\prime }.$$

$$E_{\mathrm{EM}}=\int \frac{\alpha B^2}{8\pi }\sqrt{g_3}{d}^{3}x=\frac{1}{4}\int \left[{\left(\alpha \frac{\partial G}{\partial r}\right)}^2+{\left(\frac{1}{r}\frac{\partial G}{\partial \theta }\right)}^2+{\left(\frac{S}{\alpha }\right)}^2\right]\frac{drd\theta }{\sin \theta },$$

$$H_{\mathrm{R}}=\int \vec{A}\cdot \vec{B}\sqrt{g_3}{d}^{3}x=4\pi \int \frac{GS}{{\alpha }^2}\frac{drd\theta }{\sin \theta },$$

$$\frac{{\alpha }^{-1}dr}{B_{\hat{r}}}=\frac{rd\theta }{B_{\hat{\theta }}}=\frac{\varpi d\varphi }{B_{\hat{\varphi }}}.$$

$$\Delta \varphi ={\int }_{r_1}^{r_2}\frac{rS}{{\alpha }^2\varpi G_{,\theta }}dr=-{\int }_{{\theta }_1}^{{\theta }_2}\frac{rS}{{\alpha }^2\varpi G_{,r}}d\theta .$$

$$G_{\mathrm{d}}$$

$$E_0=\frac{B_0^2{R}^3}{3}\left[1+\frac{5M}{2R}+\frac{51M^2}{10R^2}+\cdots \right].$$

$$S={\left(\frac{2\gamma }{n+1}\right)}^{1/2}{G}^{(n+1)/2}.$$

$$G=G_{\mathrm{d}}(R,\theta ).$$

$$G(r,\theta )=-\sum _{l\ge 1}{g}_l(r)\sin \theta \frac{d{P}_l(\theta )}{d\theta },$$

$$\frac{d}{dr}\left({\alpha }^2\frac{d}{dr}{g}_l\right)-\frac{l(l+1)}{r^2}{g}_l=0.$$

$$G(r,\theta )=-\sum _{l\ge 1}(a_l{g}_l^{-}(r)+b_l{g}_l^{+}(r)+q_l(r))\sin \theta \frac{d{P}_l(\theta )}{d\theta }.$$

$$q_l(r)=g_l^{-}(r){\int }_R^r\frac{g_l^{+}(r^{\prime })S_l(r^{\prime })}{W_l(r^{\prime })}d{r}^{\prime }+g_l^{+}(r){\int }_r^{r_{\infty }}\frac{g_l^{-}(r^{\prime })S_l(r^{\prime })}{W_l(r^{\prime })}d{r}^{\prime },$$

$$W_l\equiv {\alpha }^4\left[g_l^{-}\frac{d{g}_l^{+}}{dr}-g_l^{+}\frac{d{g}_l^{-}}{dr}\right]$$

$$S_l(r)=\frac{2l+1}{2l(l+1)}\left[{\int }_0^{\pi }\mathcal{S}(r,\theta )\frac{d{P}_l(\theta )}{d\theta }d\theta \right].$$

Full text

††pagerange: Axisymmetric force-free magnetosphere in the exterior of a neutron star–LABEL:lastpage

1 Introduction

The magnetic fields of pulsars are inferred from the observed spin-down rate. The intensity at the surface exceeds the QED critical field $4.4\times 10^{13}$G in a class of neutron stars conventionally called magnetars (Thompson & Duncan, 1995, 1996), which typically exhibit energetic behaviors. The observed X-ray luminosity is much higher than the spin-down luminosity in both anomalous X-ray pulsars (AXPs) and soft $\gamma$-ray repeaters (SGRs). They are very variable, and sometimes exhibit bursts, unlike normal pulsars. In giant flares of SGRs, huge amounts of energy $10^{44}$-$10^{46}$ erg are released in an instant. The activity of magnetars is widely believed to be supplied by their strong magnetic field of $10^{14}$-$10^{15}$G. Recent observations have hindered a definite neutron star classification. There are exceptional bursting sources with normal magnetic field ($<4\times 10^{13}$G), i.e, ‘low-field’ magnetar (Rea et al., 2010). Outbursts are observed in a rotation-powered pulsar(Göğüş et al., 2016). See, for example, Mereghetti (2008); Turolla et al. (2015) for a review of magnetars.

Two possibilities are proposed for the energy storage prior to magnetar outbursts to explain the relevant phenomena: storage in the magnetar crust or in the magnetosphere (Thompson et al., 2002). The latter model is discussed in terms of similarity with solar flares (Lyutikov, 2003, 2006; Beloborodov & Thompson, 2007). In the solar flare model (e.g., Aly, 1984; Mikic & Linker, 1994), the energy is quasi-statically stored by thermal motion at the surface, and is suddenly released as large-scale eruptive coronal mass ejections. The energy is dissipated via a magnetic reconnection associated with the field reconfiguration. Analogous energy buildup and release processes may be relevant to the magnetar giant flares, although the energy scale differs by many orders.

The twisted magnetosphere formed by current flow is studied in the calculations of synthetic spectra (Lyutikov & Gavriil, 2006; Fernández & Thompson, 2007; Nobili et al., 2008) and polarization (Taverna et al., 2014). In their calculations, the external global magnetosphere is described by a self-similar bipolar field. The solution of the force-free field is given by Lynden-Bell & Boily (1994); Thompson et al. (2002), and is extended to a self-similar one with higher order multipoles (Pavan et al., 2009). Recently, Akgün et al. (2016) numerically constructed the force-free magnetosphere of a magnetar, which is smoothly matched to a vacuum solution at large distance.

The magnetar magnetospheres are also studied to be matched with the interior magnetic field (e.g., Viganò et al., 2011; Glampedakis et al., 2014; Fujisawa & Kisaka, 2014; Pili et al., 2015). The exterior magnetic fields depend on the interior ones through the boundary condition. The magnetic field structure evolves on a secular timescale by Hall drift, ambipolar diffusion, or some other mechanism (e.g., Goldreich & Reisenegger, 1992; Hollerbach & Rüdiger, 2004; Kojima & Kisaka, 2012; Viganò et al., 2013; Gourgouliatos & Cumming, 2014; Wood & Hollerbach, 2015). A quasi-steady shearing motion at the footpoint of a magnetic field twists the exterior field, and the stored energy gradually increases. When the state exceeds a threshold, the energy is abruptly released on a much shorter dynamical timescale, leading to energetic flares. One important fundamental process is the accumulation of magnetic helicity in the corona before the outbursts. The magnetic energy also builds up as a natural product of the evolution of the magnetic helicity. The role of helicity accumulation in the solar corona has been discussed (Zhang et al., 2006). The same process should be considered in the framework of general relativity, since the buildup of magnetic energy occurs very near the neutron star surface.

Pili et al. (2015) obtained axially symmetric equilibrium by solving relativistic MHD equations both in the interior and in the exterior of a neutron star. The magnetosphere is modeled by filling low-density plasma, and the magnetic field is twisted by the current. Magnetic pressure is much larger than hydrodynamical pressure, and so the force-free approximation may be applicable. We here provide a systematic study of a twisted force-free magnetosphere by limiting our consideration to the exteriors. The force-free magnetosphere around a black hole has been examined (e.g., Uzdensky, 2004). Our problem is simpler because the charge density associated with plasma rotation is ignored. The magnetars are slow rotators, and the current flow is more important.

In this paper, we examine the effect of general relativity on the force-free magnetosphere in the vicinity of neutron stars. We consider static and axially symmetric magnetic fields in the exterior of a non-rotating compact star. The spacetime is described by the Schwarzschild metric, and a family of equilibrium solutions are calculated assuming that the current function is given by a simple power-law model, for which solutions in a flat spacetime have been well studied by Flyer et al. (2004); Zhang et al. (2006, 2012) in the context of solar flares. This will clarify the effect of general relativity. Equations relevant to non-rotating force-free magnetosphere in a Schwarzschild spacetime are discussed briefly in Section 2. We numerically solve the so-called Grad–Shafranov equation111 Although this name is used in many papers, a work by L$\ddot{\rm u}$st and Sch$\ddot{\rm u}$ter prior to Grad and Shafranov, was not recognized, as noted by Akgün et al. (2016). with a power-law model. In Section 3, the results are given for comparison with those obtained in flat spacetime. Finally, our conclusions are given in Section 4. We use geometric units of $c=G_{\rm N}=1$.

2 Equations

2.1 Magnetic fields

In this section, we briefly summarize the Maxwell equations in curved spacetime outside of a non-rotating compact object. The metric in the exterior of a mass $M$ in spherical coordinates is given by

[TABLE]

where

[TABLE]

We consider the static magnetic configuration. Using vector analysis in 3-dimensional curved space, the relevant equations are expressed as (e.g., Thorne et al., 1986)

[TABLE]

[TABLE]

The magnetic field for the axially symmetric case is given by two scalar functions $G(r,\theta)$ and $S(r,\theta)$:

[TABLE]

where $\vec{e}_{\hat{\phi}}\equiv\varpi^{-1}\partial_{\phi}$ is a unit vector in the azimuthal direction. The magnetic field (5) automatically satisfies the constraint equation (3), and the components are explicitly written as

[TABLE]

We introduce a vector potential ${\vec{A}}$ derived by functions $F(r,\theta)$ and $G(r,\theta)$ to express the magnetic field as ${\vec{B}}={\vec{\nabla}}\times{\vec{A}}$:

[TABLE]

There is gauge freedom in the potential ${\vec{A}}$. The electromagnetic fields are unchanged by ${\vec{A}}\to{\vec{A}}^{\prime}={\vec{A}}+{\vec{\nabla}}\Lambda$. In the axially symmetric case ($\partial_{\phi}=0$), $G(=A_{\phi})$ is gauge-invariant, and has a special meaning of magnetic flux. The poloidal part of the vector potential is determined by

[TABLE]

The force-free magnetic field ${\vec{j}}\times{\vec{B}}=0$ means that current flows along the magnetic field lines. The poloidal component of eq.(4) (the Biot-Savart law) is reduced to $\nabla S||\nabla G$, that is, $S=S(G)$, and $4\pi\alpha{\vec{j}}_{p}=S^{\prime}{\vec{B}}_{p}$, where $S^{\prime}=dS/dG$. This means that the poloidal current flows along a line of constant $S$. A factor $\alpha$ is needed in front of the current, since not ${\vec{j}}_{p}$ but $\alpha{\vec{j}}_{p}$ should be proportional to ${\vec{B}}_{p}$ from the charge conservation $\partial\rho_{e}/\partial t+{\vec{\nabla}}\cdot(\alpha{\vec{j}})=0$. Using this form, the azimuthal component is given by $4\pi\alpha j_{\phi}=S^{\prime}B_{\phi}$ $=SS^{\prime}/(\alpha\varpi)$. The azimuthal component of eq. (4) therefore reduces to the so-called Grad-Shafranov equation:

[TABLE]

This equation can be also derived in the non-rotating limit of the equation used for black hole magnetospheres (Uzdensky, 2004).

Magnetic energy is given by integrating over 3-dimensional volume

[TABLE]

where $\sqrt{g_{3}}(=\alpha^{-1}r^{2}\sin\theta)$ is a determinant of 3-dimensional space metrics. We consider the energy deposited outside the star, so that the integration range is limited to $r\geq R$. In eq. (10), a factor $\alpha$ in front of $B^{2}$ may be understood by considering the Maxwell equations in curved spacetime. Equivalently, the expression (10) can also be obtained by $\int T^{t}_{t}\sqrt{-g_{4}}d^{3}x$ in terms of energy momentum tensor $T^{t}_{t}$ and determinant of 4-dimensional spacetime metrics $\sqrt{-g_{4}}$. Helicity is obtained by integrating the product of two vectors, namely, ${\vec{A}}$ and ${\vec{B}}$. Total relative helicity $H_{\rm R}$ defined by the integration range $r\geq R$ is

[TABLE]

where the poloidal part derived by $F$ in the vector potential is eliminated by integrating by parts with eq.(8). In the limit of flat spacetime ($\alpha\to 1$), eq. (11) reduces to the equation given by Zhang et al. (2006). When $S$ is a linear function of $G$, eq. (9) is derived by minimizing the magnetic energy (10) for fixed helicity (11). Thus force-free equilibrium with constant $S^{\prime}$ is a minimum energy state. This is a relativistic extension of the work by Chandrasekhar & Woltjer (1958).

In presence of the toroidal component, the magnetic field lines deviate in the azimuthal direction. We introduce a twist angle, which is closely related to helicity. The magnetic field lines satisfy the relationship:

[TABLE]

The twist angle along the line from $(r_{1},\theta_{1})$ to $(r_{2},\theta_{2})$ is evaluated by one of the following expressions:

[TABLE]

2.2 Model and numerical methods

The potential dipole field, that is, the solution of eq. (9) for $S=0$ with magnetic moment $\mu$, is given elsewhere (e.g., Petterson, 1974; Wasserman & Shapiro, 1983; Konno et al., 1999):

[TABLE]

where the first expression is the exact solution and the second is an approximation in the weak gravity regime $M/r\ll 1$. Using magnetic dipole moment $\mu$, the surface field strength is given by $B_{0}\equiv\mu R^{-3}$, where $R$ is the stellar radius. The magnetic field at the stellar surface is $B_{\hat{r}}\approx(1+3M/(2R))\times 2B_{0}\cos\theta$ and $B_{\hat{\theta}}\approx(1+2M/R)\times B_{0}\sin\theta$ from eq. (6). The field strength at the pole is $2B_{0}$ in flat spacetime, but increases by some factor in a relativistic treatment for a fixed $\mu$. The correction depends on the direction, but is typically a factor $\approx 1.5$ for $M/R\approx 0.25$. The magnetic energy $E_{0}$ stored in the exterior $(r\geq R)$ for the solution (2.2) is expressed by

[TABLE]

The energy increases by a factor $\approx 2$ for $M/R\approx 0.25$ and a fixed $\mu$. In this paper, we use $B_{0}$ and $R$ for normalization, because $B_{0}$ is usually estimated from measurement of spin period and its time derivative of pulsars.

Our concern is how energy and helicity change by twisting the magnetic field in the exterior of a relativistic star. We adopt a power-law model for $S$ as

[TABLE]

The right-hand side of eq. (9) thus becomes $-{\mathcal{S}}\equiv-SS^{\prime}=-\gamma G^{n}$, where $\gamma$ is constant, and $n$ is assumed to be an odd integer in order to ensure that $S^{2}$ is positive definite independent of the sign of $G$. The force-free magnetosphere with the power-law model has been extensively studied in flat spacetime for the solar flare model (Flyer et al., 2004; Zhang et al., 2006, 2012). There is no solution different from that of a pure dipole in the case of index $n<5$ in the flat spacetime treatment. We calculate the models with $n=5$ and 7 in this paper. A continuous sequence of solutions is thus generated by varying a free parameter. However, the sequence may not necessarily correspond to evolutionary track. Physical evolution of the magnetosphere is not obtained until we know the evolution of footpoint motion. As a possible scenario, suppose that the magnetosphere gradually changes by current flow from the surface. The magnetic field is twisted, and the field topology changes by helicity accumulation. Not only the parameter $\gamma$ but also the functional form $\gamma G^{n}$ should be changed in the power-law description in order to describe the evolution. In spite of this fact, the sequence of solutions with a fixed index $n$ is easily constructed with an increase of magnetic twist, and provides an understanding of the structure change by helicity accumulation.

The boundary conditions for solving eq. (9) are as follows. At the stellar surface $r=R$, we assume that the magnetic function $G$ is given by that of a dipole in vacuum, that is, eq. (2.2):

[TABLE]

At asymptotic infinity ($r\to\infty$), the function is expected to decrease as $G\propto r^{-1}$. At the polar axis, the function is expected to satisfy the regularity condition, that is, $G=0$ at $\theta=0$ and $\pi$.

Equation (9) with the power-law model (16) is non-linear except for $n=1$, and an iterative method is therefore necessary for solving it. An initial guess for the function $G_{(1)}$ is assumed. Using this, the right hand side of eq. (9) is calculated, and a new function $G_{(2)}$ is solved from the source term with appropriate boundary conditions. The procedure is repeated until convergence, say, $||G_{(k+1)}-G_{(k)}||$ $<\varepsilon$, where $\varepsilon$ is a small number, and the left hand side is evaluated by discrete points $(r_{i},\theta_{j})$ inside the domain: $||G_{(k+1)}-G_{(k)}||\equiv$ $[\sum_{i,j}(G_{(k+1)}(r_{i},\theta_{j})-G_{(k)}(r_{i},\theta_{j}))^{2}]^{1/2}$ $/[\sum_{i,j}G_{(k)}(r_{i},\theta_{j})^{2}]^{1/2}$. The iteration scheme may not necessarily lead to a convergent solution, but we experimentally have a convergent solution for certain initial trials.

One of the numerical methods for solving eq. (9) is discretizing the equation (e.g., Akgün et al., 2016). Here we use another approach by polynomial expansion (e.g., Fujisawa & Kisaka, 2014). The magnetic function $G$ is expanded by Legendre polynomials $P_{l}(\theta)$:

[TABLE]

where radial functions $g_{l}$ in the case of the homogeneous equation of eq. (9) satisfy:

[TABLE]

Two independent solutions are obtained as $r^{l+1}$ and $r^{-l}$ in flat spacetime. The solutions can be expressed by hypergeometric function in a Schwarzschild spacetime. They are constructed by numerically solving eq. (19). We denote two independent solutions by $g^{-}_{l}$ and $g^{+}_{l}$. They are characterized by asymptotic behavior at large $r$: $g^{-}_{l}\to r^{-l}$ and $g^{+}_{l}\to r^{l+1}$. The general solution of eq. (9) can be written using constants $a_{l},b_{l}$ and a particular solution $q_{l}(r)$:

[TABLE]

When the source term $-{\mathcal{S}}=-SS^{\prime}$ in eq. (9) is assumed to be given by an iterative process, the inhomogeneous solution $q_{l}$ can be solved by the Green function method:

[TABLE]

where $W_{l}$ is a radial function defined by

[TABLE]

and $S_{l}$ is the $l$-th component of the source term ${\mathcal{S}}$:

[TABLE]

Although the upper limit $r_{\infty}$ of the integral should be $\infty$, $r_{\infty}/R>$ 20-30 is sufficient in the actual numerical calculation. The inhomogeneous solution in eq. (21) satisfies $q_{l}\propto r^{-l}$ in the limit of $r\to\infty$, so that all coefficients $b_{l}$ in front of growing functions should be zero. The condition at the surface, eq. (17) determines the coefficients $a_{l}$. In the numerical calculation, we truncate the polynomial expansion in eq. (18) up to $l_{\rm max}=$64-128. The numerical results are almost unchanged by further increasing it.

3 Numerical results

We first construct a magnetosphere model with a dipole in vacuum, and gradually increase the toroidal component for a given power-law index $n$ and a relativistic factor $M/R$. Calculations in flat spacetime (Flyer et al., 2004) indicate that a free parameter $\gamma$ is not suitable for specifying the sequence of models, because there are multiple solutions for a given $\gamma$. Instead, the azimuthal magnetic flux or helicity is used as the degree of twist, and the constant $\gamma$ is posteriorly determined. We follow the same method, and calculate the energy and the relative helicity along the sequence of solutions. The characteristics of the models are given in Table 1.

The increase in magnetic energy and the relative helicity are shown for four models in Fig. 1. The results in flat spacetime agree with those by Flyer et al. (2004); Zhang et al. (2006). For a better understanding of the mechanism, the energy difference $\Delta E(=E_{\rm EM}-E_{0})$ is divided into $\Delta E=\Delta E_{t}+\Delta E_{p}$, that of toroidal component and that of the poloidal component. With increasing helicity from zero in the weak regime, the energy $\Delta E_{t}$ increases monotonically, while the energy $\Delta E_{p}$ is almost zero. This phase, which corresponds to the lower branch of the sequence in Fig. 1, is characterized by increasing strength parameter $\gamma$. There is a maximum of $\gamma$, and after passing the turning point, the energy $\Delta E_{p}$ drastically increases in the upper branch in Fig. 1. This means a significant structure change from that in the dipole potential field. This fact will be confirmed later. The numerical calculation becomes complicated for finding equilibrium solutions with larger helicity. This terminates our numerical calculations. The endpoint of the sequence is mostly consistent with that of Flyer et al. (2004); Zhang et al. (2006). It is likely that a static configuration no longer exists, leading to a catastrophic transition associated with ejection of magnetic energy upon further accumulation of helicity.

In the relativistic treatment, the increase of magnetic energy is much larger. This does not merely come from the normalization of the surface field. As discussed in the previous section, the surface field in the relativistic model is a few times larger than in the flat spacetime model with the same magnetic dipole moment $\mu$. Therefore, $\Delta E$ itself is naturally scaled up. The ratio of $\Delta E$ to $E_{0}$, which does not depend on the normaliztion, is also significantly increased. The maximum of each model is given in Table 1. The curved spacetime allows higher energy and helicity to be stored in the magnetosphere. Nonlinearity with respect to the index $n$ also enhances the maximum energy. For example, the maximum in the model with $n=7$ is larger than for $n=5$. The relativistic correction, however, is much larger.

Figure 2 shows the magnetic function $G$ obtained in flat spacetime by contours in the $r$-$\theta$ plane for the inner region $(r/R\leq 5)$. The magnetic field structure with large helicity is significantly changed, as shown in the second and third panels. By increasing the helicity, the magnetic field lines become stretched outwardly and slightly squashed into the equator. Flyer et al. (2004); Zhang et al. (2006) discovered an interesting property of the highly twisted state: A local maximum of $G$ exists just outside the star at the equator. The peak, which becomes more evident for larger $n$, is located inside the stretched poloidal field. This is numerically confirmed in our model, but the peak is too small to appear in the coarse interval of the contour levels.

Results for the relativistic treatment are shown in Fig. 3. As discussed with respect to eq. (2.2), the magnetic field strength near the surface is increased by a few factors for the same magnetic dipole moment $\mu$. The structure of the potential magnetic dipole shown in the left panel is the same as obtained in flat spacetime except for the magnitude. However, there are evident loops of the field lines in models with large helicity, as shown the second and third panels in Fig. 3. This structure represents a magnetic flux rope with field line projections on the $r$-$\theta$ plane in the form of a closed curve of a constant $G$. The center of the loops ($r/R\approx 1.5$, $\theta=\pi/2$) corresponds to a global maximum of $G$, whereas the center of potential dipole field $G_{\rm d}$ is located at the stellar surface, $r=R$ and $\theta=\pi/2$. The loop structure is clearly evident in the model with $n=7$ due to the stronger nonlinearity of $n$. It should be noted that the loop structure is formed in the upper branch of the sequence in Fig. 1. It has higher energy and helicity for a fixed $\gamma$. It is therefore very important to know whether this kind of higher energy state appears evolutionarily. A helically twisted flux rope is thought to be a phase that precedes magnetar flares (Yu, 2011; Huang & Yu, 2014). The structure is modeled by a small ring current at $(r,\theta)=(r_{0},\pi/2)$. The results show that for a stable equilibrium solution the radial distance $r_{0}$ is limited to a few times the stellar radius. Our present result is an explicit example of a fixed rope.

Figure 4 shows a 3-dimensional structure for the highly twisted relativistic model with $n=5$. Twisted lines anchored to the central star are shown in only the hemisphere ($0\leq\theta\leq\pi/2$). The figure also shows that the typical twist angle $\Delta\phi$ is 1 radian, where $\Delta\phi$ denotes the angle from the stellar surface to the equator in eq. (13). It can be easily understood that the maximum of $\Delta\phi$ for a magnetosphere model is given by a line starting from the surface with an intermediate latitude $\theta_{*}$ ($0<\theta_{*}<\pi/2$). However, the exact value $\theta_{*}$ depends on the structure, and numerical work is necessary to find it by calculating $\Delta\phi$ for all possible lines. We search the maximum among calculated sequences of four models. Their values are tabulated in Table 1, and are all approximately $\Delta\phi_{\rm max}\approx 1$ radian. It should be noted that the maximum in the model sequence is not necessarily attained in the state with the largest helicity. The maximum in the G7 model is smaller, although a large amount of helicity can be stored. This fact is related to the loop formation. When this kind of detached structure forms, the lines starting from the star and returning to it do not extend radially. The twist in the azimuthal direction becomes small. The maximum twist $\Delta\phi_{\rm max}=$0.8-1.5 in our models is interesting by comparison with a different model: $\Delta\phi_{\rm max}=1.2$-1.5 in the model of Akgün et al. (2016), and $\Delta\phi_{\rm max}\leq\pi/2$ in a self-similar model (Thompson et al., 2002) 222 Our twisted angle $\Delta\phi$ is defined between a point in north-hemisphere to equator. The twisted angle $\Delta\phi_{\rm NS}$ in (Thompson et al., 2002) is $\Delta\phi_{\rm NS}=2\Delta\phi$. . Finally, it should be noted that, in a sense, the quantity $\Delta\phi$ represents the degree of twisted structure, but is not necessarily a good indicator because it is not a monotonic function along a sequence.

4 Conclusion

In this paper, we have studied how a force-free magnetosphere is changed by the accumulation of helicity. The model of current is specified by a simple power law by which the equilibria in flat spacetime are extensively examined in the context of solar flares. We extend the study to the general relativistic regime for application to magnetar flares. A family of axially symmetric static solutions is calculated. Total magnetic energy increases with increasing helicity. Numerical calculations suggest that there is a maximum of the helicity deposited in the magnetosphere. Beyond this, we could not find any static solutions. The existence of a critical state, which is characterized by magnetic helicity or twist angle in Thompson et al. (2002); Zhang et al. (2006); Akgün et al. (2016), is also found in the different models and by different numerical approaches. This fact is very interesting for the eruption of magnetic flux rope. When the total helicity injected into the magnetosphere exceeds a certain capacity, a catastrophic event such as a flare may occur. This kind of sudden event may be calculated by resistive dynamical simulation (Li et al., 2012; Parfrey et al., 2013; Kojima & Kato, 2014). It is also important to determine when an instability sets in along our equilibrium sequence: at the endpoint of our calculation or at a certain state before it.

We found that general relativistic effects are significant. A larger amount of energy is capable of being stored in the relativistic magnetosphere. In extreme cases, the increase in magnetic energy in the presence of current flow exceeds that of a current-free dipole. This marks a contrast with the models in flat spacetime, in which the magnetic energy increment is less than that of vacuum dipole. When the magnetic field structure is highly twisted, a loop of field lines in 2-dimensional meridian plane, which corresponds to an axially symmetric torus in 3-dimensional space, is formed in the vicinity of the surface. This kind of remarkable structure is also realized in flat spacetime, but the condition to the power-law index $n$ is severe. Namely, evident structure is realized only in the model with stronger nonlinearity, $n\geq 7$ (Flyer et al., 2004; Zhang et al., 2006). The loops are also found in the exterior model of a neutron star obtained by solving the relativistic MHD equation with a low-density plasma (Pili et al., 2015). A direct comparison is difficult due to the differences in both models and numerical methods, but general relativity is likely to play an important role in the formation.

Observationally, total magnetic energy in the magnetar magnetosphere is approximately $\sim 10^{48}(B_{0}/10^{15}{\rm G})^{2}(R/10{\rm km})^{3}$. The flare energy is about 1% of this energy. This means that the flare is associated with a small change between the equilibrium models. From the point of the energy scale, there are two possibilities: a small transition between twisted dipolar equilibria or a transition of a higher multipolar field with weaker strength. Long-term evolution of the internal magnetic field by Hall drift demonstrates that this kind of higher multipolar field emerges at the surface (e.g., Hollerbach & Rüdiger, 2004; Kojima & Kisaka, 2012; Viganò et al., 2013; Gourgouliatos & Cumming, 2014; Wood & Hollerbach, 2015), although the results depend on the initial conditions. A correct functional form $S=S(G)$ at the surface, which is important for constructing force-free magnetosphere, will be determined by connecting internal evolution.

Acknowledgements

This work was supported in part by a Grant-in-Aid for Scientific Research (No. 25103414) from the Japanese Ministry of Education, Culture, Sports, Science and Technology.