{'fraction_non_alphanumeric': 0.030067983756901034, 'fraction_numerical': 0.01834192635853447, 'mean_word_length': 4.872990353697749, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hydromagnetic Instabilities in Neutron Stars

We model the non-linear ideal magnetohydrodynamics of poloidal magnetic fields in neutron stars in general relativity assuming a polytropic equation of state. We identify familiar hydromagnetic modes, in particular the 'sausage/varicose' mode and 'kink' instability inherent to poloidal magnetic fields. The evolution is dominated by the kink instability, which causes a cataclysmic reconfiguration of the magnetic field. The system subsequently evolves to new, non-axisymmetric, quasi-equilibrium end-states. The existence of this branch of stable quasi-equilibria may have consequences for magnetar physics, including flare generation mechanisms and interpretations of quasi-periodic oscillations.


INTRODUCTION
Neutron stars harbour the strongest known magnetic fields in Nature. It is therefore not surprising that the magnetic field plays a key role in many aspects of neutron star physics. Most of what is known about neutron star magnetism comes from the inferred dipole moment associated with the electromagnetic stellar spin-down. In systems with prominent thermal emission, the surface magnetic field strength can be independently measured by its influence on the emission pattern. Remarkably, various astronomical observations suggest a systematic difference in field strengths between different classes of neutron stars (see Kaspi 2010, for a recent review). Aged systems, such as millisecond pulsars and accreting neutron stars in low-mass X-ray binaries, have relatively weak magnetic fields, B ∼ 10 8 G. Typical radio pulsars have B ∼ 10 12 G. The highest field strengths are associated with soft-gamma-repeaters (SGRs) and anomalous X-ray pulsars (AXPs) with magnetic fields B ∼ 10 15 G.
This last class of neutron stars, collectively known as magnetars, are the focus of this paper. These slowly-spinning objects exhibit high energy emission, which is occasionaly punctuated by bursts. In some rare occasions, powerful flares have been observed in SGRs. According to the prevailing magnetar model, as first put forward by Duncan & Thompson (1992); Thompson & Duncan (1995), the energy reservoir powering magnetar activity is the energy of the magnetic field itself. The unusually high magnetar surface temperatures are likely the result of heating produced by the Ohmic decay of the magnetic field in the neutron star crust (e.g. Pons et al. 2007). More spectacularly, giant SGR flares could be the manifestation of global magnetic field instabilities and of subsequent field reconfiguration (Thompson & Duncan 1995). Furthermore, there is strong evidence that quasi-periodic oscillations (QPOs) observed in the light curves of these giant flares are of magnetoelastic nature (i.e. the restoring force is provided by magnetic and crustal elastic stresses) (e.g. Levin lasky@tat.physik.uni-tuebingen.de 2006Sotani, Kokkotas & Stergioulas 2008;Cerdá-Durán, Stergioulas & Font 2009;Gabler et al. 2011). Unfortunately, a detailed understanding of magnetar magnetohydrodynamics is hindered by the lack of information about the properties of the interior magnetic field. Such an understanding is a vital element if we are to carry out, for example, magnetar 'asteroseismology' using the observed QPOs.
The aforementioned reasons imply it is of great theoretical and astrophysical interest to model the dynamics of magnetic fields in the interior of magnetars (and in general of neutron stars). More specifically, three key questions have attracted the most attention to date: (i) What is the nature of long-term magnetic field equilibria in neutron stars? (ii) Do hydromagnetic instabilies play a role in phenomena like the bursts and flares seen in AXPs and SGRs? (iii) Which physical mechanism(s) determines the long-term evolution and decay of magnetic fields in neutron stars?
In this letter we make contact with these first two questions. We study short-term (i.e. dynamical timescales) neutron star magnetohydrodynamics (MHD), specifically addressing the issue of global MHD instabilities in neutron stars and resulting magnetic field equilibria. The last few years have seen a flare of activity on this subject, mainly motivated by magnetar astrophysics. A large body of work has been devoted to semi-analytical studies of axisymmetric MHD equilibria (e.g. Ioka 2001;Haskell et al. 2008;Ciolfi et al. 2009;Ciolfi, Ferrari & Gualtieri 2010), building on much earlier work regarding equilibria (e.g. Monaghan 1965;Roxburgh 1966;Parker 1966) and stability analyses (Kruskal & Schwarzschild 1954;Tayler 1957Tayler , 1973Wright 1973;Markey & Tayler 1973, 1974Flowers & Ruderman 1977). It is only recently that global, fully non-linear, MHD numerical calculations have appeared studying magnetic field stability and equilibria in stellar models (Geppert & Rheinhardt 2006;Braithwaite & Spruit 2006;Braithwaite & Nordlund 2006;Braithwaite 2007Braithwaite , 2008Braithwaite , 2009Duez, Braithwaite & Mathis 2010;Kiuchi, Yoshida & Shibata 2011). The work presented herein falls into this last category. We study the dynamical evolution of purely poloidal magnetic fields in neutron stars in general relativity, investigating intrinsic instabilities and subsequent equilibria, uncovering for the first time non-axisymmetric magnetic field equilibria in barotropic stars.

NUMERICAL MODEL
We study magnetized neutron stars through time evolution of the ideal MHD equations in general relativity under the Cowling approximation. To this end, we utilise the three-dimensional GRMHD code thor (Zink, Schnetter & Tiglio 2008;Korobkin et al. 2010) and her sister GPU code horizon (Zink 2011).
The MHD portion of the code utilizes the conservative formalism outlined in Anton et al. (2006). Hyperbolic divergence cleaning is employed according to Anderson et al. (2006). We treat the exterior of the star as a low-density artificial atmosphere, allowing full evolution of the magnetic field in this region. The outer boundary of our domain, located at approximately 1.4 times the stellar radius, uses Dirichlet boundary conditions for the magnetic field. We have tested our evolutions against this condition and found no discernible difference with the nature of the instability or the endstate of the magnetic field configuration.
We employ a Cartesian grid with 120 3 grid points at single precision. We have completed numerous resolution studies including 90 3 , 150 3 and 200 3 grid points, finding consistent phenomenology with converging growth-times for the instability. Moreover, we have tested the use of double precision on these simulations and again found consistent phenomenology and timescales.

Initial Conditions
We employ the spectral code lorene 1 , which produces self-consistent solutions of the Einstein-Maxwell field equations in ideal MHD (Bocquet et al. 1995). The spectral grid is mapped to our Cartesian grid and the system is allowed to evolve. No additional perturbations are imposed, relying instead on intrinsic numerical noise.

Mode Decomposition
The modal structure of the instability is extracted by performing a Fourier decomposition of various physical quantities on a ring in the equatorial plane. We compute complex weighted averages (Zink et al. 2007) where ̟ = x 2 + y 2 = const. lies in the initial equatorial plane of the magnetic field. In equation (1), f is some quantity in the star chosen to best represent the instability. Throughout much of this article we use the φ and z components of the magnetic field, although note the velocity and density fields also exhibit the instability, however the signal is generally not as clean. Our fiducial model is a non-rotating star with polytropic equation of state relating pressure, P , and restmass density, ρ, through P = Kρ Γ , with Γ = 2 and K = 100. This gives a stellar model with gravitational mass M = 1.3 M ⊙ and equatorial radius R = 12.6 km. Our fiducial model has an average field strength inside the star ofB = 1.3 × 10 16 G, yielding an Alfvén crossing time of t A = 5.0 ms. Such models have polar magnetic field strengths of B surf = 8.8×10 15 G, which is a factor of a few larger than observations of the dipole component of the strongest magnetic fields. This is predominantly done to reduce computational costs.

RESULTS
In figure 1 we plot the evolution of m = 1, . . . , 4 modes for C m (B φ ) for our fiducial model, measured at ̟ = 0.6̟ ⋆ , where ̟ ⋆ is the stellar radius. An instability is present in all modes after the first couple of Alfvén crossing times. Each mode grows exponentially by approximately six orders of magnitude over many subsequent Alfvén crossing times. Saturation of the modes occurs after about 75 ms, at which point the simulation evolves to a pseudo-equillibrium state. We have evolved such simulations to 400 ms = 80 t A , with little variation in the equilibrium configuration following the first hundred or so milliseconds. We discuss the equilibrium state in more detail below.
In the first few milliseconds of the simulation we see the m = 4 mode in figure 1 grow by an order of magnitude over the other modes. This feature is somewhat aligned with our Cartesian grid, which we attribute to a finite differencing effect associated with the surface of the star. To understand this further we performed resolution studies with 150 3 and 200 3 grid points. In such simulations we see the initial transient that lifts the m = 4 mode above other modes slowly converging away with increasing resolution.
In figures 2 we present three-dimensional plots of the magnetic field at various instances throughout the evolution 2 . The instability acts near the neutral line, which is the line where B = 0 around the equatorial plane, located in our simulations at approximately two-thirds of the stellar radius. For clarity, we have plotted red field lines seeded near the neutral line in the equatorial plane. Also plotted are black field lines seeded on the equatorial plane interior to the neutral line. The blue volume rendering is an isopycnic surface of ρ = 0.37ρ c , where ρ c is the central rest-mass density, which lies at a radius of approximately 50% of the star; well inside the neutral line. Figure 2a shows the initial data imported from the lorene spectral code. The domain of our grid is larger than that plotted here and field lines are truncated at the surface of the star for clarity. Such a configuration is known to be unstable to the kink instability from local linear studies in Newtonian physics (Markey & Tayler 1973). Such studies have been confirmed with global, linearized, numerical evolutions (Lander & Jones 2011) and in full non-linear studies (Braithwaite 2007). This is the first time such a configuration has been studied in general relativity. Figure 2b shows the evolution after 25 ms = 5 t A . The onset of the 'sausage' or 'varicose' mode (Markey & Tayler 1973), involving changes in the crosssectional area of a flux tube around the neutral line, is clearly visible. This is strongest in the m = 4 mode, which is a result of the transient excitation at the beginning of our simulation. Whilst this transient reduces with increasing grid resolution, the presence of the varicose mode is an inherent characteristic of the system (as evidenced from its almost constant quasi-equilibra state that is resolution independent).
The varicose mode is visually present in our simulation for approximately eight or nine Alfvén crossing times before the 'kink' instability appears and begins to dominate. In figure 2c we show the onset of this mode acting perpendicular to the gravitational field, in accordance with the prediction of Markey & Tayler (1973). This figure is after 50 ms = 10 t A , at which point we can still see the presence of the varicose mode; i.e. the red field lines still exhibit varying cross-sectional area around the star. It is worth noting that this mode will have been excited from the beginning of the simulation (as seen in figure 1), however it is only now that the exponential growth has reached a point at which it is visually obvious. In many ways this represents the non-linear development of the instability where the change in field structure is of similar order to the background field. It is also worth mentioning here that our choice of the Cowling approximation is well justified considering the dominant instability acts on equipotential surfaces.
The modal analysis presented in figure 1 does not distinguish between the varicose or kink modes, implying exponential growth of these quantities does not necessarily indicate an instability in one or the other mode. Indeed, figures 2b and 2c suggest that, if the varicose mode is unstable it grows slower than the kink mode. It is the kink mode that causes a cataclysmic reconfiguration of the topology of the magnetic field. The linear analysis regarding the stability of the varicose mode is also unclear. Tayler (1957) showed the varicose mode is unstable in cylindrical geometry, however utilising a system that is not equivalent to a full neutron star. Lander & Jones (2011) looked at the stability of poloidal fields in neutron stars by way of numerical evolutions of the linearized equations, however did not identify the difference between the varicose and kink modes. Summarily, we believe the stability of the varicose mode to be an open question. Figure 2d shows the simulation after 195 ms = 39 t A . This is a typical snapshot many Alfvén timescales after the non-linear saturation of the unstable modes. We discuss the end-state in more detail in section 3.2.

Instability Growth-Times
We are interested in extracting instability growthtimes for the individual modes. We define the growthtime for each mode as which we measure during the near exponential growth of the instability. We have performed multiple simulations with different magnetic field strengths to examine the scaling behaviour of the instability. In figure 3 we plot the growth-times as a function of mode number for each model. The scaling behaviour is consistent with what is anticipated, i.e. models with stronger magnetic fields and hence faster Alfvén velocities have faster growing instabilities. This scales almost linearly; the model with average field strengthB 15 = 4.4 (whereB 15 =B/10 15 G) has an Alfvén crossing time of t A = 13.9 ms and growthtime of approximately τ ∼ 12 ms, whereas the crossing time for theB 15 = 25.6 model is t A = 2.4 ms and growthtime τ ∼ 2 ms. Whilst this coarse measurement of growth-times scales correctly with magnetic field strength, for some of our models the individual modes do not scale in the manner necessarily expected. Kruskal & Schwarzschild (1954) and Tayler (1957) showed the instability grows faster for increasing m (our m corresponds to their wavenumber k), akin to what is seen for our two simulations with the weakest magnetic fields, i.e. models withB 15 = 4.4 and 8.5 in figure 3. This further agrees with the analysis of Lander & Jones (2011), who showed growth-times converging soon after the m = 4 mode. We do note that our simulations with stronger magnetic fields show different behaviour in that the growth-time as a function of mode number is almost constant.
We provide three possibilities for the cause of the aforementioned discrepancy. 1) We are possibly seeing here mode coupling between the lower order modes, which may scale with magnetic field strength. A more detailed analysis of this is required. 2) We are extracting a superposition of both the varicose and kink modes, which potentially have significantly different growth-times for the same mode number. 3) Our mode signal is 'dirty', in that when evaluating equation (2), we pick a region of the exponential growth and perform a linear regression (in log space) to determine the growth-time. Choosing a different region of the exponential growth can give significantly varied growth-times. This is displayed by way of error bars on the fiducial model in figure 3, which are of similar order to the size of variation one expects between the modes. We expect this to be an inherent difficulty with extracting such growth timescales from non-linear simulations.

