# Poloidal-Field Instability in Magnetized Relativistic Stars

## Abstract

We investigate the instability of purely poloidal magnetic fields in nonrotating neutron stars by means of three-dimensional general-relativistic magnetohydrodynamics simulations, extending the work presented in Ciolfi et al. (2011). Our aim is to draw a clear picture of the dynamics associated with the instability and to study the final configuration reached by the system, thus obtaining indications on possible equilibria in a magnetized neutron star. Furthermore, since the internal rearrangement of magnetic fields is a highly dynamical process, which has been suggested to be behind magnetar giant flares, our simulations can provide a realistic estimate of the electromagnetic and gravitational-wave emission which should accompany the flare event. Our main findings are the following: (i) the initial development of the instability meets all the expectations of perturbative studies in terms of the location of the seed of the instability, the timescale for its growth and the generation of a toroidal component; (ii) in the subsequent nonlinear reorganization of the system, % of magnetic energy is lost in few Alfvén timescales mainly through electromagnetic emission, and further decreases on a much longer timescale; (iii) all stellar models tend to achieve a significant amount of magnetic helicity and the equipartition of energy between poloidal and toroidal magnetic fields, and evolve to a new configuration which does not show a subsequent instability on dynamical or Alfvén timescales; (iv) the electromagnetic emission matches the duration of the initial burst in luminosity observed in giant flares, giving support to the internal rearrangement scenario; (v) only a small fraction of the energy released during the process is converted into -mode oscillations and in the consequent gravitational-wave emission, thus resulting in very low chances of detecting this signal with present and near future ground based detectors.

###### Subject headings:

stars: neutron — gravitational waves — magnetohydrodynamics (MHD) — methods: numerical## 1. Introduction

Neutron stars (NSs) are endowed with very intense, long-lived, large-scale magnetic fields, reaching strengths which are estimated to be of the order of G at the magnetic pole for ordinary NSs, and around G in the case of magnetars. Such extreme magnetic fields play a crucial role in the physics of NSs, affecting their structure and evolution. They are involved in the processes through which NSs are observed, like the pulsar magnetic dipole radiation and the magnetically-powered burst activity of magnetars, and they have been recently recognised as essential to explain the quasi-periodic oscillations detected in the aftermath of magnetar giant flares [see, e.g. , Gabler et al. (2012) and references therein]. Moreover, they are responsible for deformations which may cause a significant emission of gravitational waves (Bonazzola & Gourgoulhon, 1996; Cutler, 2002) and precession (Wasserman, 2003), they influence the thermal evolution of the star (Pons et al., 2009), to list a few.

All these processes depend on the magnetic field configuration inside the NSs, whose geometry is basically unknown. From observations of the spindown in pulsars the exterior magnetic field appears to be purely poloidal and mainly dipolar, but substantially different internal geometries can reproduce this external appearance. The importance of obtaining such information has motivated a significant effort in studying possible equilibrium models of magnetized NSs, at first with simple geometries, e.g. purely poloidal or purely toroidal fields, and recently with mixed poloidal-toroidal fields. The latest models, built in Newtonian and general-relativistic framework, include Tomimura & Eriguchi (2005); Lander & Jones (2009, 2012); Ciolfi et al. (2009, 2010); Fujisawa et al. (2012), where the so-called ‘twisted-torus’ configuration is considered. This particular geometry has been found as a result of the evolution of initial random fields in Newtonian magnetohydrodynamic (MHD) simulations by Braithwaite & Nordlund (2006).

Once a magnetic-field geometry is chosen, building a corresponding equilibrium configuration is not sufficient to assess whether this represents a good description of the NS interior. The magnetic field, in fact, should also be long-lived and thus stable on timescales which are much longer than the dynamical timescale. Assessing the stability of a given magnetic field configuration is not trivial and most of the work done on the subject concerns simple field geometries and nonrotating stars. Until very recently, the problem has been only addressed with a perturbative analytic approach. These calculations established that both a purely poloidal field and a purely toroidal field are unstable in nonrotating stars, giving important predictions about the onset of the instability, but they could not predict the following evolution of the system. Only recently, thanks to the progress in numerical simulations, it has become possible to study these hydromagnetic instabilities by performing the fully three-dimensional (3D) MHD evolutions of magnetized relativistic stars. These simulations represent a very powerful tool, allowing to confirm the predicted features of the instability and to obtain information about the nonlinear dynamics of the process. In addition, the end-state of simulations can provide important hints about the preferred magnetic field configuration in magnetized stars.

There is an additional and important motivation for studying hydromagnetic instabilities in NSs. The induced global rearrangement of magnetic fields is a violent, strongly dynamical process, and soon after the magnetar model was proposed (Duncan & Thompson, 1992), this kind of process was suggested as a trigger mechanism of giant flares (Thompson & Duncan, 1995, 2001). Nowadays, this internal rearrangement scenario represents one of the leading models to explain the phenomenology observed in magnetars, the other one involving a large-scale rearrangement of magnetic fields in the magnetosphere surrounding the star (Lyutikov, 2003, 2006; Gill & Heyl, 2010). Since both mechanisms can be present in a giant flare, the main question becomes whether most of the magnetic energy powering the flare is stored inside the star or in its exterior magnetosphere. The basic tests on these models rely on the comparison of the predicted timescales and energies involved with the giant flare observations. Hence, determining self-consistently the dynamics associated to this kind of instability can provide important hints on the underlying mechanism.

