# Instability of Magnetic Equilibria in Barotropic Stars

###### Abstract

In stably stratified stars, numerical magneto-hydrodynamics simulations have shown that arbitrary initial magnetic fields evolve into stable equilibrium configurations, usually containing nearly axisymmetric, linked poloidal and toroidal fields that stabilize each other. In this work, we test the hypothesis that stable stratification is a requirement for the existence of such stable equilibria. For this purpose, we follow numerically the evolution of magnetic fields in barotropic (and thus neutrally stable) stars, starting from two different types of initial conditions, namely random disordered magnetic fields, as well as linked poloidal-toroidal configurations resembling the previously found equilibria. With many trials, we always find a decay of the magnetic field over a few Alfvén times, never a stable equilibrium. This strongly suggests that there are no stable equilibria in barotropic stars, thus clearly invalidating the assumption of barotropic equations of state often imposed on the search of magnetic equilibria. It also supports the hypothesis that, as dissipative processes erode the stable stratification, they might destabilize previously stable magnetic field configurations, leading to their decay.

###### keywords:

(magnetohydrodynamics)MHD-stars; stars: magnetic fields; stars: neutron; stars: white dwarfs; stars^{†}

^{†}pagerange: Instability of Magnetic Equilibria in Barotropic Stars–LABEL:lastpage

^{†}

^{†}pubyear: 2014

## 1 Introduction

There are several kinds of stars, namely upper main sequence stars, white dwarfs, and neutron stars, that, contrary to the Sun, have magnetic fields that are organized on large scales and persist unchanged over long time scales. Much of their interior is stably stratified, due to (a) gradients of entropy in white dwarfs and radiative envelopes of main sequence stars; and (b) gradients of composition (relative abundances of particle species) in the case of neutron stars. Clearly, no dynamo action is taking place in these stars, so their magnetic fields must be in stable hydromagnetic equilibrium states in which the Lorentz force is balanced by small perturbations to the otherwise spherical pressure and gravitational forces (with solid shear stresses in neutron star crusts possibly also playing a role).

It has long been known that purely toroidal (azimuthal) and purely poloidal (meridional) magnetic field configurations are unstable (Tayler, 1973; Markey & Tayler, 1973; Wright, 1973). However, 3-dimensional magneto-hydrodynamics (MHD) simulations of self-gravitating balls of highly conducting gas (Braithwaite & Spruit, 2004; Braithwaite & Nordlund, 2006) have shown their magnetic field to evolve naturally from an initially random configuration to a nearly axisymmetric structure in which poloidal and toroidal components of comparable amplitudes stabilize each other, persisting for many Alfvén times. Further work has suggested that the stable stratification of the fluid plays an important role in stabilizing these structures (Braithwaite, 2009; Akgün et al., 2013), and Reisenegger (2009) has conjectured that there are no stable magnetic configurations in stars that are barotropic, since they are not stably stratified.

On the other hand, attempting to construct plausible axisymmetric magnetic equilibria for these stars, several authors (Tomimura & Eriguchi, 2005; Yoshida et al., 2006; Haskell et al., 2008; Akgün & Wasserman, 2008; Kiuchi & Kotake, 2008; Ciolfi et al., 2009; Lander & Jones, 2009; Duez & Mathis, 2010; Pili et al., 2014) have imposed a barotropic equation of state, forcing the magnetic field structure to satisfy the non-linear and thus highly non-trivial Grad-Shafranov equation (Grad & Rubin, 1958; Shafranov, 1966). This assumption is not only very restrictive and unjustified for any of the stars likely to contain hydromagnetic equilibria (only very massive, radiation-dominated stars or extremely cold white dwarfs might come close to being barotropic); it might even be incompatible with having stable equilibria if the conjecture mentioned above is correct.

On the other hand, stable stratification will in the long term be eroded by non-ideal MHD processes such as heat diffusion (in main sequence and white dwarf stars), beta decays, and ambipolar diffusion (the latter two in neutron stars; see Reisenegger 2009). If the above conjecture is true, these processes would destabilize the magnetic field on these long time scales, making it decay (unless stabilized, e.g., by the solid neutron star crust).