Equilibrium Configurations
We now return to figure 2d, particularly concerning our equilibrium configurations. In some ways, figure 2d is consistent with the 'twisted-torus' configurations seen in the non-linear evolutions of Braithwaite & Nordlund (2006) and Braithwaite (2009), and in the semi-analytic equilibrium derivations of  and Ciolfi et al. (2009). This figure shows about half of the star is well approximated by a twisted-torus, however the remainder of the star also exhibits large nonaxisymmetric structures. This is typical of the evolution following the saturation of the instability. Whilst Braithwaite (2008) has also found non-axisymmetric equilibria as a result of various non-linear simulations, we believe ours are the first such simulations to do so with a barotropic equation of state. This therefore extends existing results in the literature (e.g. Ciolfi et al. 2009), in that we provide evidence for a new branch of stable, equilibrium solutions with a barotropic equation of state.
A way of classifying variation from the initial state is through calculating the relative energy in the toroidal and poloidal components of the magnetic field. The initial field obviously contains E p /E = 1.0, where E p and E are respectively the poloidal and total magnetic field energy's. Throughout the evolution of the instability we see a transfer of energy between the toroidal and poloidal energies, although it is interesting to note that this does not take place until we visually see the onset of the kink instability (i.e. ∼ 50 ms). Our equilibrium configuration is then characterised by a magnetic field with E p /E ∼ 0.75, which is similar to the upper-end of the stability window found in stably stratified models of Braithwaite (2009).
It can further be seen from figure 2d and the full evolution of the system that, despite large-scale restructuring of the interior field, the center of the star is still threaded by a dominantly poloidal component. Moreover, the magnetic field exterior to the star is also dominantly poloidal, with little visual change to its initial structure. This is also important for stability issues associated with poloidal fields, particularly the Flowers & Ruderman (1977) instability which relies on minimising the energy in the exterior region of the star (see recent result of Marchant, Reissenegger & Akgün 2010).