Moreover, magnetar flares (and in particular giant flares) are likely to be accompanied by a significant excitation of NS oscillations, in particular in the -mode, which can then lead to a strong emission of gravitational waves (GWs). This possibility has motivated recent searches for GWs in connection to magnetar flares, published by the LIGO and Virgo collaboration [see, e.g. Abadie et al. (2011)]. Semi-analytic efforts have been devoted to establishing the maximum amount of magnetic energy released in a magnetar flare, which, in turn, provides an upper limit on the energy emitted in GWs (Ioka, 2001; Corsi & Owen, 2011). These upper limits are based on analytical calculations and simplified models, and can only provide rough, order-of-magnitude estimates. The assumption that all the available energy (which is at most of the order of the total magnetic energy) is converted into GWs, leads to the optimistic conclusion that the signal would be detectable with the next-generation ground-based detectors (Corsi & Owen, 2011). This result has been questioned in Levin & van Hoven (2011), where a simple perturbative analysis is employed to show that only a small fraction of the magnetic energy involved in a giant flare event can be actually be converted into -mode oscillations and that the consequent GW emission is expected to be very weak. Again, referring to the internal magnetic field rearrangement scenario of giant flares, MHD simulations of hydromagnetic instabilities can provide a realistic picture of the GW signal produced, together with estimates of the fraction of energy which can be pumped into the -mode and of the signal detectability.

In this work we focus on the instability of purely poloidal fields,
with the goal of shedding some light on all of the points made
above^{4}

The organization of the paper is as follows. In Sect. 2 we reconsider the setup of the system, improving in particular our treatment of the atmosphere. Within the new setup, we confirm that our evolutions meet all the expectations on the onset of the instability, in agreement with the previous perturbative studies. In Sect. 3.1 we examine in more detail the nonlinear rearrangement of the system, performing much longer simulations and gaining new substantial insight on the final state reached by the system. Section 3.2 is dedicated to a discussion of the implications for the most-likely magnetic-field configurations in magnetized NSs and of the role played by magnetic helicity. We also study the emission properties of the system, relevant for the internal field rearrangement scenario of giant flares, estimating in Sect. 3.3 the timescale of the process and its electromagnetic luminosity, and discussing the detectability of the GW signal in Sect. 3.4. Both for electromagnetic and GW emissions, our estimates rely on a good agreement with the dependence on the magnetic field strength expected in the perturbative limit of weak magnetic fields, which allows us to extrapolate our results to lower and more realistic values than those actually considered in the simulations. Our conclusions are finally presented in Sect. 4. Unless specified differently, we adopt units in which , .

## 2. Physical system and numerical setup

Our physical system of interest is represented by a nonrotating
isolated neutron star surrounded by vacuum, initially endowed with a
purely poloidal magnetic field permeating the star and extending to
the exterior. The initial configurations are fully-relativistic
self-consistent solutions generated with the multi-domain
spectral-method code LORENE, developed at the Observatoire de
Paris-Meudon (Bocquet et al., 1995). The stars are modeled as composed of
a barotropic fluid obeying a polytropic equation of state , with and . The reference
unmagnetized star has a gravitational mass of
and a radius of 12.2 km. Magnetic field
strengths at the pole vary in the range^{5}

We perform fully 3D general-relativistic MHD simulations of the system adopting the Cowling approximation, i.e. neglecting changes in the spacetime metric. Evolutions are obtained with the WhiskyMHD code. Most of the details on our mathematical and numerical setup are discussed in depth in Pollney et al. (2007); Giacomazzo & Rezzolla (2007); Giacomazzo et al. (2009, 2011). One aspect worth stressing is that, to guarantee the divergence-free character of the MHD equations, we use as an evolution variable the vector potential instead of the magnetic field, as described in Giacomazzo et al. (2011). Our standard grid setup consists of three refinement levels using the Carpet driver (Schnetter et al., 2004), with the finest one covering the entire star and having resolution ( m). Our computational box extends to km and we impose no symmetries.

To shorten the time for the development of the instability and thus reduce computational costs, a small perturbation is added to the initial velocity of the fluid. In particular, we add a -component of fluid velocity in the region surrounding the neutral line, where the instability is expected to be triggered, and with azimuthal distribution. The strength of the perturbation corresponds to relative changes in the magnetic field of the order of . We have verified that the instability occurs even in absence of initial perturbations and that no appreciable differences (besides that of reducing computational costs) are introduced in the dynamics by the perturbation.

In our numerical setup, special care is paid to the treatment of the atmosphere. We recall that, as customary in finite-volume relativistic hydrodynamics simulations, our star is surrounded by a fluid at much lower densities (i.e. the “atmosphere”), obtained by imposing a minimum rest-mass density [see, e.g. Baiotti et al. (2005) and Baiotti et al. (2008)]. In the atmospheric region, where the density is equal to the imposed floor value, we set the fluid velocity to zero to avoid the spurious accretion of the atmosphere onto the star. Within the common assumption of ideal MHD, also adopted in the WhiskyMHD code, this prescription would imply that magnetic fields do not evolve in the atmosphere, which represents a serious limitation from both the physical and numerical point of view. In particular, this rapidly leads to errors at the stellar surface and, for the magnetic field strengths we consider, ends the simulations prematurely.

A self-consistent solution to this problem would require the implementation in general-relativity of the equations of resistive MHD, along the lines of the work reported by Palenzuela et al. (2009) and which has seen recent progress in the work of Dionysopoulou et al. (2012). Lacking for the time being this more systematic solution, a reasonable first approximation can be obtained by adding a resistive term to the evolution equation for the vector potential, i.e.

(1) |

where [see Antón et al. (2006) for the notation], is the scalar resistivity and, for simplicity, we take to represent the Laplacian in a flat spacetime. We impose that resistivity is zero inside the star (thus retaining an ideal-MHD behaviour), up to a transition layer, where it continuously increases towards the atmospheric value . This is shown in Fig. 1, where we plot the initial resistivity profile along the radial direction for our fiducial simulation with . The plot shows two different radial profiles, where we change the width of the transition region inside the star, whose surface is indicated with a vertical solid line. As we will discuss below, the two choices have a different impact on the overall results. The resistivity profile is set to be an explicit function of the rest-mass density, namely , where is a Fermi-like function, with in the atmosphere and inside the stellar core. In this way, the resistivity will mimic the evolution of the rest-mass density, following the stellar surface in its evolution.

