# Stability of axially symmetric magnetic fields in stars

Laura Becerra,1 Andreas Reisenegger,2 Juan Alejandro Valdivia,3,4 Mikhail
Gusakov5
1Instituto de Astrofísica, Facultad de Física, Pontificia Universidad Católica
de Chile, Av. Vicuña Mackenna 4860, Macul, Santiago, Chile
2Departamento de Física, Facultad de Ciencias Básicas, Universidad
Metropolitana de Ciencias de la Educación, Av. José Pedro Alessandri 774,
Ñuñoa, Santiago, Chile
3Departamento de Física, Facultad de Ciencias, Universidad de Chile, Las
Palmeras 3425, Ñuñoa, Santiago, Chile
4Centro para el Desarrollo de la Nanociencia y Nanotecnología, CEDENNA,
Santiago, Chile
5Ioffe Institute, 26 Politekhnicheskaya Street, St. Petersburg 194021, Russia
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The magnetic fields observed in Ap-stars, white dwarfs, and neutron stars are
known to be stable for long times. However, the physical conditions inside the
stellar interiors that allow these states are still a matter of research. It
has been formally demonstrated that both purely toroidal and purely poloidal
magnetic fields develop instabilities at some point in the star. On the other
hand, numerical simulations have proved the stability of roughly axisymmetric
magnetic field configurations inside stably stratified stars. These
configurations consist of mutually stabilizing toroidal and poloidal
components in a twisted torus shape. Previous studies have proposed rough
upper and lower bounds on the ratio of the magnetic energy in the toroidal and
poloidal components of the magnetic field. With the purpose of mapping out the
parameter space under which such configurations remain stable, we used the
Pencil Code to perform 3D magnetohydrodynamic simulations of the evolution of
the magnetic field in non-rotating, non-degenerate stars in which viscosity is
the only dissipation mechanism, both for stars with a uniform (barotropic) and
radially increasing (stably stratified) specific entropy. Furthermore, we
considered different conditions regarding the degree of stable stratification
and the magnetic energy in each component, roughly confirming the previously
suggested stability boundaries for the magnetic field.

###### keywords:

MHD – stars: magnetic field – stars: massive – stars: neutron – stars: white
dwarfs – software: simulations

††pubyear: 2022††pagerange: Stability of axially symmetric magnetic fields in
stars–A

##  1 Introduction

The magnetic fields observed in stars can be grouped into two general types
based on their phenomenology. The first group is formed by the magnetic fields
observed in low-mass non-degenerate stars. These have a complex structure, are
variable on time scales that go from weeks to decades, and are probably
produced by a contemporary dynamo in the star’s convective envelopes (Reiners,
2009). In the second group are the magnetic fields detected in some
intermediate and high-mass main-sequence stars, white dwarfs, and neutron
stars, which are long-lived and globally organized. The partly (or completely)
stably stratified interiors of these stars are not compatible with a dynamo
origin for their magnetic field. These fields are slowly varying or stable,
likely having been set up in an earlier stage of the evolution of these stars
(Donati & Landstreet, 2009; Sikora et al., 2019).

It has yet to be fully understood what conditions need to be satisfied by the
latter group of stars and their magnetic field configuration so that they are
able to survive over the star’s lifetime. Since Chandrasekhar & Fermi (1953),
many works have been devoted to building magnetic field configurations by
solving the hydromagnetic equilibrium equation numerically or analytically and
to prove the stability of these configurations over dynamical timescales.
Tayler (1973) found general conditions for the stability of purely toroidal
fields111Within the Cowling approximation (that is, neglecting perturbations
of the gravitational potential; Cowling 1941), these conditions are necessary
and sufficient for stability. Relaxing this approximation, the conditions
remain necessary, but not sufficient. and used them to show that any such
field having a finite current density along the axis will be unstable in this
region. The application of these conditions to more general field
configurations (Goossens et al., 1981; Akgün et al., 2013) proved that all
toroidal fields are unstable. Similarly, purely poloidal magnetic fields were
also found to be subject to dynamical instabilities (Markey & Tayler, 1973;
Wright, 1973; Flowers & Ruderman, 1977), suggesting that a combination of
poloidal and toroidal components is required for the stability of axially
symmetric magnetic fields (Prendergast, 1956).

Later, the numerical simulations of Braithwaite & Nordlund (2006) (see also
Braithwaite & Spruit, 2004) found that roughly axisymmetric stable equilibria,
consisting of both toroidal and poloidal components of comparable strength in
a twisted-torus shape, formed spontaneously from initially random magnetic
field configurations in stably stratified stars. In Becerra et al. (2022)
(hereafter Paper I), we showed that the formation of these axially symmetric
magnetic fields is strongly influenced by the use of a high-order diffusion
scheme introduced for numerical stability purposes, particularly by the
presence of a substantial magnetic diffusivity. Based on our simulations, we
argued that, under more realistic physical conditions, the magnetic field
configuration is likely to evolve to a stable, but non-axially symmetric
state. However, due to their simplicity, the axially symmetric field
configurations are the most widely studied among the community, and
observations seem to support the existence of such equilibrium configurations
in at least some stars (Aurière et al., 2007).

In such configurations, the poloidal and toroidal magnetic field components
need to stabilize each other. For the stabilization of a poloidal field, the
toroidal field must be at least of a similar strength, as first suggested by
Markey & Tayler (1973) and later confirmed by Braithwaite (2009), whose
simulations yielded the approximate condition that the energy in the poloidal
component, $E_{\rm pol}$, can be at most $\sim 80\%$ of the total magnetic
energy $E_{\rm mag}$, for the field to be stable.

For the stability of the toroidal field, it is not enough to invoke the
presence of a poloidal component, but it is also necessary for the stellar
matter to be stably stratified, with a positive squared Brunt-Väisälä (or
buoyancy) frequency,

$N^{2}\equiv\left({1\over\gamma}-{1\over\Gamma}\right){\rho g^{2}\over p}>0,$
(1)

where $\rho$, $p$, and $g$ are the local values of the fluid density,
pressure, and acceleration of gravity, respectively. $\Gamma$ is the adiabatic
index,

$\Gamma=\left(\frac{\partial{\rm ln}\,p}{\partial{\rm ln}\,\rho}\right)_{X},$
(2)

where $X$ represents one or more conserved quantities characterizing adiabatic
perturbations, such as specific entropy in the case of main sequence stars and
white dwarfs, and relative particle abundances in the case of neutron stars
(Reisenegger, 2009), while $\gamma$ is the conventional polytropic index:

$\gamma=\frac{d\,{\rm ln}\,p}{d\,{\rm ln}\,\rho}\,,$ (3)