Thus, it is important to clarify whether stable magnetic equilibria can exist in barotropic stars. A previous study by Lander & Jones (2012), based on a perturbative approach, gives a negative answer. The present work explores the same question more extensively, using full MHD simulations, in which we have evolved initially disordered magnetic fields as well as the ordered, axially symmetric, twisted-torus magnetic equilibria written analytically by Akgün et al. (2013), in both barotropic and (for comparison) stably stratified stars, in order to see whether stable equilibria can be reached.

Section 2 gives a brief introduction of the effect that stratification or its absence is expected to have on the stability of magnetic equilibria. Section 3 explains the numerical scheme used for our simulations. In Section 4, we explain how the disordered fields are created, and show the results of their evolution. In Section 5, we present the axially symmetric equilibria we use, as well as their evolution. The conclusions from our results can be found in Section 6.

## 2 Stable stratification vs. barotropy.

The physical effect of stable stratification can be illustrated by the thought experiment in which a fluid element is taken from a given position and moved vertically against gravity to a new position. Then, it is allowed to achieve pressure equilibrium with its new surroundings, without allowing it to change its entropy or chemical composition (which generally occur on much longer time scales). If the entropy or chemical composition vary with depth in the star, they will be different inside and outside the displaced fluid element, causing its density to be different, and thus creating a restoring force that is quantified by the squared Brunt-Väisälä (or buoyancy) frequency (force per unit mass per unit displacement)

(1) |

Here, is the local acceleration of gravity, is the fluid mass density, is its pressure, and the subscripts “eq” and “ad” refer to the equilibrium profile existing in the star and to the changes produced in the adiabatic displacement of the fluid element, respectively. If , the restoring force tends to move the fluid element back to its initial position, causing the fluid to be stably stratified, whereas causes a runaway fluid element, resulting in a convective instability. If, on the other hand, the fluid is barotropic, there is a one-to-one relation between pressure and density, (i. e. because the specific entropy and/or composition are uniform throughout the fluid or adjust to an equilibrium as fast as the fluid element can move), the two partial derivatives will be equal, and the fluid will be neutrally stable.

In a hydromagnetic equilibrium state, the forces on a given fluid element must be balanced,

(2) |

where is the effective gravitational potential (including a centrifugal contribution in a rotating star), is the electric current density, and is the magnetic field. From this, it follows that

(3) |

In the barotropic case, and must be parallel everywhere, so the left-hand side vanishes, forcing the right-hand side to vanish as well, and thus imposing a strong constraint (three scalar equations) on the magnetic field configuration. This constraint is much weaker for the stably stratified case. In order to understand the latter, we make the astrophysically well justified assumption that the last term in Eq. (2) is much smaller than the other two, so the thermodynamic variables can be written as and , where and are their values in an unmagnetized star, and and are small perturbations caused by the Lorentz force, which can be regarded as independent because of the third, relevant thermodynamic variable entering their relation (small perturbations of specific entropy or chemical composition; see also Reisenegger 2009; Mastrano et al. 2011). In this case, , so the left-hand side of Eq. (3) is generally non-zero, only being constrained to be a horizontal vector field (perpendicular to ). Thus, only the vertical component of the curl on the right-hand side is required to vanish (a single scalar equation), a much weaker constraint. The particular, axially symmetric case is discussed in Section 5.

In realistic stars, the characteristic Alfvén frequency associated with magnetically induced motions (e. g., the growth rate of magnetic instabilities) is much smaller than , thus the buoyancy will produce a powerful restoring force, largely impeding magnetically induced vertical (radial) motions, and thus strongly constraining possible magnetic instabilities. If, on the other hand, we had (or ), there would be no such constraints on the motion, and the magnetic field could decay much more easily. This general idea, as well as more detailed arguments based on the same physics (Braithwaite, 2009; Akgün et al., 2013) motivate the conjecture that there will be stable magnetic equilibria only in stably stratified stars (Reisenegger, 2009).

## 3 The Models

In this section we describe the setup of the model, starting with a description of the numerical code.

### 3.1 The numerical scheme.