Since the velocity is set to zero outside the star, eq. (1)
reduces there to . A few
remarks are worth doing about this limit and we start by considering
the special relativistic case for simplicity. Maxwell’s equations in
vacuum are simple wave equations, e.g. , so that a non-zero Laplacian of the magnetic field
would be simply radiated away in electromagnetic waves. In the
quasi-static limit, that is, in the limit in which the timescale
associated with magnetic-field variations is large compared to the
light travel time (where is the typical lenghtscale of field
variations), any episodic time variation in the magnetic field
produced by the dynamics in the stellar interior would be rapidly
radiated away, leading to an exterior magnetic field which is again
with . In our system, where the quasi-static
limit represents a good approximation, a similar behaviour outside the
star can be obtained by evolving the magnetic field according to a
diffusion equation^{6}

This logic is valid as long as the external evolution takes place in the quasi-static limit, and we fulfil this requirement by suitably choosing the timescale of magnetic diffusion via the choice of the resistivity . More precisely, we set to be always high enough (hence with sufficiently short associated resistive timescales) so as to keep the exterior of the star always close to the condition . Although only an approximation, the approach discussed above allows us to have a dynamical exterior magnetic field, which adjusts itself to the changes triggered in the interior by the development of the instability. The price to pay is that, in addition to the one lost because of the numerical resistivity, the magnetic energy of the system is not conserved. Per-se, this would not be particularly problematic, since in a realistic system one would expect the magnetic energy to be radiated away, but it does introduce a dependence of the problem on the profile chosen for the resistivity, since different profiles will be responsible for energy losses.

As we will discuss in more detail the following Section, in fact, during the initial development of the instability, the magnetic-field modifications are still very small compared to the background field, and the system is not expected to show visible changes, with the magnetic energy being essentially conserved. Yet, because of the resistive layer inside the star introduced via eq. (1), magnetic-energy losses will be present even in the early stages of the evolution, simply because the magnetic field will be divergence free but not Laplacian-free in the outer layer of the star.

In our previous work (Ciolfi et al., 2011), we had chosen a rather wide resistive transition layer (cf. red dashed line in Fig. 1), which resulted in significant magnetic-energy losses from the very beginning of the simulation. To minimize these losses, we have here considered a thinner resistive transition layer (cf. blue solid line in Fig. 1), thus involving a smaller stellar volume of the star where resistivity can act. As a result, the magnetic-energy losses during the exponential-growth phase (i.e. during the first ms) go from in the case of a wide layer, to below in the case of a thin layer. In the subsequent stages of the instability, the differences between the two prescriptions are much smaller and this is because the dynamics of the field is much less influenced by the properties of the resistive layer at the surface of the star. Yet, since we have shown that different layers lead to different energy losses, it is reasonable to wonder whether the properties of the transition region can be important also for the subsequent evolution of the system. As we will argue in detail in the Appendix, results are effectively independent of as long as a suitable value for is chosen for the different magnetic field strengths; more specifically, if the resistivity is chosen to scale linearly with as (see Appendix).

## 3. Results

### 3.1. Poloidal field instability

We next describe the evolution of the system for our fiducial set of
initial models, from the onset of the instability to the nonlinear
rearrangement of magnetic fields. First of all, we briefly recall the
basic expectations of the perturbative studies on
the onset of the instability for a purely poloidal field in nonrotating
magnetized stars (Markey & Tayler, 1973; Wright, 1973): (i) the instability first
develops in the region of closed-field lines surrounding the neutral
line; (ii) toroidal magnetic fields are generated in this region and
grow exponentially until their local intensity is comparable to the
poloidal one, which corresponds to the saturation of the instability;
(iii) the instability saturates in about one Alfvén time^{7}

A first visual overview of the dynamics of the system is given in
Fig. 2, where we present snapshots of the star and its
magnetic field in the meridional (top row) and equatorial
planes (bottom row) at three different stages of the
evolution. In this example . The left panels refer
to the time of saturation of the instability, when the magnetic field
starts to show visible modifications ( ms). From the equatorial
view we note that the initial axisymmetric geometry is lost and
replaced by a non-axisymmetric structure. Around 7.5 ms (central
panels) the instability has fully developed. As expected, the closed
line region is filled with toroidal fields of strength comparable to
the background, resulting in vortex-like structures of magnetic field
lines. This stage is the most dynamical and the rapid modifications of
the field lead to the expulsion of matter from the surface of the star
( in this example). Up to this point in
the evolution the instability has affected only the region of closed
magnetic field lines, in accordance with the perturbative predictions.
However, as the nonlinear rearrangement of the magnetic field
proceeds, the whole star is involved, with changes encompassing also
the open field lines and the exterior field. Of course, the extent of
these modifications will depend on the strength of the magnetic field,
being larger for more violent field dynamics. The last panels refer to
the end state of our simulation ( ms). At this stage there is no
trace of the initial geometry and the system has lost most of its
magnetic energy (see discussion below). To make sure that the external
magnetic field is always close to a potential one (even if its
geometry is unusual, as in this case) we have computed its curl and
found it is always smaller than , while inside the star
and near the surface it reaches . The same result holds
for all the cases considered in this work^{8}

In Fig. 3 we plot the evolution of the poloidal,
toroidal and total magnetic energies. In practice, the evolution of
the total magnetic energy (left panel) provides information on the
various stages of the evolution: the initial phase up accompanied by
very small variations, the sudden drop of magnetic energy with a loss
of % of magnetic energy in few ms to tens of ms depending on
the initial magnetic field strength, and the final stage with a slower
evolution as the energy decreases to few percent of the initial value
and the system reaches a quasi-stationary equilibrium^{9}

The overall evolution does not change with smaller initial field strengths, apart from the expected increase of timescales. Figure 3 also reports the poloidal and toroidal field energies in a logarithmic scale (middle panel). Note that the the toroidal energy grows exponentially in the initial stages of the evolution and that the growth rate depends on the initial magnetic-field strength, the slope being steeper for larger magnetic fields. The poloidal component, on the other hand, remains almost unchanged during this initial stage, until the nonlinear rearrangement of the magnetic field starts. Then, the system evolves losing most of the poloidal energy, while the toroidal energy experiences smaller variations, thus resulting in a growing ratio of toroidal and poloidal energies. In the last phase of the evolution, the ratio of the poloidal-to-toroidal magnetic energies stabilizes in the range , indicating that at this stage the toroidal magnetic field provides a dynamically important contribution to the final quasi-stationary balance (right panel of Fig. 3).