characterizing the equilibrium inside the star. In the astrophysical examples
mentioned above, the typical buoyancy frequency is much larger than the Alfvén
frequency, $\omega_{A}=B/(\sqrt{\mu_{0}\rho}R_{s})$, where $B$ is the magnetic
field strength, $R_{s}$ is the stellar radius and $\mu_{0}$ is the magnetic
vacuum permeability. In this case, the stable stratification strongly
suppresses radial motions. In particular, instabilities of the toroidal field
require displacements with a non-zero radial component, but this must be much
smaller than the horizontal components if the star is strongly stratified
(Akgün et al., 2013), and this helps the poloidal component to stabilize such
motions. Braithwaite (2009) used heuristic arguments to write the condition
for the stabilization of the toroidal component as $E_{\rm pol}/E_{\rm
mag}\gtrsim aE_{\rm mag}/|E_{\rm grav}|$, where $E_{\rm grav}$ is the star’s
gravitational energy and $a$ is a dimensionless constant that is inversely
proportional to the difference $\Gamma-\gamma$, and which he found numerically
to be $\sim 10$ for the case of main-sequence stars, for which $\Gamma\approx
5/3$ and $\gamma\approx 4/3$. Through a mostly analytic calculation, Akgün et
al. (2013) suggested $a\approx 1.8/(\Gamma/\gamma-1)$, in rough agreement with
Braithwaite’s result.

The two conditions for the mutual stabilization of the toroidal and poloidal
field components discussed above can be written together as222Here it is
assumed that $N\gg\omega_{A}$, which implies that the term on the right-hand
side is large, as is the case in the stars of interest and in all the
simulations used to obtain the stability limits.

$0.25\lesssim\frac{E_{\rm tor}}{E_{\rm pol}}\lesssim
0.5\sqrt{\left(\frac{\Gamma}{\gamma}-1\right)\frac{|E_{\rm grav}|}{E_{\rm
pol}}},$ (4)

showing the two essential ingredients for the stability of the magnetic field
in the stellar interior: the relative strengths of the toroidal and poloidal
components of the magnetic field, and the star’s stable stratification
($\Gamma>\gamma$). The importance of the latter effect was confirmed by
simulations that failed to find any stable configuration in barotropic (non-
stratified, $\Gamma=\gamma$) stars (Lander & Jones, 2012; Mitchell et al.,
2015).

We note that, even disregarding the issue of stability, there is an important
difference between hydromagnetic equilibria in barotropic stars (with $N=0$)
and those with strong stable stratification ($N^{2}\gg\omega_{A}^{2}$). In the
latter, pressure and density can be treated as independent variables, so the
fluid has two scalar degrees of freedom, and the only condition for an axially
symmetric magnetic field to correspond to a possible equilibrium state is that
the toroidal component of the Lorentz force vanishes, leading to the condition
$\beta=\beta(\alpha)$ for the scalar functions introduced in Equation (5). In
the former, there is a one-to-one relation between pressure and density, so
the fluid has only one scalar degree of freedom, and the functions $\alpha$
and $\beta$ are additionally required to satisfy the Grad-Shafranov equation
(Grad & Rubin, 1958; Shafranov, 1958). This is a highly nonlinear equation,
whose general solution space has not yet been characterized. Various authors
have explored numerical solutions, generally finding only configurations with
$E_{\rm tor}/E_{\rm pol}<0.1$ (see for example Yoshida & Eriguchi, 2006;
Lander & Jones, 2009; Ciolfi et al., 2009; Fujisawa et al., 2012; Armaza et
al., 2015, among others), which was explained through an approximate
analytical argument in Armaza et al. (2015). If the magnetic field
configuration is totally confined inside the star, the $E_{\rm tor}/E_{\rm
pol}$ fraction could become larger (Haskell et al., 2008; Duez & Mathis,
2010). However, Lander & Jones (2012) studied the stability of these
configurations, finding that they were all unstable in barotropic stars.

Equation (4) also has important astrophysical consequences. For example, since
$E_{\rm mag}\ll|E_{\rm grav}|$, a star can store a (internal) toroidal
component much stronger than the poloidal field measured on its surface. If
this is the case, the anisotropic pressure of the internal magnetic field
deforms the star into a prolate shape (Chandrasekhar & Fermi, 1953). This
could cause the precession (Wasserman, 2003) inferred in some pulsars, such as
PSR B1828-11 (Stairs et al., 2000; Ashton et al., 2017), which has also been
suggested to be relevant for fast radio bursts (FRBs) from magnetars
(Wasserman et al., 2022). It would also lead to the emission of continuous
gravitational waves by rotating neutron stars (Dall’Osso et al., 2009),
contributing to their spin-down (Cutler, 2002; Akgün et al., 2006). On the
other hand, a strong internal toroidal field could provide an energy source
for the violent activity of magnetars (Thompson & Duncan, 1996; Mereghetti,
2008) and the large pulsed fraction in the X-ray emission of central compact
objects (Shabaltas & Lai, 2012; Viganò & Pons, 2012).

In this paper, we address the problem of the stability of axially symmetric
magnetic fields in stably stratified stars by simulating the dynamical
evolution of many different configurations, in order to verify the accuracy of
Equation (4). For this, we use the Pencil Code (Pencil Code Collaboration et
al., 2021), a high-order finite-difference code for compressible hydrodynamic
flows with magnetic fields. For all the simulations presented in this paper,
the only dissipative mechanism considered is ordinary shear viscosity, while
all hyper-diffusion coefficients are taken to vanish. As argued in Paper I,
this approach is expected to most closely mimic the physics of the real stars.
Section 2 gives a brief description of the numerical setup and initial
conditions, and section 3 is devoted to testing the stability of different
axially symmetric magnetic field configurations as a function of the ratio of
magnetic to gravitational energy, the fraction of the magnetic energy
contained in the poloidal and toroidal components, and the stable
stratification of the stellar matter. Our conclusions are presented in section
4.

##  2 Initial setup

We test the stability of axially symmetric magnetic field configurations by
simulating their dynamical evolution inside the star with the Pencil Code333
https://github.com/pencil-code/ (Pencil Code Collaboration et al., 2021),
which uses sixth-order centered spatial derivatives and a third-order Runge-
Kutta time-stepping scheme to solve the magneto-hydrodynamics (MHD) equations
(given in Appendix A). All the simulations are performed in a cubic box of
side $L_{\rm box}=4.5R_{\rm s}$, centered on the star, with a Cartesian grid
with $128^{3}$ equally spaced points and periodic boundary conditions.