We use a three-dimensional Cartesian grid-based MHD code developed by Nordlund & Galsgaard (1995) (see also Gudiksen & Nordlund 2005), which has been used in many astrophysical contexts such as star formation, stellar convection, sunspots, and the intergalactic medium (e.g. Padoan et al. (2007); Braithwaite (2012); Collet et al. (2011); Padoan & Nordlund (1999)). It has a staggered mesh, so that different variables are defined at different locations in each grid box, improving the conservation properties. The third-order predictor-corrector time-stepping procedure of Hyman (1979) is used, and interpolations and spatial derivatives are calculated to fifth and sixth order respectively. The high order of the discretization is a bit more expensive per grid point and time step, but the code can be run with fewer grid points than lower order schemes, for the same accuracy. Because of the steep dependence of computing cost on grid spacing (4th power for explicit 3D) this results in greater computing economy. We model the star as a “star in a box”, by modeling a self-gravitating ball of gas inside of the cubic computational domain. The simulations described here, unless stated otherwise, are run at a resolution of , and the star has a radius of grid spacings (see below).

For stability, high-order “hyper-diffusive” terms are employed for thermal, kinetic, and magnetic diffusion. These are an effective way of smoothing structure on small, badly-resolved scales close to the spatial Nyquist frequency whilst preserving structure on larger length scales. This results in a lower “effective” diffusivity compared to standard diffusion, and is more computationally efficient than achieving the same by increasing the resolution. In this study we are interested in physical processes taking place on a dynamic timescale, and so the diffusion coefficients are simply set to the lowest value needed to reliably prevent unpleasant numerical effects (‘zig-zags’), so that the diffusion timescale is as long as possible. All three diffusivities are equal, i.e. the Prandtl and magnetic Prandtl numbers are both set to . The code contains a Poisson solver to calculate the self-gravity.

### 3.2 Timescales, numerical acceleration scheme

The presence of several very different timescales in stellar MHD problems presents computational difficulties. There is the sound crossing time (governing deviations from pressure equilibrium), the Alfvén crossing time (on which a magnetic field evolves towards equilibrium), the time scale on which the magnetic field evolves under magnetic diffusion, and the rotation period of the star. We simplify the problem by assuming a non-rotating star, and by making sure the magnetic diffusion inherent in the numerical method is small enough that it does not significantly affect the evolution of the field, while still allowing for magnetic reconnection. The timescales to be included are then the sound travel time and the Alfvén time. For realistic Ap star field strengths of about a few kG, is of order a few years, which is 5 orders of magnitude longer than the sound crossing time. To cover a range like this, we make use of the fact that in a star close to hydrostatic equilibrium the evolution of the magnetic field of a given configuration depends on the field amplitude only through a factor in time scale. That is to say, if a field has the initial state , a field with initial amplitude, , where is a constant, evolves approximately as

(4) |

In other words, the time axis scales in proportion to the Alfvén crossing time. This approximation is valid as long as the Alfvén crossing time is sufficiently long compared with the sound crossing time that the evolution takes place close to hydrostatic equilibrium, and at the same time sufficiently short compared to the magnetic diffusion time scale.

We make use of this scaling to speed up the computation. The field strength of the configuration is multiplied by a time dependent factor , chosen so as to keep the Alfvén crossing time approximately constant between time steps. With Eq. (4) the resulting (unphysical) field is related to the physical field by , where is related to by . For the numerical implementation see the Appendix. This acceleration scheme (Braithwaite & Nordlund 2006) makes it possible to follow the decay of an unstable configuration to very low field strengths, by maintaining an artificially high amplitude magnetic field, and consequently shorter Alfvén time, considerably decreasing computational needs and making such studies feasible.

Since the Alfvén speed is not uniform within the star, it is necessary to decide on some suitable average for the definition of ; we use

(5) |

where , , , , and are the radius, mass, average density, average magnetic field, and total magnetic energy of the star. Similarly, the sound crossing time is defined as

(6) |

where is the total thermal energy content of the star and is the adiabatic index, further discussed in Section 3.3.

The initial field strength of the configuration , or equivalently the value of the initial Alfvén time , is an adjustable parameter. Its choice involves a compromise. A lower value of the initial field strength increases the separation between and . The approximation made is then better, but computationally more expensive. The value used here corresponds to .

### 3.3 Implementing stably stratified vs. barotropic stellar models

The code assumes a chemically uniform, monatomic, classical ideal gas, whose specific entropy is , and

(7) |

This partial derivative as given corresponds to an adiabatic change as discussed in § 2.

Simulations so far (e. g., Braithwaite & Spruit 2004; Braithwaite & Nordlund 2006; Braithwaite 2009) have taken an initial density profile inside the star corresponding to a polytrope , with

(8) |