From a closer look at the details of the dynamics, we notice a
qualitative difference between stronger and weaker magnetic fields,
the latter having a much smoother evolution. For the system undergoes indeed less dramatic
modifications, where the initial magnetic field geometry is partially
maintained^{10}

In Fig. 5 we show instead the growth time of the toroidal field energy in the exponential phase of the evolution, as a function of the initial magnetic field strength. The predicted linear scaling is very well satisfied in the full range of magnetic field strengths considered (the red dashed line represents a linear fit to the data).

### 3.2. Final magnetic field configuration

Since the new simulations have been carried out on timescales which are much longer than the ones investigated in Ciolfi et al. (2011), we are in a much better position to discuss the properties of the final magnetic-field configuration. Also important in this context is the evolution of purely-hydrodynamical quantities, which can provide important additional information on the development of the instability and on the quasi-stationary state approached by the system in the final stages. One of such quantities is the central (i.e. maximum) rest-mass density. We recall that the initial purely poloidal magnetic field is responsible for a quadrupolar deformation of the star, whose resulting shape is slightly oblate (i.e. with a positive ellipticity ), despite being nonrotating. This deformation lowers the central rest-mass density of the star with respect to the unmagnetized case, and the difference is of the order of . Since the ellipticity scales as , when the system loses a factor in poloidal magnetic energy, should be also reduced by about the same factor. In addition, if there is a significant toroidal component, this tends to deform the star in a prolate shape, contrasting the effect of the poloidal field. As a consequence, in our simulations we expect the final to be much smaller than the initial one, i.e. the star basically recovers the spherical shape of the unmagnetized case.

In Fig. 6 we show the evolution of the central rest-mass density for different initial magnetic field strengths. Note that our different models are all built with the same initial central density, which gives a percent difference in the central density of the unmagnetized star corresponding to each model. Clearly, also for the rest-mass density the evolution is much more dramatic for stronger initial magnetic fields, with excursions in the central density that become larger and more rapid as the initial magnetic field strength is increased. In all cases, however, the central density mimics the evolution of the magnetic energy and reaches an approximately constant value, with an overall jump which scales roughly as . The absence of additional evolution in the rest-mass density after the instability has fully developed provides the indication that the new configuration is, at least hydrodynamically, stable, and that, if the new magnetic-field configuration is about to lead to a new unstable evolution, it will do so on much larger timescales. For initial magnetic field strengths of and , the change is so violent and rapid that it leads to high-frequency oscillations in the central density, while we do not observe the same effect with smaller fields. As we will further discuss in Sect. 3.4, a rapid variation of central density can also affect the emission properties of the star, e.g. introducing a significant modulation of the gravitational-wave signal.

Other useful information on the final state can be obtained by looking
at the evolution of the total magnetic helicity . This quantity
is associated with the topological properties of a given magnetic
field geometry, and measures the degree of linkage of magnetic field
lines (Moffatt, 1969). The initial purely poloidal configuration
has because only regions with a mixed poloidal-toroidal field
contribute to the magnetic helicity^{11}

We recall that although in ideal MHD the total magnetic helicity is conserved (Woltjer, 1958), a highly-conducting fluid star in vacuum with a thin resistive layer does not represent an ideal-MHD system, and magnetic helicity conservation is therefore not guaranteed [of course helicity is also produced because of the intrinsic nonzero numerical resistivity, but the latter is much smaller than the one introduced in eq. (1)]. We compute the total magnetic helicity as

(2) |

where is the magnetic helicity 4-current and is the spatial hypersurface at a given time [see Bekenstein (1987) and references therein]. Since our initial helicity is , we normalize to the helicity of an axisymmetric twisted-torus configuration in which the poloidal field is arranged as in our initial condition and the toroidal field, confined in the region of closed field lines around the neutral line, is uniform and with an energy equal to the poloidal one, i.e. . Indicating with such reference helicity, we can estimate it as , where (8 km) is the radius of the neutral line and is set to coincide with the magnetic energy of the stellar model under exam. Note that we are not concerned here with the sign of (or of ), since the latter depends on whether toroidal fields are directed along the positive or negative -direction and a global transformation would not change the properties of the system.

We can now monitor the evolution of . As long as the total magnetic helicity is small, while is an indication that there is a significant magnetic helicity. In Fig. 7 (top panel) we show the evolution of the normalized helicity for the different initial magnetic field strengths. Note that in all cases grows in time, and around 60 ms it has reached a substantial fraction of 1. This fraction is higher for weaker fields, where the evolution is smoother, i.e. for , while for . This behaviour supports the idea that configurations with a significant amount of magnetic helicity are more stable, as already suggested in the literature [see Braithwaite & Spruit (2006) and references therein]. Such indication is compatible with the observed instability of the initial purely poloidal field and the production of a toroidal component, which could be interpreted as an attempt of the system to produce magnetic helicity.

We have seen that the system can achieve a significant amount of magnetic helicity, but in the meanwhile most of its magnetic energy is lost. In order to weigh the absolute variation of magnetic helicity, we can compare with the initial value of (bottom panel of Fig. 7). The ratio in this case is much smaller, , indicating that the magnetic helicity produced would represent a small amount for the initial star, and becomes significant only because the system has lost most of its magnetic energy.

On the basis of these results and bearing in mind the limits of our approach, we conjecture that magnetic helicity could play an important role as a stabilizing element for a fluid magnetized star. If this is the case, we can give a natural interpretation of the evolution of our system. (i) A purely poloidal magnetic field is unstable and a significant modification of its topology is necessary for stabilization. (ii) While the system is evolving the rate of change of magnetic helicity is small compared to the one of magnetic energy, thus the system has to lose most of its magnetic energy in order to significantly alter its topology and reach a more stable state.