The fluid inside the simulation box obeys the ideal gas equations of state
with an adiabatic index $\Gamma=5/3$ (see Equation [2]), with specific entropy
as the stabilizing quantity, $X=s$. As done in Paper I, to prepare the initial
condition for the simulations, we first built a non-rotating, unmagnetized
spherical star model. Thus, we solved the (non-magnetic) hydrostatic
equilibrium equations, adopting, inside the star, a polytropic relation
between the gas pressure and density: $p=K\rho^{\gamma}$, where $K$ and
$\gamma$ are constants. When $\gamma<\Gamma$, the star is stably stratified
(the specific entropy is an increasing function of the radius), while for
$\gamma=\Gamma$, it is barotropic (no entropy gradient). For the initial setup
of the barotropic case ($\Gamma=\gamma=5/3$) and the slightly stratified case
($\Gamma=5/3$ and $\gamma=1.61$), it was not possible to apply the polytropic
relation to the whole stellar interior, because this made the density near the
surface very low, causing problems to the numerical scheme. Thus, for these
cases, the specific entropy was taken to increase for $r>(0.97-0.98)R_{s}$,
creating a thin, stably stratified surface layer, as seen in Figure 1 of Paper
I.

Outside the star, we placed a low-density atmosphere with a uniform
temperature and a low electrical conductivity, to make the magnetic field
quickly relax to a potential field. Inside the star, the electrical
conductivity was assumed to be very large. So, we set the magnetic diffusivity
equal to zero inside the star and large and constant outside it (for a
complete account of the initial setup, we refer the reader to Paper I, ).

For this spherically symmetric equilibrium configuration, we calculated the
gravitational potential, which was kept constant in all our simulations. We
note that by doing this, we apply the Cowling approximation (Cowling, 1945) in
two different ways. First, we neglect the perturbation of the gravitational
potential produced by the presence of the initial magnetic field, since the
magnetic equilibria of the star can be seen as a small perturbation of its
non-magnetic spherical background. And second, we also neglect the changes in
the gravitational potential along the simulations, due to the evolution of the
magnetic field (The latter is the “standard” Cowling approximation used by
Tayler (1973) and others.)

Once this was done, we introduced the desired initial magnetic field and
established the initial magnetic equilibrium by running the code for a few
sound crossing timescales while freezing the evolution of the magnetic field,
allowing the hydrodynamical forces to balance the Lorentz force. The latter
was only done for the stably stratified models, since the initial magnetic
field configurations we chose do not satisfy the Grad-Shafranov equation, so
they do not correspond to equilibria for barotropic stars (see the discussion
in the Introduction). Thus, for this case, it was not attempted to make the
fluid relax to an equilibrium before starting the simulation.

Finally, to accelerate the appearance of eventual instabilities, if not
specified otherwise, a small random perturbation was added to the density
field, and the simulation was continued, now using the full set of equations
described in Appendix A, so the magnetic field was allowed to evolve.

In all the simulations, the viscosity coefficient is $\nu=8.24\times
10^{-4}R_{s}/\tau_{s}$. We note that our results are independent of the
dissipation mechanism for the kinetic or the magnetic energy, which we
verified by also running simulations with the hyper-diffusion scheme studied
in Paper I, and obtaining the same results for the stability of each magnetic
field configuration considered. We refer the reader to Paper I for more
details of the initial setup.

###  2.1 Magnetic field configuration

(a) Field I. (b) Field II.

Figure 1: View of the two choices of axially symmetric magnetic field
configuration in the star’s meridional plane, with the symmetry axis in the
vertical direction: (a) Field I: analytic expression from Akgün et al. (2013);
(b) Field II: obtained numerically from a random initial field (Paper I; see
text for details).

In order to verify the validity of Equation (4), we follow the evolution of
ordered axially symmetric initial magnetic field configurations in the stellar
interior, under different initial conditions for the star’s stable
stratification (or its absence), total magnetic energy, ratio between the
energies in the poloidal and toroidal components of the magnetic field, and
two different functional forms for the magnetic field structure, “Field I” and
“Field II”, shown in Figure 1.

For Field I, we use the configuration built analytically in Akgün et al.
(2013). In spherical coordinates $(r,\theta,\phi)$, an axially symmetric
magnetic field can be written as the sum of a poloidal and a toroidal
component:

$\vec{B}=\vec{B}_{\rm pol}+\vec{B}_{\rm tor}=a_{\rm
pol}\vec{\nabla}\alpha(r,\theta)\times\vec{\nabla}\phi+a_{\rm
tor}\beta(r,\theta)\vec{\nabla}\phi$ (5)

(Chandrasekhar & Prendergast, 1956), where $\alpha(r,\theta)$ and
$\beta(r,\theta)$ are scalar functions further discussed below,
$\vec{\nabla}\phi=\hat{\phi}/(r\sin\theta)$, and $a_{\rm pol}$ and $a_{\rm
tor}$ are adjustable coefficients to set the strengths of the poloidal and
toroidal components. In equilibrium, the azimuthal component of the Lorentz
force must vanish, implying $\beta=\beta(\alpha)$. For the poloidal component
of the magnetic field, Akgün et al. (2013) adopted a simple form that makes
the external field a pure vacuum dipole:

$\alpha(x,\theta)=f(x)\sin^{2}\theta\,,$ (6)

with $x=r/R_{s}$ and