where (here ) is the usual “polytropic index”, and the value was chosen to roughly match the radiative zones of Ap stars. For these values, since , we have and , so the star is stably stratified. We use this model as a reference against which to compare the new barotropic model for which we choose an initial profile with (), so now , , and .

The profiles of the specific entropies of both models can be seen in Figure 1(a). The star in either model is surrounded by a low-density, poorly conducting atmosphere, which causes the magnetic field outside of the star to relax to a potential field. Figure 1(b) shows the magnetic diffusivity profile for the stably stratified model. The atmosphere has a high temperature, so that its density does not become too small towards the edges of the computational domain.

Note, however, that the numerical code contains heat diffusion. Since the chosen profiles correspond to a radially decreasing temperature, the diffusion will reduce the temperature gradient, increasing the specific entropy gradient and thus making the fluid increasingly stably stratified. In order to counteract this effect and keep the fluid barotropic, we introduce an ad-hoc term in the evolution equation for the internal energy per unit volume , namely,

(9) |

where is the local temperature, is the initial value of the specific entropy inside the star, and is a timescale (to be chosen) at which this “entropy term” forces the star back to its initially isentropic structure.

## 4 Disordered initial fields

The setup of the disordered initial field is done in the same way as in Braithwaite & Spruit (2004) and Braithwaite & Nordlund (2006). This consists in assigning random phases to locations in three-dimensional wavenumber space to wavenumbers from that which corresponds to the size of the star up to wavenumbers corresponding to about four grid spacings. The amplitudes are scaled as a function of wavenumber as :

(10) |

where is the amplitude for wavenumber , and and are two randomly created phases. A reverse Fourier transformation is performed to obtain a scalar field. Two more scalar fields are produced in the same way and these become the three components of a vector potential. The vector potential is then multiplied by , where was set equal to roughly a quarter of the radius of the star. This was done to concentrate the field in the inner quarter of the star, so that the field dies off sufficiently quickly outside the star.

From this potential field, the magnetic field is then computed by taking the curl, and then its magnitude is scaled so as to obtain the desired total magnetic energy, which was times the thermal energy.

We used four different barotropic models, each including different initially disordered field configurations and an entropy timescale value roughly equal to the sound crossing time. We also include one model where a disordered magnetic field was put into a stably stratified model as a comparison. The works of Braithwaite & Spruit (2004) and Braithwaite & Nordlund (2006) have shown that disordered fields in a stably stratified model can reach stable equilibria. The evolution of the total magnetic energy of all of these models can be found in Figure 2. The conclusion is that the stably stratified model reaches a stable equilibrium in a few Alfvén timescales after which, the field decays only on the long diffusive timescale, which in a main-sequence star is of the order yr. On the other hand, none of the barotropic models reach such a stable equilibrium. Figure 3 contains 3-D renderings of the magnetic field for both a stably stratified model and a barotropic model after a time of 70. In the stably stratified model, there is a visible twisted torus structure, while the barotropic model still has a disordered magnetic field configuration. This suggests that stable stratification is a crucial ingredient for reaching a stable MHD equilibrium from a disordered field.

We have also investigated the effect of using different values of on the evolution of a disordered initial field. Simulations were carried out with four different values of , as well as a simulation where the entropy term was not used, which is akin to being equal to infinity. Figure 4 shows the evolution of the total magnetic energy for these models, as well as the stably stratified one. It is obvious that increasing the value of has little effect on the evolution, although after a few , models with a larger entropy timescale decay slightly more slowly. The reason for this is that the longer timescale allows for a slight positive entropy gradient to evolve, thus a small buoyant restoring force will act to slow down the rise of the magnetic field. It should be noted, however, that even in the case where the entropy term is not utilized in the code, the formation of a stable stratification does not occur quickly enough for a stable equilibrium to be created.

## 5 Axisymmetric equilibria

Akgün et al. (2013) wrote down analytic expressions of linked poloidal-toroidal magnetic field configurations that correspond to hydromagnetic equilibria in stably stratified stars. (A particular form of these was also used by Mastrano et al. (2011)) They have not been shown to be stable, and they are not equilibria in the barotropic case, but they might be stable equilibria (or close to these) in the stably stratified case and a good approximation to the best candidates for stability in the barotropic case.

In order to motivate their functional form, consider a general, axisymmetric field written as a sum of poloidal and toroidal components,

(11) |