As a complement to the information obtained from our nonlinear relativistic calculations, we should remark that the possibility of a stable equilibrium has been questioned in the case of a barotropic fluid star (Reisenegger, 2009; Lander & Jones, 2012). A barotropic equation of state does not take into account the effect of stable stratification associated to composition gradients, which could play an important role in determining the stability of magnetic equilibria in relativistic NSs (Reisenegger & Goldreich, 1992; Reisenegger, 2009). Stable stratification, indeed, offers an additional degree of freedom which favours the hydrostatic balance of fluid and magnetic forces. A strong indication in support of stable equilibria in stratified magnetized stars has been obtained for main-sequence stars, where stratification is provided by entropy gradients (Braithwaite & Nordlund, 2006), while there is no such evidence for barotropic NSs.

### 3.3. Electromagnetic emission

As discussed in the Introduction, one of the leading explanations of magnetar giant flares assumes that the event is due to a sudden, global rearrangement of the internal magnetic field, i.e. the same kind of process that we are studying. A giant flare starts with an initial burst of duration, followed by a long pulsating tail lasting hundreds of seconds. In such scenario, the energy emitted in the initial burst is a significant fraction of the total magnetic energy of the star, while the duration of the burst is dictated by the timescale of the internal field rearrangement. We have also remarked repeatedly that in all cases considered, our magnetized stars lose most of the initial magnetic energy within few tens of ms. This energy is converted in very small part into fluid motions, e.g. -mode oscillations of the star with consequent emission of gravitational waves (see Sect. 3.4), while most of it is dissipated in the resistive layer inside the star and then in the atmosphere, thus mimicking the radiative losses that would be measured if the stars were actually in vacuum. We can therefore use the information about the dissipated magnetic energy to deduce, within the approximation of our approach, an order-of-magnitude estimate of the electromagnetic luminosity associated with the hydromagnetic instability. In addition, we can compare the luminosity and duration of the emission computed in this way with the observations of giant flares.

In practice, most of the emission comes from the initial fast drop of magnetic energy by , corresponding to an initial spike in . Taking that as a reference, we measure the duration of the spike and compute its peak and (time) average luminosities for the different initial magnetic field strengths considered. Note that the expected scaling for the luminosity is , since the timescale of the emission goes as , while the magnetic energy scales as . Our estimate for the average luminosity, obtained by fitting for the different magnetic field strengths, gives

(3) |

Similarly, for the reference value of a magnetar-type magnetic field of , we obtain a duration of the initial spike of and a peak luminosity of .

We can now compare with the observations. In the brightest giant flare detected, the one in 2004 from SGR 1806-20, the initial spike lasted , had a peak luminosity of , where is the distance of the source, and accounting for the total isotropic energy radiated (which was % of the total energy emitted in the whole giant flare event, including the s tail), the average luminosity was [see Mereghetti (2008) and references therein]. The duration of the initial burst from SGR 1806-20 is clearly comparable with our estimate. The difference in average luminosity (our estimate is larger), then, essentially reflects a difference in the electromagnetic energy released, which does not constitute a significant limitation for the scenario: while in our system the magnetic energy is almost completely lost, it is perfectly reasonable that in a more realistic situation the magnetic field would migrate from one configuration to another with a smaller jump in magnetic energy. Note that the timescale of the process, on the other hand, is essentially controlled by the initial internal magnetic field strength and poorly depends on the overall jump in magnetic energy. Interestingly, also the ratio between the estimated peak and average luminosities is in good agreement with the observations.

If the comparison with the phenomenology of SGR 1806-20 leads to a reasonably good match in the duration and energetics, thus providing a substantial support to the internal field rearrangement scenario, the comparison with other observations is not as striking. In particular, the average luminosity of the initial spike in the case of the other two giant flares observed (1979 from SGR 0526-66 and 1998 SGR 1900+14) was , with a comparable duration of . In these cases, therefore, the mismatch in the energy losses is far larger, although still acceptable if the field rearrangement is to a new configuration with comparable magnetic field strengths.

The luminosity rise time at the beginning of the initial spike constitutes an additional relevant timescale in a giant flare. Being of the order of ms (Palmer et al., 2005), it can be easily associated with the Alfvén propagation time in the magnetosphere, while it can be hardly explained by the internal magnetic field rearrangement, which acts on longer timescales. Therefore, in order to have the internal scenario fully compatible with the observations one probably needs to assume that the instability also triggers an initial, less energetic but sudden reorganization of the magnetospheric field, with a mechanism similar to the one proposed in the alternative external rearrangement scenario (Lyutikov, 2003, 2006; Gill & Heyl, 2010). However, this feature can not be captured by our present modelling of the NS exterior.

In summary, although the dynamics of our system is oversimplified when compared to the complex phenomenology shown by a realistic magnetar giant flare, the overall agreement in the duration of the burst and in its energetics when compared with the observations of the giant flare from SGR 1806-20, lends support to the suggestion that the basic phenomenology is that of an internal-field readjustment. To the best of our knowledge, this is the first time that a comparison between general-relativistic MHD simulations and the magnetar phenomenology has been attempted. Clearly, a lot more work needs to be done to improve our self-consistent but crude modelling.

### 3.4. Gravitational wave emission

In the Introduction we have pointed out the potential relevance of our study for GW observations in connection to magnetar giant flares. The basic idea, already presented in Ciolfi et al. (2011), is that the instability triggers stellar oscillations at the -mode frequency, with the consequent emission of GWs, similarly to what should happen in association with a magnetar giant flare. In what follows we provide a systematic assessment of this emission for the different initial magnetic field strengths considered and an estimate of the detectability of the GW signal for the planned advanced GW detectors, i.e. advanced LIGO and advanced Virgo, and the new generation ones, e.g. Einstein Telescope (Punturo et al., 2010).

The GW signals produced in our simulations are shown in the different panels of Fig. 8, where we report the amplitudes in the and polarizations for the various initial magnetic field strengths , ranging from to . Besides the differences in the overall amplitude, that naturally grows with (see also discussion below), the waveforms show also other differences with the magnetic field strength. These include variations in the transient stage of the instability and low-frequency modulations. For example, for we can observe a significant modulation at , which could be associated with the jump in central density (cf. Fig. 6), whose timescale is compatible with the period of the modulation. A similar modulation appears also for . Rather irregular transients characterize instead the initial amplitudes for and . In both cases, in fact, we can notice sudden jumps in the signal, which correspond to the highly dynamical phases in the evolution of the system. This is indeed confirmed by the evolution of the total magnetic energy in Fig. 3, where bumps appear around and for and , respectively, while they are not present in the other cases.