CONCLUSION
We have presented the first, non-linear GRMHD simulations of hydromagnetic instabilities inherent to purely poloidal magnetic fields. In particular, we have evolved initially self-consistent solutions of the Einstein-Maxwell field equations under the assumption of ideal MHD and the Cowling approximation. Purely poloidal magnetic field configurations are particularly vulnerable to two modes -the varicose/sausage and kink modes. The varicose mode, a cross-sectional change in a toroidal tube of magnetic field around the neutral line, was initially excited in our simulations due to the Cartesian grid. This mode evolved to a pseudo-equilibrium (figure 2b) before it was swamped by the kink instability (figure 2c), which acts perpendicular to the gravitational field. The kink instability saturated after approximately 15 Alfvén crossing times, at which point the system evolved towards a new, non-axisymmetric equilibrium (figure 2d).
The equilibrium end-states achieved from our simulations partially resemble twisted-torus configurations seen from other non-linear evolutions (Braithwaite & Nordlund 2006;Braithwaite 2008). The center of the star is threaded by a dominantly poloidal component, which is consistent with the instability dominating near the neutral line of the star. In some regions of the star we see the presence of twisted flux tubes, with other regions displaying highly non-axisymmetric portions of the magnetic field. These are the first non-axisymmetric equilibria attained with a barotropic equation of state. An interesting extension to this work is to study the evolution of the exterior magnetic field as a result of the interior rearrangement. Observations of the spin down of neutron stars only infers the dipole component of the field, implying the presence of extra degrees of freedom can only be inferred at this stage through detailed modeling. However, a better understanding of the dynamics of the stellar surface in our simulations is required to accomplish this goal.
There are many further extensions of this work worth pursuing towards more realistic neutron star models.
Such models should feature composition stratification, which is likely to have an impact on MHD equilibrium configurations (e.g.
The end-state quasi-equilibria models discovered herein add to the possible solution space of barotropic equilibria.
Axisymmetric equilibrium models have been explored by Ciolfi et al. (2009);Ciolfi, Ferrari & Gualtieri (2010), which would be interesting to test for stability. Indeed, any pseudo-stable equilibrium could serve as possible starting points for studying the secular magnetic field evolution in magnetars and its connection to activity seen in these objects. Moreover, the theoretical modeling of magnetar QPOs (which has so far been done using symmetric magnetic field configurations), should be extended to account for the non-axisymmetric configurations akin to those encountered here. We hope to address some of these exciting issues in the near future.
This work is supported by the Transregio 7 'Gravitational Wave Astronomy', financed by the Deutsche Forschungsgemeinschaft (DFG). PL and KG are supported by the Alexander von Humboldt Foundation. We are grateful to Jerome Novak and Erik Schnetter for assistance with Lorene, and particularly grateful to Jon Braithwaite for useful discussions.