where and are characteristic values of the two components, and and are dimensionless scalar functions that depend on , the radial coordinate, and , the polar angle in the star. Since there cannot be an azimuthal component of the Lorentz force (which, in axial symmetry, cannot be balanced by pressure gradients and gravity), and must be parallel, thus one can write .

The latter condition (a special case of the condition of the vanishing vertical component of discussed in Section 2) is the only condition required to be satisfied by and in order to yield a hydromagnetic equilibrium in a stably stratified star. Besides this, these functions are arbitrary and need not satisfy, e. g., the commonly assumed Grad-Shafranov equation, which is physically justified only for barotropic fluids and otherwise represents an arbitrary, additional constraint (e. g., Reisenegger 2009). Of course, in neither case there is a guarantee for the stability, and thus for the astrophysical relevance, of the constructed equilibria.

For the poloidal part, Akgün et al. (2013) choose

(12) |

where x is the dimensionless radial coordinate, , for the stellar radius, and

(13) |

where the are constants. Outside the star, this gives an exact dipole. For the internal field, three of the four constants can be obtained from the condition that all components of the magnetic field are continuous across the stellar boundary (implying and , where the primes denote derivatives with respect to ) and that the current density vanishes at the surface (implying ). The fourth constant, here chosen as , is a free parameter, which can be tuned to change the size of the closed field line region of the equilibria, as seen in Figure 5.

The closed field line region is important because the toroidal component of the field, described by the function , can only exist in the volume where the poloidal field lines are closed within the star. This motivates the choice:

(14) |

where the exponent , so that the current due to the toroidal field decreases smoothly to zero at the stellar surface.

### 5.1 Stability tests in stably stratified models.

We do not know a priori whether these equilibria will be stable in any type of star. As a first step, we verified their stability in a stably stratified model. If the magnetic field configurations are not stable in the stably stratified model, it seems likely that they will not be stable in a barotropic model.

We utilize the aforementioned equilibria of Akgün et al. (2013) with the toroidal exponent of Eq. 14 set equal to 2, as was done in their analytic work, and from Eq. 13 equal to -100. The reason for choosing this value of was two-fold. The equilibria found to evolve from random field configurations in the stably stratified models were found to have larger tori than the equilibria of Akgün et al. (2013) in which the value of was set to zero. Secondly, we chose to be a large amplitude, so as to be spread out the toroidal flux across a larger area of the star. This would have the effect of decreasing the local magnetic energy density in the torus, as compared to a configuration with a smaller closed field line region with the same fraction of the toroidal to total magnetic energy. In spreading out the toroidal energy, the local value of , defined as , in the torus would be larger, and it would thus be less buoyant.

We ran a number of simulations with various initial poloidal field energies relative to total magnetic energy , ranging between 0.0015 and 0.93. Such a large range was used, because at low values of this parameter, the Tayler instability (Tayler, 1973) is expected to set in, and at high values, instability along the “neutral line” may occur (Markey & Tayler, 1973; Wright, 1973). Models with initial between 0.008 and 0.8 were found to be stable, and to decay on a diffusive timescale, while all others decayed more quickly due to instabilities. Figure 6 shows the evolution of the magnetic energy for some of the models. To try to diagnose the source of instabilities, we analyzed the evolution of the =0-3 azimuthal modes, which were determined by performing a root-mean-square integration of the component of the velocity in the meridional plane as:

(15) | ||||

where is the amplitude of the theta component of the velocity for a particular -mode, is the area of the star in the meridional plane, and is the component of the velocity. We then calculated the component of the kinetic energy in each of these modes, and plotted these values normalized to the initial total magnetic energy versus time in Figure 7. Here, it is evident that the model with of 0.5 is stable to all of the modes, as the curves for each mode evolve to a steady value. The model with initial of 0.87 was unstable to m=2 which can be seen by the sharp peak in the m=2 curve at a time of roughly 2-5 after which this model calms to a new equilibrium as the m=2 mode flattens out. The model with initial of 0.0025 was unstable to the m=1 mode which begins to set in at a time of roughly 4, as can be seen by the peak in the figure. It was then found that all configurations with initial above 0.8 were unstable to the m=2 mode, while those with initial below 0.008 were unstable to the m=1 mode. In the m=1 instability, the toroidal component of the field can be thought of as a stack of rings on top of one another; as the instability occurs, the rings slip with respect to each other. In the m=2 instability, the region along the neutral line, which is the space where the poloidal field goes to zero, experiences tension, which results in kinking of the neutral line, stretching any toroidal field that may be present along the neutral line. This kinking of the neutral line bends the toroidal field lines into a shape that is similar to the seam of a tennis ball. Both the =1 and =2 unstable configurations end up reaching another stable equilibrium after some time. In the case of the =1 unstable model, the equilibrium is reached at a much later time. However, it can be seen that the =2 unstable model reaches a new “tennis ball” shaped equilibria after roughly 8 when it begins to decay on a diffusive timescale, and the kinetic energy amplitude in the =2 mode drops off and becomes stable in Figure 7. These results are similar to those of Braithwaite (2009), who studied the stability of twisted torus configurations in stably stratified stars.