Despite these differences, all signals are essentially dominated by oscillations at the -mode frequency, as it is clear from Fig. 9, where we report the Fourier transform of the GW amplitude in both the and polarization for and . The -mode we can read off the plot is about for the lowest field strength, with a shift to lower frequencies for higher magnetic fields (this shift is of for ). We recall that since our evolutions are in the Cowling approximation, the corresponding -mode frequencies are notoriously larger of than the correct ones computed in full general relativity (Dimmelmeier et al., 2006; Takami & Rezzolla, 2011). This property will be taken into account when considering the detectability of the signal.

Given the complexity of the signal, with large rapid and secular variations, the determination of the signal-to-noise ratio (SNR) cannot be done by simply looking at the maximum-minimum amplitudes, but rather by computing a strain which represents a suitable time average. This is indeed what can be done through the root-sum-square amplitude , defined as

(4) |

where . The amplitudes in the two polarizations are computed considering a source within the Galaxy, i.e. at a fiducial distance of . If, as in our case, the signal is dominated by a single frequency , then the SNR can be estimated using the simple expression, , where is the strain-noise amplitude of the detector at the frequency . When computing we also need to specify the duration of the signal (in addition to the distance of the source). Typical estimates of the -mode damping time range between and (Andersson & Kokkotas, 1998; Benhar et al., 2004), and here we assume .

The resulting amplitudes are reported in Fig. 10, where we also show with dotted lines the strain-noise amplitude of advanced LIGO and Virgo and of the Einstein Telescope, and where we have assumed to correct for the Cowling approximation. Furthermore, the “noise” is estimated to be for advanced LIGO and and advanced Virgo, and for the Einstein Telescope [see Andersson et al. (2011), where the same strain-noise amplitudes were considered].

Because in our simulations we have considered magnetic-field strengths that are at least an order of magnitude larger than those typically associated with magnetars, we need to rescale our results for to smaller field strengths in order to assess the GW signal that could be expected if the instability is associated to a magnetar giant flare. This is rather straightforward to do since the expectation from perturbative analyses [see, e.g. , Levin & van Hoven (2011)] is that, independently on the mechanism considered, the amplitude of the excited -mode by the flare should scale as the magnetic energy, i.e. as . This prediction is valid as long as the magnetic field strengths are sufficiently small so that the corresponding magnetic energy can be considered as a small perturbation to the total binding energy. The results shown in Fig. 10 confirm that a quadratic scaling is a rather good approximation to the computed data for . This is highlighted in Fig 10 by the red solid line, which represents a quadratic fit to the GW amplitude in the case of low-magnetic fields. We find it very reassuring that our data does show the expected scaling behaviour in the weak-field limit, which instead appears to be absent in the scaling reported by Zink et al. (2011) and subsequently by Lasky et al. (2012) for rotating stars. Of course, this scaling is then broken as the magnetic field is increased and a steeper dependence is then expected and found for . The characteristic amplitude then appears to saturate for . We note that the very rapid increase in for is somewhat surprising but not completely unexpected. We recall, in fact, that the root-sum-square amplitude refers to a secondary quantity which is calculated in a Newtonian approximation. Because of the high time-derivatives of the fluid variables that contribute to its measure, this quantity is very sensitive to the dynamics of the system and hence to the large velocities that develop for large magnetic fields. By using different prescriptions of the resistivity we have verified that this rapid change is robust, although probably amplified by the breaking of the Newtonian approximation.

The importance of the quadratic fit is that it allows us to extrapolate back to even smaller values of the magnetic fields and obtain the following predictions for the SNRs for the different detectors:

(5) | ||||

(6) |

We conclude that the GW signal produced by the -mode oscillations triggered by the hydromagnetic instability in association with a magnetar giant flare will be undetectable for realistic values of the magnetic field, i.e. and marginally detectable by third generation detectors for unrealistic magnetic fields, . Because a longer damping time, say of , would yield a gain of only a factor , we do not expect that the prospects of a detection would improve if the -mode oscillations would last a factor ten more in time.

Understanding why the detectability of this GW signal is so hard in practice can be made easier if we take into account how much energy is actually lost to gravitational waves and compare it with the amount of dissipated magnetic energy. More specifically, if is the root-sum-square GW amplitude, the energy emitted in GWs can be computed as [see, e.g. Levin & van Hoven (2011)]

(7) |

where we retained factors of and , and is the source
distance^{12}

## 4. Summary and conclusion

We have performed 3D general-relativistic MHD simulations of nonrotating magnetized NSs endowed with a purely poloidal magnetic field and studied the development of the hydromagnetic instability that develops dynamically. This work represents an extension of our previous study on the subject (Ciolfi et al., 2011), where we have improved our treatment of the atmosphere, drastically reducing the undesired energy losses due to the transition between the ideally conducting stellar interior and the resistive exterior, and we have performed simulations on much longer timescales, which have allowed us to gain essential information about the final configuration reached by the system. The resulting overall picture is much clearer and conclusive.

As expected from perturbation analyses, the instability is first triggered in the region of close field lines and is accompanied by the production of a toroidal magnetic field. When the growth of the latter saturates in about one Alfvén timescale, the toroidal field has reached a local strength which is comparable to the poloidal one and the axisymmetry of the initial configuration is lost. At this point, the most dynamical phase of the nonlinear evolution takes place, with major modifications of the magnetic field leading to a strong electromagnetic emission carrying away of the magnetic energy in few Alfvén timescales. At the same time, a small fraction of magnetic energy is converted into stellar oscillations, mostly at the -mode frequency, which cause the emission of gravitational waves. The subsequent evolution proceeds on longer timescales, with further loss of magnetic energy.