$f(x)=\left\\{\begin{aligned}
f_{2}x^{2}+f_{4}x^{4}+f_{6}x^{6}+f_{8}x^{8}&\quad&\mathrm{for}&\quad x\leq
1,\\\ x^{-1}&\quad&\mathrm{for}&\quad x>1.\end{aligned}\right.$ (7)

The coefficients $f_{i}$ are set by applying boundary conditions at the
stellar surface (continuity of the magnetic field and the current density),
which yield $f_{2}=\frac{35}{8}-f_{8}$, $f_{4}=-\frac{21}{4}+3f_{8}$,
$f_{6}=\frac{15}{8}-3f_{8}$, with $f_{8}$ as a free parameter. The function
$\beta$ must be chosen so that the toroidal field is non-zero only on poloidal
field lines that close within the star. In Akgün et al. (2013), it is taken to
be

$\beta(\alpha)=\left\\{\begin{aligned}
\left(\alpha-1\right)^{2}&\quad&\mathrm{for}&\quad\alpha\geq 1,\\\
0&\quad&\mathrm{for}&\quad\alpha<1.\end{aligned}\right.$ (8)

We choose the free parameter $f_{8}=-1000$, so that the torus containing the
toroidal field lines extends over a large fraction of the stellar volume.

For Field II, we use a magnetic equilibrium obtained numerically from the
evolution of a random initial configuration. Specifically, we take the
magnetic field from Model V of Paper I at $t=240~{}\tau_{A,0}$, from which we
build an axially symmetric configuration by taking the azimuthal average of
each spherical component of the magnetic field ($B_{r},B_{\theta},B_{\phi}$)
around the magnetic axis, $\vec{M}$, defined as the direction that minimizes
the asymmetry parameter (see also Paper I, ):

$\mathcal{A}=\min_{\theta,\phi}\frac{\int_{\rm star}|\vec{B}-\vec{B}^{\rm
axial}_{\theta,\phi}|^{2}\,dV}{\int_{\rm star}|\vec{B}|^{2}\,dV}.$ (9)

Here, $\vec{B}^{\rm axial}_{\theta,\phi}$ is built by taking the azimuthal
average of the magnetic field components in spherical coordinates around a
certain axis oriented in the direction given by the angles $\theta$ and
$\phi$.

We define the Alfvén speed as

$v_{A}\equiv\frac{B_{\rm rms}}{\sqrt{\mu_{0}\rho_{\rm rms}}},$ (10)

and the Alfvén and sound crossing time scales as

$\tau_{A}\equiv\frac{R_{s}}{v_{A}}\quad{\rm
and}\quad\tau_{s}\equiv\frac{R_{s}}{c_{\rm s,rms}},$ (11)

respectively, where $c_{s}$ is the sound speed, and we use $A_{\rm rms}$ to
denote the root-mean-square volume averaged over the simulation box of any
quantity $A$.

##  3 Stability conditions: Numerical results

###  3.1 Dependence on the poloidal and toroidal magnetic energy

Figure 2: Time evolution of (a) the total magnetic energy (normalized to its
initial value) and (b) the ratio between total kinetic energy and magnetic
energy for the initial magnetic field configuration labeled as ‘Field I’ and
different initial values of $E_{\rm tor}/E_{\rm pol}$. For all cases, the
total initial magnetic energy is $E_{\rm mag}=2.8\times 10^{-3}|E_{\rm
grav}|$. Figure 3: Snapshots of the magnetic field configurations at the
indicated times for the initial ‘Field I’ configuration and different initial
values of $E_{\rm tor}/E_{\rm pol}$. The $z$-axis of the simulation box was
aligned with the instantaneous magnetic axis $\vec{M}$, so the upper panels
show a meridional cut through the star, while the lower panels show an
equatorial cut. The color scale corresponds to the strength and direction of
the component of the magnetic field perpendicular to the plane shown, while
the magnetic field lines correspond to the magnetic field parallel to that
plane. The blue circle represents the star’s surface. The value of the
asymmetry parameter, $\mathcal{A}$, is specified at the top of the upper
panels.  Figure 4: Time evolution of the amplitudes of the azimuthal modes
($m=0$ to $m=5$) for the radial (first row), polar (second row), and azimuthal
(third row) components of the velocity field. The different columns correspond
to different initial values of $E_{\rm tor}/E_{\rm pol}$.

Figure 5: $E_{\rm tor}/|E_{\rm grav}|$–$E_{\rm pol}/|E_{\rm grav}|$ plane for
(a) Field I and (b) Field II magnetic field configurations in stably
stratified stellar models with $\gamma=4/3$. Red points represent stable
configurations, while light and dark blue points correspond to unstable
configurations. The dot-dashed and solid lines represent the upper and lower
limits in Equation (4), respectively. Gray dashed lines are lines of constant
ratio $E_{\rm mag}/|E_{\rm grav}|$, with values as indicated.

In our first set of simulations, we constructed a stably stratified star with
$\gamma=4/3$ (and $\Gamma=5/3$), using the Field I magnetic configuration.
Figure 2 shows the time evolution of the magnetic energy (left panel) and the
ratio of the kinetic energy to the magnetic energy (right panel) for some
representative initial ratios between the energies in the toroidal and
poloidal magnetic field components (adjusting the coefficients $a_{\rm pol}$
and $a_{\rm tor}$ of Equation [5]). In all these cases, the initial Alfvén
travel time is $\tau_{\rm A,0}=9.02\,\tau_{s}$ and the initial magnetic energy
is $E_{\rm mag}=2.8\times 10^{-3}|E_{\rm grav}|$, for which the stability
limits of Equation (4) can be written as $0.25\lesssim E_{\rm tor}/E_{\rm
pol}\lesssim 23$. Figure 3 shows snapshots of the magnetic field configuration
after several Alfvén time scales for some of the instances shown in Figure 2.

For the cases within the stability boundaries, namely $E_{\rm tor}/E_{\rm
pol}=18.99$, and $1.0$ (and several others not shown here), the magnetic
energy decays only very slowly, and, for $t\gtrsim{\textbf{0.1}}\tau_{\rm
A,0}$, the ratio $E_{\rm kin}/E_{\rm mag}$ also drops (almost) monotonically,
suggesting a small adjustment to a nearby stable equilibrium. For these cases,
the asymmetry parameter remains small: $\mathcal{A}<10^{-3}$ (see central
panel of Figure 3).

For the mostly poloidal cases with $E_{\rm tor}/E_{\rm pol}\lesssim 0.2$,
after $\sim 10$ Alfvén times, there is a fast decay of the magnetic energy
accompanied by a peak in the kinetic energy, signaling an instability.
Afterwards, the magnetic energy becomes almost constant and the kinetic energy
drops. As seen in Figure 3, for the cases with small, but non-zero $E_{\rm
tor}$, we find the magnetic field to evolve to a new non-axisymmetric magnetic
field configuration in a ‘tennis ball’ shape, as in Braithwaite (2009) and
Mitchell et al. (2015).

For the other extreme of mostly and purely toroidal fields, $E_{\rm
tor}/E_{\rm pol}\gtrsim 150$, the magnetic energy at some time
($t\sim\mathrm{few}\times 10\tau_{\mathrm{A},0}$) starts to decay
continuously, accompanied by an increase in the kinetic energy, while the
field structure becomes increasingly asymmetric. Although the ratio $E_{\rm
kin}/E_{\rm mag}$ decreases at late times, $E_{\rm mag}$ keeps decreasing as
far as it was possible to run the simulations.444In Figure 2, the purely
toroidal field and the one with a very small poloidal component (which should
be stabilizing) seem to become unstable at essentially the same time. We have
checked this result with different initial perturbations in the velocity field
and found that the former configurations generally become unstable before the
latter, confirming our expectations. It might be worth noticing that in stably
stratified stars the instability sets in substantially later ($t\sim
8-20\tau_{A,0}$) than in barotropic stars ($t\sim\tau_{A,0}$; see below and
Paper I).

As done in Braithwaite (2006), we can separately follow the evolution of each
azimuthal mode, $m$, by performing a Fourier decomposition in the
$\phi$-direction on each component of the velocity field:

$\vec{u}(r,\theta,\phi,t)=\sum_{m=0}^{\infty}\vec{u}_{m}(r,\theta,t)e^{im\phi}\,,$
(12)

where

$\vec{u}_{m}(r,\theta,t)={1\over
2\pi}\int_{0}^{2\pi}d\phi\,\vec{u}(r,\theta,\phi,t)e^{-im\phi}.$ (13)

Figure 4 shows the amplitudes, $U_{m}^{i}(t)$, of the spherical components
($i=r,\theta,\phi$) of the velocity field (defined as the rms of
$|u_{m}^{i}(r,\theta,t)|$ over the stellar volume), for modes from $m=0$ to
$m=5$ for the same simulations of Figure 3. In all panels, we see an initial
noisy decay of the velocity amplitude. For the stable case ($E_{\rm
tor}=E_{\rm pol}$), this behavior continues until essentially no motion
occurs. For all the other cases, we later see an increase of (at least) some
of the modes, signaling an instability. As expected from the stable
stratification of the star (e. g., Tayler 1973; Akgün et al. 2013), the radial
component of the velocity field associated with these instabilities is much
smaller than the horizontal components ($U_{m}^{r}\ll U_{m}^{\theta}\sim
U_{m}^{\phi}$).

Goossens et al. (1981) showed that all purely toroidal magnetic fields are
dynamically unstable in the relatively small regions that simultaneously
satisfy $B_{\phi}=0$ and $\partial B_{\phi}^{2}/\partial\theta>0$, with the
fastest growing instabilities being those with $m=1$ azimuthal dependence,
which they find to have linear growth rates $\sim 1/t_{A}$ (or even faster for
particular configurations; Goossens & Veugelen 1978). For the purely or mostly
toroidal fields (first two columns of Figure 4), we confirm that the most
unstable mode is $m=1$, in agreement with the predictions of Tayler (1973) and
Goossens et al. (1981) and the simulation results of Braithwaite (2009).
Motivated by this result, we also ran simulations with an initially small
velocity field $\delta\vec{u}\propto\vec{\nabla}Y_{11}\times\vec{r}$, where
$Y_{\ell m}$ is the standard spherical harmonic function (see also Lander et
al., 2010). We found that this makes the instability appear at somewhat
earlier times ($\sim 10-20\tau_{A,0}$), and, for simulations with smaller
values of the viscosity coefficient, this time decreases to around
$9~{}\tau_{A,0}$ (we note, in any case, that the exponential growth time
inferred from the initial increase of the kinetic energy is several times
smaller than the time at which the unstable mode appears above the noise
level, as read from Figure 4). We checked for the places where the instability
starts to grow but they did not coincide with the ones predicted by Goossens
et al. (1981). Moreover, the instability seen in our simulations appears to be
global, involving a motion of the whole star (see Figure 2), clearly different
from the one studied by Goossens & Veugelen (1978) and Goossens et al. (1981).
A general study of the toroidal instability will be pursued in a forthcoming
work.

For the mostly poloidal field ($E_{\rm tor}/E_{\rm pol}=0.1$; fourth column of
Figure 4), the most unstable mode is $m=2$, while for the purely poloidal
field ($E_{\rm tor}/E_{\rm pol}=0$; last column of Figure 4), all azimuthal
modes appear to be similarly unstable, in agreement with the linear analysis
of Lander & Jones (2011).

We ran similar simulations (always for the ‘Field I’ magnetic field
configuration) changing the initial strength of the magnetic field to $E_{\rm
mag}/|E_{\rm grav}|=1.2\times 10^{-3}$, $0.011$, and $0.025$, which correspond
to $\tau_{\rm A,0}/\tau_{s}=13.5$, $4.5$, and $3.0$, respectively. Figure 5(a)
summarizes the results in the $E_{\rm tor}/|E_{\rm grav}|$–$E_{\rm
pol}/|E_{\rm grav}|$ plane. The types of behaviour found are as described in
the previous paragraphs, which allowed us to distinguish quite clearly which
of the initial magnetic field configurations were stable and which were not.
We can see that, for this set of simulations, the lower boundary of the stable
region defined by Equation (4) fairly accurately represents our numerical
results, whereas the upper boundary only gives a very rough approximation.

The same study was done for the Field II magnetic field configurations, and
the results are summarized in Figure 5(b). We ran simulations with magnetic
energies $E_{\rm mag}/|E_{\rm grav}|=2.2\times 10^{-3}$, $8.9\times 10^{-3}$,
$0.023$, and $0.053$. In all these cases, Equation (4) describes both
boundaries of the stability region reasonably well.

Unfortunately, simulations with higher ratios $E_{\rm mag}/|E_{\rm grav}|$
suffer from numerical instabilities that do not allow us to follow their
evolution and check the stability of very strong magnetic fields.

###  3.2 Dependence on the star’s stable stratification

Finally, we explored the accuracy of Equation (4) for different degrees of
stable stratification of the star. Figure 6 shows the evolution of the
magnetic energy for $\Gamma=5/3$ and an initial Field I configuration with
$E_{\rm tor}=E_{\rm pol}$, $\tau_{\rm A,0}\approx 9~{}\tau_{s}$, and $E_{\rm
mag}\approx 2.8\times 10^{-3}~{}|E_{\rm grav}|$ for different values of
$\gamma$. For the cases when $\Gamma>\gamma$, even for the very slightly non-
barotropic one ($\gamma=1.61$), the magnetic energy decays on a much longer
timescale than the Alfvén timescale, indicating stability, consistent with
Equation (4).

Figure 6 also shows the evolution of the magnetic energy for Fields I and II
when $\gamma=\Gamma=5/3$ (barotropic case). In both cases, the magnetic energy
quickly drops to half its initial value at about $10~{}\tau_{A,0}$. As seen in
Figure 7, by this time the magnetic field has moved to the star’s surface,
roughly preserving its axial symmetry ($\mathcal{A}<0.1$). Unfortunately,
beyond this point, the numerical inaccuracies increase considerably (the
energy balance equation is no longer satisfied with good accuracy) and can no
longer be neglected, making the simulation unreliable. We expect that, due to
the magnetic buoyancy, a large portion of the magnetic field will get expelled
from the star and eventually be dissipated by the magnetic diffusivity of the
atmosphere. We believe that these numerical problems originate near the
surface of the star, where the transition to the uniform temperature
atmosphere happens. They do not appear in the stably stratified cases
($\Gamma>\gamma$), even in the slightly stratified ones, probably because the
magnetic field does not move to the star’s surface, as it does in the
barotropic case. Moreover, we recall that, due to our numerical setup (see
Section 2.1), the star is not completely barotropic, but has a thin stably
stratified surface layer that could be preventing the magnetic field from
being completely expelled.

Figure 8 shows the time evolution of the amplitudes of the azimuthal modes in
a barotropic star with the initial Field I configuration with $E_{\rm
tor}/E_{\rm pol}=\infty$. In contrast with the stably stratified case (first
column of Figure 4), a clear instability appears at $\sim\tau_{A,0}$, with all
velocity components of similar amplitude, $U_{1}^{r}\sim U_{1}^{\theta}\sim
U_{1}^{\phi}$.

We ran additional simulations for different values of $E_{\rm tor}/E_{\rm
pol}$. Figure 9 summarizes the results of these simulations on the $(E_{\rm
tor}/E_{\rm pol})$ vs. $(\Gamma/\gamma-1)(|E_{\rm grav}|/E_{\rm mag})$ plane
for the Field I and Field II magnetic field configurations, showing Equation
(4) to be more accurate for the latter.

Figure 6: Evolution of the total magnetic energy for different specific
entropy gradients. The initial magnetic field configuration corresponds to
Field I with $E_{\rm tor}=E_{\rm pol}$, and $E_{\rm mag}=2.8\times
10^{-3}~{}|E_{\rm grav}|$. We also show for reference the case of Field II for
$\gamma=5/3$.  Figure 7: Snapshots of the magnetic field for a barotropic star
($\Gamma=\gamma$) at the indicated times for the initial ‘Field I’(left panel)
and ‘Field II’ (right panel) configurations with $E_{\rm tor}=E_{\rm pol}$,
and $E_{\rm mag}=2.8\times 10^{-3}~{}|E_{\rm grav}|$. The colour scale and
field lines are the same as in Fig. 3.  Figure 8:  Same as Figure 4 for a
barotropic star ($\Gamma=\gamma=5/3$) with the initial Field I configuration
with $E_{\rm tor}/E_{\rm pol}=\infty$.

Figure 9: Plane of $x\equiv(E_{\rm tor}/E_{\rm pol})$ vs.
$(\Gamma/\gamma-1)(|E_{\rm grav}|/E_{\rm mag})$ for (a) Field I and (b) Field
II and different values of the polytropic index $\gamma=1+1/n$, always with
adiabatic index $\Gamma=5/3$ and $E_{\rm mag}=2.8\times 10^{-3}~{}|E_{\rm
grav}|$. Filled circles represent configurations found numerically to be
stable, while crosses signal unstable configurations. The dash-dotted and
solid lines correspond to the lower and upper limit, respectively, for
stability established in Equation (4).

##  4 Discussion and conclusions

We did a parameter-space study to explore the validity of the proposed
stability condition for axially symmetric magnetic fields in stars given by
Equation (4). It is worth noticing that we have used the Cowling approximation
to speed up the simulations, which has been shown to lead to an overestimation
of the system’s stability (Tayler, 1973). We verified the stability of axially
symmetric magnetic fields even in very slightly non-barotropic stars
($\Gamma=5/3$ and $\gamma=1.61$), but we did not find any stable magnetic
field configurations in barotropic stars ($\Gamma=\gamma=5/3$). These
simulations confirm previous suggestions that the star’s stable stratification
is a crucial ingredient for the stability of magnetic fields inside stars
(Reisenegger, 2009; Mitchell et al., 2015). Of course, neither the previous
simulations nor those reported here can prove that there are no stable
equilibria in barotropic stars. If there were any, they should satisfy the
Grad-Shafranov equation (Grad & Rubin, 1958; Shafranov, 1958). As mentioned in
the Introduction, when the magnetic field is not totally confined inside the
star, the numerical solutions of this equation found by various authors all
satisfy $E_{\rm tor}<0.1E_{\rm pol}$, so they lie below the lower limit of the
stability region defined by equation 4. We saw in the simulations that, in
stably stratified stars, such configurations evolve into non-axisymmetric
stable equilibria, while for barotropic stars the magnetic field moves to the
star’s surface, where our numerical setup becomes unreliable, so we cannot
follow the further evolution of the magnetic field.

We verified that the lower limit for $E_{\rm tor}/E_{\rm pol}$ in Equation (4)
is quite robust, a toroidal field with $E_{\rm tor}\gtrsim 0.25E_{\rm pol}$ is
needed to stabilize the poloidal field, independent of the magnetic field
configuration or the star’s stratification, in agreement with the earlier
analytical work of Wright (1973); Markey & Tayler (1973) and the simulations
of Braithwaite (2009).

On the other hand, the upper bound for $E_{\rm tor}/E_{\rm pol}$ in Equation
(4) appears to be somewhat less robust, as we found some numerically stable
configurations beyond it. In addition to the star’s stable stratification and
gravitational energy (considered in Equation [4]), this limit appears to
depend on the geometry of the magnetic field. Additionally, many works have
established that rotation reduces the growth rate of the Tayler instability
(see, e.g., Pitts & Tayler, 1985; Spruit, 1999) and might have a similar
effect on others instabilities.

Extrapolating Equation (4) to very strong magnetic fields (where its
applicability is questionable) implies an upper bound for the total magnetic
energy:

$E_{\rm mag}\lesssim 5(\Gamma/\gamma-1)|E_{\rm grav}|,$ (14)

with comparable fractions in the poloidal and toroidal components, $E_{\rm
pol}^{\rm max}\sim 4E_{\rm tor}^{\rm max}$. For the Ap/Bp stars, the fluid can
be approximated as a classical monatomic ideal gas with $\Gamma=5/3$ and
$\gamma=4/3$ (MacGregor & Cassinelli, 2003), so $E_{\rm mag}\lesssim{5\over
4}|E_{\rm grav}|$, very similar to the rigorous upper limit implied by the
virial theorem, $E_{\rm mag}\leq|E_{\rm grav}|$. From this condition, the
maximum allowed value for $B_{\rm rms}$ in Ap/Bp stars is $\sim 5\times
10^{7}$ G. For degenerate stars, since $\Gamma/\gamma-1\ll 1$, Equation (14)
is more constraining than the virial theorem. In white dwarfs,
$\Gamma/\gamma-1\sim T_{7}/500$ with $T=T_{7}\times 10^{7}$ K (Reisenegger,
2009), implying $B_{\rm rms,max}\sim 10^{11}T_{7}^{1/2}$ G, whereas in neutron
stars $\Gamma/\gamma-1\sim Y$ ($Y$ is the proton fraction; Reisenegger &
Goldreich 1992), so $B_{\rm rms,max}\sim 3\times
10^{17}\left(Y/0.1\right)^{1/2}$ G. The observed surface magnetic fields in
all these stars are far below these limits.

We note that our application of the present results to neutron stars assumes
that they are composed of a “normal” fluid, that is, ignoring the formation of
a solid crust and the transition of protons and neutrons to a superconducting
and superfluid state, respectively. Most of the crust (its inner regions)
freezes after a few years (Krüger et al., 2015); while protons are expected to
become superconducting in a few months (or even faster, if direct URCA
processes are allowed; Gnedin et al. 2001), and neutrons to become superfluid
in a few hundred years (Shternin et al., 2011). The crust has a stabilizing
effect on the magnetic field, which can be important unless the field strength
is extremely high (e. g., Lander & Gourgouliatos 2019 and references therein).
The transition of the protons to a superconducting state causes the
confinement of the magnetic field into magnetic flux tubes (for type-II
superconductors) or other small-scale structures (for type-I superconductors),
increasing the magnetic energy and thus the Lorentz force for a given magnetic
flux (Easson & Pethick, 1977). In the absence of rotation, the superfluid
transition of the neutrons decouples them from the charged particles, and thus
from the magnetic field (Gusakov et al., 2020). In low-density regions of the
core, where protons and electrons are the only charged particles, they act as
a barotropic fluid, which cannot by itself stabilize the magnetic field.
However, at higher densities, there are additional charged particles, such as
muons, whose relative abundance is a function of density, making the charged-
particle fluid stably stratified. In this case, the magnetic field could be
stabilized by this stably stratified fluid, making the results of this paper
at least qualitatively applicable. The magnetic field should be non-zero also
in the barotropic region, but there it will be constrained to satisfy the
Grad-Shafranov equation. In a rotating neutron star, the quantized neutron
vortices interact with the magnetic field, providing a coupling between
neutrons and charged particles whose effect remains to be understood in detail
(Dommes & Gusakov, 2021).

## Acknowledgements

The authors are grateful to J. Braithwaite, F. Castillo, N. Moraga, D.
Ofengeim, and H. Spruit for useful discussions, and to the Referee for
stimulating comments and suggestions. This work was supported by FONDECYT
projects 3190172 (L.B.), 1201582 (A.R.), and 1190703 (J.A.V.), and the Center
for Astrophysics and Associated Technologies (CATA; ANID Basal project
FB210003). J.A.V. thanks for the support of CEDENNA under ANID/CONICYT grant
FB0807. M.E.G. acknowledges support from the Russian Science Foundation [Grant
No. 22-12-00048]. The calculations presented in this work were performed on
the Geryon computer at the UC Center for Astro-Engineering, part of the BASAL
PFB-06 and AFB-170002 projects, which received additional funding from QUIMAL
130008 and Fondequip AIC-57 for upgrades.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Akgün et al. (2006) Akgün T., Link B., Wasserman I., 2006, MNRAS, 365, 653
  * Akgün et al. (2013) Akgün T., Reisenegger A., Mastrano A., Marchant P., 2013, MNRAS, 433, 2445
  * Armaza et al. (2015) Armaza C., Reisenegger A., Valdivia J. A., 2015, ApJ, 802, 121
  * Ashton et al. (2017) Ashton G., Jones D. I., Prix R., 2017, MNRAS, 467, 164
  * Aurière et al. (2007) Aurière M., et al., 2007, A&A, 475, 1053
  * Becerra et al. (2022) Becerra L., Reisenegger A., Valdivia J. A., Gusakov M. E., 2022, MNRAS, 511, 732
  * Braithwaite (2006) Braithwaite J., 2006, A&A, 453, 687
  * Braithwaite (2009) Braithwaite J., 2009, MNRAS, 397, 763
  * Braithwaite & Nordlund (2006) Braithwaite J., Nordlund Å., 2006, A&A, 450, 1077
  * Braithwaite & Spruit (2004) Braithwaite J., Spruit H. C., 2004, Nature, 431, 819
  * Chandrasekhar & Fermi (1953) Chandrasekhar S., Fermi E., 1953, ApJ, 118, 116
  * Chandrasekhar & Prendergast (1956) Chandrasekhar S., Prendergast K. H., 1956, Proceedings of the National Academy of Science, 42, 5
  * Ciolfi et al. (2009) Ciolfi R., Ferrari V., Gualtieri L., Pons J. A., 2009, MNRAS, 397, 913
  * Cowling (1941) Cowling T. G., 1941, MNRAS, 101, 367
  * Cowling (1945) Cowling T. G., 1945, MNRAS, 105, 166
  * Cutler (2002) Cutler C., 2002, Phys. Rev. D, 66, 084025
  * Dall’Osso et al. (2009) Dall’Osso S., Shore S. N., Stella L., 2009, MNRAS, 398, 1869
  * Dommes & Gusakov (2021) Dommes V. A., Gusakov M. E., 2021, Phys. Rev. D, 104, 123008
  * Donati & Landstreet (2009) Donati J.-F., Landstreet J. D., 2009, ARA&A, 47, 333
  * Duez & Mathis (2010) Duez V., Mathis S., 2010, A&A, 517, A58
  * Easson & Pethick (1977) Easson I., Pethick C. J., 1977, Phys. Rev. D, 16, 275
  * Flowers & Ruderman (1977) Flowers E., Ruderman M. A., 1977, ApJ, 215, 302
  * Fujisawa et al. (2012) Fujisawa K., Yoshida S., Eriguchi Y., 2012, MNRAS, 422, 434
  * Gnedin et al. (2001) Gnedin O. Y., Yakovlev D. G., Potekhin A. Y., 2001, MNRAS, 324, 725
  * Goossens & Veugelen (1978) Goossens M., Veugelen R., 1978, A&A, 70, 277
  * Goossens et al. (1981) Goossens M., Biront D., Tayler R. J., 1981, Ap&SS, 75, 521
  * Grad & Rubin (1958) Grad H., Rubin H., 1958, in Proceedings of the Second United Nations International Conference on the Peaceful Uses of Atomic Energy. p. 190
  * Gusakov et al. (2020) Gusakov M. E., Kantor E. M., Ofengeim D. D., 2020, MNRAS, 499, 4561
  * Haskell et al. (2008) Haskell B., Samuelsson L., Glampedakis K., Andersson N., 2008, MNRAS, 385, 531
  * Krüger et al. (2015) Krüger C. J., Ho W. C. G., Andersson N., 2015, Phys. Rev. D, 92, 063009
  * Lander & Gourgouliatos (2019) Lander S. K., Gourgouliatos K. N., 2019, MNRAS, 486, 4130
  * Lander & Jones (2009) Lander S. K., Jones D. I., 2009, MNRAS, 395, 2162
  * Lander & Jones (2011) Lander S. K., Jones D. I., 2011, MNRAS, 412, 1730
  * Lander & Jones (2012) Lander S. K., Jones D. I., 2012, MNRAS, 424, 482
  * Lander et al. (2010) Lander S. K., Jones D. I., Passamonti A., 2010, MNRAS, 405, 318
  * MacGregor & Cassinelli (2003) MacGregor K. B., Cassinelli J. P., 2003, ApJ, 586, 480
  * Markey & Tayler (1973) Markey P., Tayler R. J., 1973, MNRAS, 163, 77
  * Mereghetti (2008) Mereghetti S., 2008, A&ARv, 15, 225
  * Mitchell et al. (2015) Mitchell J. P., Braithwaite J., Reisenegger A., Spruit H., Valdivia J. A., Langer N., 2015, MNRAS, 447, 1213
  * Pencil Code Collaboration et al. (2021) Pencil Code Collaboration et al., 2021, The Journal of Open Source Software, 6, 2807
  * Pitts & Tayler (1985) Pitts E., Tayler R. J., 1985, MNRAS, 216, 139
  * Prendergast (1956) Prendergast K. H., 1956, ApJ, 123, 498
  * Reiners (2009) Reiners A., 2009, in Strassmeier K. G., Kosovichev A. G., Beckman J. E., eds, Vol. 259, Cosmic Magnetic Fields: From Planets, to Stars and Galaxies. pp 339–344 (arXiv:0901.1719), doi:10.1017/S1743921309030701
  * Reisenegger (2009) Reisenegger A., 2009, A&A, 499, 557
  * Reisenegger & Goldreich (1992) Reisenegger A., Goldreich P., 1992, ApJ, 395, 240
  * Shabaltas & Lai (2012) Shabaltas N., Lai D., 2012, ApJ, 748, 148
  * Shafranov (1958) Shafranov V. D., 1958, Soviet Journal of Experimental and Theoretical Physics, 6, 545
  * Shternin et al. (2011) Shternin P. S., Yakovlev D. G., Heinke C. O., Ho W. C. G., Patnaude D. J., 2011, MNRAS, 412, L108
  * Sikora et al. (2019) Sikora J., Wade G. A., Power J., Neiner C., 2019, MNRAS, 483, 3127
  * Spruit (1999) Spruit H. C., 1999, A&A, 349, 189
  * Stairs et al. (2000) Stairs I. H., Lyne A. G., Shemar S. L., 2000, Nature, 406, 484
  * Tayler (1973) Tayler R. J., 1973, MNRAS, 161, 365
  * Thompson & Duncan (1996) Thompson C., Duncan R. C., 1996, ApJ, 473, 322
  * Viganò & Pons (2012) Viganò D., Pons J. A., 2012, MNRAS, 425, 2487
  * Wasserman (2003) Wasserman I., 2003, MNRAS, 341, 1020
  * Wasserman et al. (2022) Wasserman I., Cordes J. M., Chatterjee S., Batra G., 2022, ApJ, 928, 53
  * Wright (1973) Wright G. A. E., 1973, MNRAS, 162, 339
  * Yoshida & Eriguchi (2006) Yoshida S., Eriguchi Y., 2006, ApJS, 164, 156

##  Appendix A Equations

The MHD equations solved in the simulations presented in this paper are:

$\displaystyle\frac{\partial({\rm ln}\,\rho)}{\partial t}$ $\displaystyle=$
$\displaystyle-\vec{\nabla}\cdot\vec{u}-\vec{u}\cdot\vec{\nabla}({\rm
ln}\,\rho)$ (15) $\displaystyle\frac{\partial\vec{u}}{\partial t}$
$\displaystyle=$
$\displaystyle-\vec{u}\cdot\vec{\nabla}\vec{u}-\frac{\vec{\nabla}p}{\rho}-\vec{\nabla}\Phi+\frac{\vec{j}\times\vec{B}}{\rho}+\vec{f}_{\rm
visc}$ (16) $\displaystyle\frac{\partial\vec{A}}{\partial t}$ $\displaystyle=$
$\displaystyle\vec{u}\times\vec{B}-\eta\mu_{0}\vec{j}$ (17)
$\displaystyle\frac{\partial s}{\partial t}$ $\displaystyle=$
$\displaystyle-\vec{u}\cdot\vec{\nabla}s+\frac{\eta\mu_{0}|\vec{j}|^{2}}{\rho
T}+\frac{2\nu\mathbf{S}^{2}}{T}\,,$ (18)

where $\vec{A}$ is the magnetic vector potential,
$\vec{B}=\vec{\nabla}\times\vec{A}$ is the magnetic field,
$\vec{j}=\mu_{0}^{-1}\vec{\nabla}\times\vec{B}$ is the current density, and
$\mu_{0}$ the magnetic vacuum permeability. Here, $\Phi$ is the gravitational
potential, $\rho$ is the fluid mass density, $p$ is its pressure, $\vec{u}$ is
the fluid velocity, and $s$ is the entropy per unit mass. The viscous force is
given by

$\vec{f}_{\rm
visc}=\rho^{-1}\vec{\nabla}\cdot\left(2\rho\nu{\mathbf{S}}\right)\,,$ (19)

with $\nu$ as the kinematic viscosity, and $\mathbf{S}$ as the rate-of-shear
tensor whose components are

$S_{ij}=\frac{1}{2}\left[\frac{\partial u_{j}}{\partial x_{i}}+\frac{\partial
u_{i}}{\partial
x_{j}}-\frac{2}{3}\delta_{ij}(\vec{\nabla}\cdot\vec{u})\right]\,.$ (20)

We assume an ideal equation of state:

$p(\rho,s)=(\mathcal{R}/\mu)\rho T(\rho,s)\,,$ (21)

with,

$T(\rho,s)=T_{c}\left(\frac{\rho}{\rho_{c}}\right)^{\Gamma-1}\exp{[(s-s_{c})/c_{V}]}\,,$
(22)

where $s_{c}$, $\rho_{c}$ and $T_{c}$ are the entropy, density, and
temperature at the center of the star, respectively.