### 5.2 Stability tests in barotropic models

With the knowledge that these equilibria are stable in a stably stratified star for initial values between about 0.008 and 0.8, we now investigate whether such stable equilibria can exist in barotropic stars.

#### 5.2.1 Existence of stable equilibria?

A series of models have been run with the same equilibria used in barotropic models. For these models, was set to be equal to the sound crossing time of the star. The initial values used were between 0.03 and 0.96. The motivation for using values of above 0.9 comes from the works by Tomimura & Eriguchi (2005); Ciolfi et al. (2009); Lander & Jones (2009); Armaza et al. (2014), where magnetic equilibria have been calculated in axial symmetry for barotropic stars, and none of which obtained 0.9. The rest of the spectrum of values was chosen based on the results of Section 5.1, where for a stably stratified star, the equilibria were stable for initial values of 0.008 to roughly 0.8. In addition, Ciolfi & Rezzolla (2013) have calculated magnetic equilibria in barotropic stars and found equilibria with values as low as roughly 0.11.

Figure 8 shows the evolution of the total magnetic energy for barotropic models, as well as a stably stratified model with 0.5 for comparison. It is evident from the figure that none of the equilibria are stable in the barotropic models. An interesting point, however, is that fields with a initial values between 0.33 and 0.2 decay much more slowly at earlier times than all other barotropic models. The reason for the slower decay of the magnetic field for these models can be seen in Figure 9, where the effective magnetic radius, , defined as:

(16) |

is plotted. It is evident that for magnetic configurations with an initial 0.33, increases, caused by the torus rising in the star. For values below 0.33, the magnetic tension in the torus results in its contraction, causing a decrease of . However, these models that initially experience a contraction of the torus still show that the value will increase after roughly 6 , as the torus rises out of the star.

#### 5.2.2 Effect of

We also investigated the effect of varying for a magnetic configuration with an initial value of 0.5 in Figure 10. It is evident that a longer allows the magnetic field to decay more slowly. As the entropy timescale becomes larger, the star, due to numerical diffusion, will evolve an entropy gradient, which will impede the rise of the torus. However, even without the entropy adjustment term, the equilibria are still not stable.

#### 5.2.3 Variations of the free parameters in the axially symmetric equilibria.

It was shown in Section 5.2.1 that, regardless of the initial fraction used, all the magnetic configurations were unstable. However, the initial fraction of the total magnetic energy in the poloidal component is just one of the free parameters that can be investigated in the axisymmetric configurations. The effect of varying other free parameters, namely the value of in Eq. 14, which affects the distribution of the toroidal component, the value of in Eq. 13, which controls the size of the closed field line region, and the radius of the “neutral line” were also studied. We found that variations to these parameters had no effect on the final state of the simulations, as all configurations were found to be unstable to the same effect seen in Section 5.2.1, where the torus flowed radially out of the star, subsequently decaying in the atmosphere.

### 5.3 Instabilities

We have analyzed the instability of the -modes in the barotropic models, with the same method as the stably stratified models (See Section 5.1). Figure 11 shows the evolution of the kinetic energies for each of the m=0-3 modes for models with four different initial values. In all cases, the m=0 mode is the dominant instability, quickly rising to a value of a few hundredths of a percent of the total initial magnetic energy, at which point it remains at this range. Even in the cases where and , which in the stably stratified star were found to be unstable to the m=2 and m=1 modes respectively, the m=0 mode is the dominant instability in the barotropic models. To help visualize what the m=0 mode looks like, Figure 12 shows that the cause of the m=0 mode is the rise of the torus, driven by a combination of its own buoyancy and the pressure of the enclosed poloidal field. Thus, the stable stratification seems to play a very strong role in suppressing the m=0 mode, as it always dominates in the absence of stable stratification.