Since in the post-instability phase the magnetic field is continuously changing, losing strength because of resistive dissipation, it is difficult to determine whether the corresponding configuration is a stationary one or not. The only robust evidence is that the variations in the hydrodynamical and electromagnetic quantities are much smaller than those during the instability, taking place on larger timescales. We therefore interpret this behaviour as evidence that the new magnetic-field configuration has reached a quasi-stationary state or that, if still intrinsically unstable, the growth time of the instability is much larger than the one that can be possibly investigated numerically.

Because in this quasi-stationary configuration the ratio of the toroidal and poloidal magnetic energies tends to unity, and because the growth of toroidal magnetic field at the expense of the poloidal one is also accompanied by an increase in the magnetic helicity, we are led to conclude that an equilibrium magnetic field configuration with a significant amount of magnetic helicity and comparable poloidal and toroidal magnetic field energies could be a preferred one for stability. Bearing this in mind, we remark that it is still unclear if stable equilibria exist for a simple barotropic fluid star, and that other stabilizing contributions, such as a stable stratification, may be necessary to obtain long-lived magnetic field configurations.

The violent reorganization of magnetic fields induced by the development of a hydromagnetic instability has been proposed as a possible mechanism to explain giant flares in magnetars. Using our simulations and in particular the information about the dissipated magnetic energy, we have deduced, within the approximation of our resistive approach, an estimate of the electromagnetic luminosity associated with the hydromagnetic instability. More specifically, we found that the average luminosity scales as with the magnetic-field strength, in good agreement with the expectation that the radiated energy should scale as , while the duration of the emission should scale as . In this way we were able to perform a direct comparison with the observations of SGR 1806-20, finding a very good agreement with the duration of the burst and an emitted luminosity which is about an order of magnitude larger than the measured one. As a result, although our modelling is oversimplified, the overall agreement with the observations from SGR 1806-20, lends support to the suggestion that the basic phenomenology is that of an internal-field readjustment.

Finally, we have discussed the gravitational-wave emission which should be expected as a result of the -mode oscillations triggered by the instability. Also in this case, our calculations reproduce the expected scaling behaviour of the root-mean-square gravitational-wave amplitude in the limit of weak magnetic fields, i.e. . This important result allows us to extrapolate with confidence our estimates to smaller magnetic fields than those covered by our simulations and conclude that the gravitational-wave signal from -mode oscillations will be undetectable for realistic values of the magnetic field, i.e. and marginally detectable by third generation detectors for unrealistically large magnetic fields, i.e. .

As a self-consistent solution to this problem in fully resistive regime is close to be within reach (Dionysopoulou et al., 2012), we plan to extend the investigation reported above and thus remove many of the approximations that our present analysis had to sustain. It will then be possible to set even more precise connections between the complex phenomenology observed in magnetar flares and the violent dynamics of hydromagnetic instabilities in relativistic stars.

We are grateful to D. Radice, A. Harte, S. K. Lander, S. Mereghetti, W. Kastaun, B. Giacomazzo, E. Bentivegna, G. M. Manca, R. De Pietri and S. Bernuzzi for useful discussions. R. Ciolfi is supported by the Humboldt Foundation. Support comes also from “CompStar”, a Research Networking Programme of the European Science Foundation and by the DFG grant SFB/Transregio 7. The calculations have been performed on the supercomputing clusters at the AEI.

## Appendix A The role of resistivity in the atmosphere

In this Appendix we discuss how the results of our simulations depend on the atmospheric value of resistivity and how a suitable choice for can considerably limit the influence of this free parameter. We recall that if the value of is too small, we expect that the exterior field is not evolving rapidly enough relative to the interior dynamics, thus accumulating magnetic field distortions in the external layer of the star. As is increased, the timescale for the exterior field evolution decreases until it becomes comparable or shorter than the timescale of the internal evolution. At that point the overall dynamics is less sensitive to a further increase of . Finally, if we further increase , it eventually will become too large and the effects of resistivity in the outer layer of the star will become dominant, significantly influencing also the interior dynamics.

In Fig. 11 we show the evolution of the total magnetic energy and of the ratio of toroidal and poloidal energies, obtained with different values of ranging from to , with an initial magnetic field strength . We first focus on the total magnetic energy (left panel of Fig. 11) and note that in the range to , the differences among different evolutions become smaller and smaller as the value of is increased. Furthremore, while a value of still gives comparable results, an additional increase leads again to significant differences. As a result, for these magnetic field strengths, represents a reasonable value for the resistivity.

Let us now consider the toroidal-to-poloidal energy ratio (right panel of Fig. 11), which is particularly sensitive to the value of resistivity. When considering ranging from to , we have a significant increase of the toroidal field production. Around , however, the ratio becomes almost independent of , with no significant differences for between and . However, if we further increase above , we again find a significant dependence, this time with a decrease in the production of toroidal-field.

In conclusion, for a resistivity value of yields optimal results and the scaling of this value for different magnetic field strengths can be done simply as to ). In the limit of low magnetic field, a change in accompanied by this rescaling of would give exactly the same dynamics, only rescaled in time.

### Footnotes

- affiliation: Max-Planck-Institut für Gravitationsphysik, Albert-Einstein-Institut, Potsdam, Germany
- affiliation: Max-Planck-Institut für Gravitationsphysik, Albert-Einstein-Institut, Potsdam, Germany
- affiliation: Department of Physics and Astronomy, Louisiana State University, Baton Rouge, Louisiana, USA
- The instability of purely-toroidal fields has several analogies with the one considered here for purely-poloidal fields and has been investigated by Kiuchi et al. (2008, 2011).
- For simplicity, hereafter we will use a more compact notation for magnetic-field strengths in scales of , i.e. .
- Note that the same equation holds for the vector potential, .
- An estimate of the Alfvén time is given by , where is the average rest-mass density. Hence, ms for our fiducial model with .
- An essentially potential external magnetic field is not surprising given that the field in the stellar exterior is to a good precision with zero divergence and zero Laplacian.
- Note that we refer to a “quasi-stationary equilibrium” and not to a stable-equilibrium because, strictly speaking, the latter is impossible to prove with numerical simulations over a finite amount of time.
- In particular, the open field lines and the exterior field are not strongly affected as in the case of stronger fields (see Fig. 2). This result is more similar to what obtained in Lasky et al. (2011).
- We note, however, that even a mixed field can be constructed such that , e.g. in axisymmetry, through the superposition of a poloidal magnetic field and of two toroidal ones of opposite polarity.
- Note that does not depend on the source distance, while alone does.