In order to see how robustly this mode dominates in de-stabilizing the magnetic equilibria, we have investigated what happens when an initial non-axially symmetric perturbation is added to the system. To do so, we introduced perturbations to the component of the velocity along the equatorial plane within the star. One case included an initial =1 perturbation for a model with initial =0.03, a second model contained an =2 perturbation with initial =0.74, while the third model contained perturbations to =1, 2, and 3 modes for an initial =0.5. The strength of all perturbations was the same, with each mode containing an initial =.01%, roughly of the same order that the =0 mode was seen to reach in the non-perturbed models. We again found that regardless of the initial perturbations added, the =0 mode always became the dominant mode of instability.

In addition to non-axially symmetric perturbations, we have run two simulations with a non-equatorially symmetric perturbation. The perturbation was created by offsetting the magnetic field configuration by 0.016 in the Northern direction, for two initial values of 0.5 and 0.24. Snapshots of the time evolution of the poloidal and toroidal field lines can be seen in Figure 13. The evolution in these northerly perturbed models is similar to their non-perturbed counterparts, with the only difference being that these models have a slight tendency to evolve away from the equator.

## 6 Conclusions

We have conducted two different kinds of numerical experiments.

In the first kind, we started with disordered initial fields. In the previously studied case of stably stratified stellar models, we confirmed that the magnetic field evolved into an ordered, stable, roughly axially symmetric configuration with comparable poloidal and toroidal fields. In the barotropic cases we studied, this never happened, and the magnetic energy decayed much faster than expected from diffusion.

In the second kind, we started with a smooth, axially symmetric magnetic field with poloidal and toroidal components. We confirmed that, in the stably stratified star, there are stable hydromagnetic equilibria for a fairly wide range of values of the initial fraction of poloidal to total magnetic energy, . Outside this range, the field decayed through non-axisymmetric modes, in the toroidally dominated and in the poloidally dominated case. In contrast, in the barotropic case, we found no stable configurations (exploring a fairly large parameter space), and the field always decayed through an axially symmetric () instability in which the toroidal flux rose radially out of the star, being dissipated in the atmosphere.

Although far from a rigorous mathematical proof, this provides strong evidence (added to that previously provided by Lander & Jones 2012) that there are no stable equilibria in barotropic stars. It also strongly supports the intuition that, in a stably stratified star, the buoyancy force strongly constrains the potential instabilities by impeding any substantial radial motions, whereas such motions are not hindered in the barotropic case, and in fact these motions destabilize the magnetic equilibrium.

As argued above, the stabilization provided to the magnetic field by the buoyancy force can be roughly quantified by comparing the Brunt-Väisälä frequency to the Alfvén frequency . It will only be effective if , whereas in the opposite case the star would behave as if it were barotropic, not being able to contain a stable magnetic field. If there are dissipative processes effectively eroding the stable stratification on long time scales, they will lead to a destabilization and eventual decay of the magnetic field, unless it can be stabilized, e. g., by the solid crust of a neutron star. In fact, in the neutron star case, this effect has been argued to act on astrophysically relevant time scales, which does not seem to be the case in white dwarfs and upper main sequence stars (Reisenegger, 2009), with the possible exception of very massive O stars, which are radiation-dominated and thus only weakly stratified.

Note that the numerical models used here, which assume a chemically uniform, classical ideal gas, stably stratified by an entropy gradient, clearly do not directly apply to degenerate stars, such as white dwarfs and neutron stars. However, we have no reason to doubt that the essential physics, in particular the competition between magnetic and buoyancy forces, is the same in these cases, and the (very rough) condition is still required to stabilize a hydromagnetic equilibrium in these cases, even in the neutron star case, where is due to a composition gradient. Additional changes in the neutron star case are the presence of a stabilizing solid crust, as well as superconducting and superfluid regions in the interior, which modify the form of the Lorentz force and the dynamical equations and thus do not allow a direct application of the present results, although some of their features might carry over to this regime.

## 7 Acknowledgements

This work was supported by the DFG-CONICYT International Collaboration Grant DFG-06, FONDECYT Regular Project 1110213, and Proyecto de Financiamiento Basal PFB-06/2007. Some of the figures were made with VAPOR (www.vapor.ucar.edu).

## Appendix A Numerical acceleration scheme

The scaling of magnetic fields that evolve from different initial amplitudes (but otherwise identical configuration) with time is described in Section 3.2. The practical implementation in the MHD code is as follows. We need to make a distinction between the field and time in the accelerated code, and the physical field and time reconstructed from it. The MHD induction equation

(17) |

is first evolved over a time step , where , to yield an intermediate update . The magnetic energy over the numerical volume is measured from the result of this time step:

(18) |

and the evolution time scale (in the numerical time unit) is calculated:

(19) |

The intermediate update is modified:

(20) |

which brings the numerical Alfvén crossing time back to its initial value . The new value of the physical magnetic energy has changed by the amount

(21) |

and the physical time coordinate by the amount

(22) |

## References

- Akgün et al. (2013) Akgün T., Reisenegger A., Mastrano A., Marchant P., 2013, MNRAS, 433, 2445
- Akgün & Wasserman (2008) Akgün T., Wasserman I., 2008, MNRAS, 383, 1551
- Armaza et al. (2014) Armaza C., Reisenegger A., Valdivia J. A., Marchant P., 2014, in IAU Symposium Vol. 302. pp 419–422
- Braithwaite (2009) Braithwaite J., 2009, MNRAS, 397, 763
- Braithwaite (2012) Braithwaite J., 2012, MNRAS, 422, 619
- Braithwaite & Nordlund (2006) Braithwaite J., Nordlund Å., 2006, A&A, 450, 1077
- Braithwaite & Spruit (2004) Braithwaite J., Spruit H. C., 2004, Nature, 431, 819
- Ciolfi et al. (2009) Ciolfi R., Ferrari V., Gualtieri L., Pons J. A., 2009, MNRAS, 397, 913
- Ciolfi & Rezzolla (2013) Ciolfi R., Rezzolla L., 2013, MNRAS, 435, L43
- Collet et al. (2011) Collet R., Hayek W., Asplund M., Nordlund Å., Trampedach R., Gudiksen B., 2011, A&A, 528, A32
- Duez & Mathis (2010) Duez V., Mathis S., 2010, A&A, 517, A58
- Grad & Rubin (1958) Grad H., Rubin H., 1958, Int. Atom. Ener. Ag. Conf. Proc., 31, 190
- Haskell et al. (2008) Haskell B., Samuelsson L., Glampedakis K., Andersson N., 2008, MNRAS, 385, 531
- Hyman (1979) Hyman J. M., 1979, in Advances in Computer Methods for Partial Differential Equations - III A method of lines approach to the numerical solution of conservation laws. pp 313–321
- Kiuchi & Kotake (2008) Kiuchi K., Kotake K., 2008, MNRAS, 385, 1327
- Lander & Jones (2009) Lander S. K., Jones D. I., 2009, MNRAS, 395, 2162
- Lander & Jones (2012) Lander S. K., Jones D. I., 2012, MNRAS, 424, 482
- Markey & Tayler (1973) Markey P., Tayler R. J., 1973, MNRAS, 163, 77
- Mastrano et al. (2011) Mastrano A., Melatos A., Reisenegger A., Akgün T., 2011, MNRAS, 417, 2288
- Nordlund & Galsgaard (1995) Nordlund Å., Galsgaard K., 1995, http://www.astro.ku.dk/~aake/papers/95.ps.gz
- Padoan & Nordlund (1999) Padoan P., Nordlund Å., 1999, ApJ, 526, 279
- Padoan et al. (2007) Padoan P., Nordlund Å., Kritsuk A. G., Norman M. L., Li P. S., 2007, ApJ, 661, 972
- Pili et al. (2014) Pili A. G., Bucciantini N., Del Zanna L., 2014, MNRAS, 439, 3541
- Reisenegger (2009) Reisenegger A., 2009, A&A, 499, 557
- Shafranov (1966) Shafranov V. D., 1966, Reviews of Plasma Physics, 2, 103
- Tayler (1973) Tayler R. J., 1973, MNRAS, 161, 365
- Tomimura & Eriguchi (2005) Tomimura Y., Eriguchi Y., 2005, MNRAS, 359, 1117
- Wright (1973) Wright G. A. E., 1973, MNRAS, 162, 339
- Yoshida et al. (2006) Yoshida S., Yoshida S., Eriguchi Y., 2006, ApJ, 651, 462