### References

- Abadie, J., et al. 2011, ApJ, 734, L35
- Andersson, N., & Kokkotas, K. D. 1998, MNRAS, 299, 1059
- Andersson, N., et al. 2011, Gen. Relativ. Gravit., 43, 409
- Antón, L., Zanotti, O., Miralles, J. A., Martí, J. M., Ibáñez, J. M., Font, J. A., & Pons, J. A. 2006, ApJ, 637, 296
- Baiotti, L., Hawke, I., Montero, P. J., Löffler, F., Rezzolla, L., Stergioulas, N., Font, J. A., & Seidel, E. 2005, Phys. Rev. D, 71, 024035
- Baiotti, L., Giacomazzo, B., & Rezzolla, L. 2008, Phys. Rev. D, 78, 084033
- Benhar, O., Ferrari, V., & Gualtieri, L. 2004, Phys. Rev. D, 70, 124015
- Bekenstein, J. D. 1987, ApJ, 319, 207
- Bocquet, M., Bonazzola, S., Gourgoulhon, E., & Novak, J. 1995, A&A, 301, 757
- Bonazzola, S., & Gourgoulhon, E. 1996, A&A, 312, 675
- Braithwaite, J., & Nordlund, Å. 2006, A&A, 450, 1077
- Braithwaite, J., & Spruit, H. C. 2006, A&A, 450, 1097
- Braithwaite, J. 2007, A&A, 469, 275
- Ciolfi, R., Ferrari, V., Gualtieri, L., & Pons, J. A. 2009, MNRAS, 397, 913
- Ciolfi, R., Ferrari, V., & Gualtieri, L. 2010, MNRAS, 406, 2540
- Ciolfi, R., Lander, S. K., Manca, G., & Rezzolla, L. 2011, ApJ, 736, L6
- Corsi, A., & Owen, B. J. 2011, Phys. Rev. D, 83, 104014
- Cutler, C. 2002, Phys. Rev. D, 66, 084025
- Dimmelmeier, H., Stergioulas, N., & Font, J. A. 2006, MNRAS, 368, 1609
- Duncan, R. C., & Thompson, C. 1992, ApJ, 392, L9
- Dionysopoulou, K., Alic, D., Palenzuela, C., Rezzolla, L., & Giacomazzo, B. 2012, arXiv:1208.3487
- Fujisawa, K., Yoshida, S., & Eriguchi, Y. 2012, MNRAS, 422, 434
- Gabler, M., Cerdá-Durán, P., Stergioulas, N., Font, J. A., & Müller E. 2012, MNRAS, 421, 2054
- Geppert, U., & Rheinhardt, M. 2006, A&A, 456, 639
- Giacomazzo, B., & Rezzolla, L. 2007, Class. Quantum Grav., 24, S235
- Giacomazzo, B., Rezzolla, L., & Baiotti, L. 2009, MNRAS, 399, L164
- Giacomazzo, B., Rezzolla, L., & Baiotti, L. 2011, Phys. Rev. D, 83, 044014
- Gill, R., & Heyl, J. S. 2010, MNRAS, 407, 1926
- Ioka, K. 2001, MNRAS, 327, 639
- Kiuchi, K., Shibata, M., & Yoshida, S. 2008, Phys. Rev. D, 78, 024029
- Kiuchi, K., Yoshida, S., & Shibata, M. 2011, A&A, 532, A30
- Lander, S. K., & Jones, D. I. 2009, MNRAS, 395, 2162
- —. 2011, MNRAS, 412, 1730
- —. 2012, arXiv:1202.2339
- Lasky, P. D., Zink, B., Kokkotas, K. D., & Glampedakis, K. 2011, ApJ, 735, L20
- Lasky, P. D., Zink, B., Kokkotas, K. D. 2012, arXiv:1203.3590
- Levin, Y., & van Hoven, M. 2011, MNRAS, 418, 659
- Lyutikov, M. 2003, MNRAS, 346, 540
- Lyutikov, M. 2006, MNRAS, 367, 1594
- Markey, P., & Tayler, R. J. 1973, MNRAS, 163, 77
- Mereghetti, S. 2008, A&A Rev., 15, 225
- Moffatt, H. K. 1969, Journal of Fluid Mechanics, 35, 117
- Palenzuela, C., Lehner, L., Reula, O., & Rezzolla, L. 2009, MNRAS, 394, 1727
- Palmer, D. M., et al. 2005, Nature, 434, 1107
- Pollney, D., et al. 2007, Phys. Rev. D, 76, 124002
- Pons, J. A., Miralles, J. A., & Geppert, U. 2009, A&A, 496, 207
- Punturo, M., et al. 2010, Class. Quantum Grav., 27, 084007
- Reisenegger, A., & Goldreich, P. 1992, ApJ, 395, 240
- Reisenegger, A. 2009, A&A, 499, 557
- Schnetter, E., Hawley, S. H., & Hawke, I. 2004, Class. Quantum Grav., 21, 1465
- Takami, K., & Rezzolla, L. 2011, MNRAS, 416, L1
- Thompson, C., & Duncan, R. C. 1995, MNRAS, 275, 255
- Thompson, C., & Duncan, R. C. 2001, ApJ, 561, 980
- Tomimura, Y., & Eriguchi, Y. 2005, MNRAS, 359, 1117
- Wasserman, I. 2003, MNRAS, 341, 1020
- Woltjer, L. 1958, Proceedings of the National Academy of Science, 44, 489
- Wright, G. A. E. 1973, MNRAS, 162, 339
- Zink, B., Lasky, P. D., & Kokkotas, K. D. 2012, Phys. Rev. D, 85, 024030