# Magnetohydrodynamic stability of magnetars in the ultrastrong field regime
II: The crust

Peter B. Rau1 and Ira Wasserman2
1Institute for Nuclear Theory, University of Washington, Seattle, WA, U.S.A.
2Cornell Center for Astrophysics and Planetary Science, Cornell University,
Ithaca, NY, U.S.A E-mail<EMAIL_ADDRESS>
###### Abstract

We study the stability of Hall MHD with strong magnetic fields in which Landau
quantization of electrons is important. We find that the strong-field Hall
modes can be destabilized by the dependence of the differential magnetic
susceptibility on magnetic field strength. This instability is studied using
linear perturbation theory, and is found to have typical growth time of order
$\lesssim 10^{3}$ yrs, with the growth time decreasing as a function of
wavelength of the perturbation. The instability is self-limiting, turning off
following a period of local field growth by a few percent of the initial
value. Finite temperature is also shown to limit the instability, with
sufficiently high temperatures eliminating it altogether. Alfvén waves can
show similar unstable behaviour on shorter timescales. We find that Ohmic
heating due to the large fields developed via the instability and magnetic
domain formation is not large enough to account for observed magnetar surface
temperatures. However, Ohmic heating is enhanced by the oscillatory
differential magnetic susceptibility of Landau-quantized electrons, which
could be important to magneto-thermal simulations of neutron star crusts.

###### keywords:

stars: neutron – stars: magnetars – stars: magnetic fields – instabilities –
MHD

††pubyear: 2022††pagerange: Magnetohydrodynamic stability of magnetars in the
ultrastrong field regime II: The crust–References

##  1 Introduction

Magnetars are a class of strongly-magnetized neutron stars with typical
surface fields $\sim 10^{14}$–$10^{15}$ G distinguished by their unique X-ray
pulsations and repeating bursts of gamma rays (Mereghetti et al., 2015;
Turolla et al., 2015; Kaspi & Beloborodov, 2017). One class of proposed
mechanisms behind the bursts and more powerful but rarer giant flares involves
magnetohydrodynamic (MHD) instabilities inside the magnetar leading to sudden
release of energy into the surrounding magnetosphere through the excitation of
Alfvén waves along the magnetic field lines anchored to the stellar crust
(Thompson & Duncan, 1995; Heyl & Hernquist, 2005). A second class involves a
slower build-up of magnetospheric twisting leading to reconnection and the
creation of a magnetically-confined pair plasma (Lyutikov, 2003, 2006). The
exact mechanism for triggering these instabilities or generating the
magnetospheric twists is not known, but both are associated with magnetic
field evolution and crust motions or mechanical failures. As significant
cracking and slipping of the crust is forbidden by the high pressures inside a
neutron star (Jones, 2003; Levin & Lyutikov, 2012), the failure mechanism of
the crust cannot be analogous to earthquakes. A more plausible possibility is
that the crust undergoes plastic deformation once the magnetic stresses exceed
a particular yielding value– this has been investigated in the context of
thermoplastic waves (Beloborodov & Levin, 2014) and Hall wave avalanches (Li
et al., 2016) which have been suggested as the source of magnetar outbursts.
MHD instabilities are also invoked in many classes of fast radio burst (FRB)
models with magnetar progenitors (e.g., Metzger et al. (2019); Lu et al.
(2020) and references therein), where they are responsible for excitation of
Alfvén waves which travel into the magnetosphere or for releasing ejecta
during flaring events.

An additional aspect of magnetar crust physics that is not fully understood is
the anomalously high surface temperatures inferred from the thermal luminosity
of many magnetars. They are observed to be significantly hotter than other
isolated neutron stars, with typical luminosities $\mathscr{L}\approx 10^{35}$
erg/s, roughly two orders of magnitude larger than those of standard rotation-
powered pulsars. Many detailed simulations of their magneto-thermal evolution,
including transport accounting for Landau quantization of fermions (Hernquist,
1984; Potekhin, 1999; Potekhin & Yakovlev, 2001) and different accreted
envelopes (oceans) have been performed to try to understand this problem. Pons
et al. (2009) found increased heating in magnetars because of enhanced Joule
heating from decay of crustal currents at strong magnetic field strengths.
Viganò et al. (2013) found that the high luminosities of magnetars could be
explained by high conductivity light element envelopes. Simulations by
Potekhin & Chabrier (2018) including transport by Landau-quantized electrons
were able to explain the temperatures of about half the known magnetars, but
discrepancies between the theoretical luminosity prediction and observation of
around an order of magnitude remain for the other half. One possible
explanation for this disagreement is an unaccounted-for heat source from
dissipation of the magnetic field. A phenomenological calculation by Kaminker
et al. (2006) of interior heating in a magnetar argued that a heat source
located within the outer crust of the star would be consistent with the
observed high surface temperatures, as heat energy from a source deeper in the
star would be cooled by neutrino emission before warming the surface.
Beloborodov & Li (2016) reviewed four different heating mechanisms that could
be responsible for these hotter temperatures, but found none that were
completely satisfactory in explaining both the observed surface luminosities
and magnetar lifetimes $\approx 10^{4}$ yrs. Akgün et al. (2018) proposed that
Ohmic dissipation of currents in weakly-conducting magnetar envelopes could be
sufficient to explain observed surface temperatures. More recently, Chamel et
al. (2021) have proposed electron capture and pycnonuclear fusion reactions
triggered by magnetic field evolution within the crust as a viable heating
source to explain the theory-observation discrepancy.

In an idealized neutron star crust, only the electrons (and in the inner
crust, the nonconducting superfluid neutrons) are free to move. The electrons
thus carry the electric current, and the magnetic field moves with them with
respect to the neutralizing nuclear lattice. Under these conditions, the ideal
magnetohydrodynamics (MHD) applicable in the core of a neutron star must be
supplanted with Hall MHD, and the Hall term included in Ohm’s Law. Hall MHD
and Hall equilibria, the equilibrium field configurations when Hall MHD
governs the evolution of the field, have been studied for neutron star crusts
in many previous calculations (e.g., Jones (1988); Goldreich & Reisenegger
(1992); Rheinhardt & Geppert (2002); Cumming et al. (2004); Gourgouliatos et
al. (2013); Gourgouliatos & Cumming (2015); Lyutikov (2015); Li et al.
(2016)). A notable difference between ideal MHD and Hall MHD is that the
latter lacks a (canonical) energy principle to determine stability: this was
shown for electron MHD by Lyutikov (2013), and this conclusion can be extended
to the closely-related Hall MHD. Numerical MHD simulations including the Hall
term have demonstrated that, by generating small-scale spatial structures, it
greatly enhances dissipation of magnetic field energy during the early stages
of a neutron star’s lifetime (Viganò et al., 2013).

In Rau & Wasserman (2021) (henceforth Paper I), we studied MHD stability in
ultrastrong magnetic fields relevant to magnetar cores using ideal MHD and a
canonical energy principle (Friedman & Schutz, 1978; Glampedakis & Andersson,
2007). We included the magnetic field-dependence of the internal energy due to
Landau quantization of fermions, showing how this can lead to a fast-growing,
but spatially limited, instability. In this paper, we extend that study to
magnetar crusts, replacing our ideal MHD formalism with Hall MHD. This is the
first study including Landau quantization effects in Hall MHD. As mentioned,
we cannot use the canonical energy principle computed in Paper I, Eq. 8 to
analyse stability in the crust. Instead we use perturbation theory to look for
a strong-field Hall MHD instability. This is done through a linear mode
analysis with coupling to the crust following Cumming et al. (2004). We first
look at this instability at zero temperature, then include finite temperature
and examine how temperature can limit the instability. In a subsequent paper
we will examine Hall MHD including Landau quantization numerically to examine
nonlinear effects. One of our goals is to examine the effect that Landau
quantization has on stability and the potential to trigger magnetar outbursts.

In Section 2 we briefly outline the thermodynamics of a strongly magnetized
outer crust, referring to Paper I where appropriate. The results in this
section are used later in the stability analysis. Section 3 describes Hall MHD
and its linear perturbative modes, showing the mechanisms by which unstable
Hall waves (whistler modes) could develop in the crust. In Section 4 we study
the temperature dependence of this instability and determine the potentially
unstable region in field–density–temperature parameter space. In Section 5 we
examine the modification to Ohmic dissipation due to Landau quantization of
electrons, and estimate the Ohmic heating resulting from the dissipation of
the unstably growing magnetic field during magnetic domain formation. The
astrophysical relevance of the instability is discussed in that section and in
the conclusion. We work in Gaussian units and employ the Einstein summation
convention using Latin letters as spatial indices $i=1,2,3$.

##  2 Thermodynamics of a strongly magnetized crust

We consider a neutron star crust of nuclei and electrons. In the inner crust,
we make the simplifying assumption that the dripped neutron superfluid moves
with the nuclear lattice– this is equivalent to assuming strong entrainment
between the superfluid neutrons and the nuclei (Chamel, 2017). To determine
the MHD stability, we must compute certain thermodynamic derivatives which are
given below. This section is organized similarly to, and with analogous
content to, Section 4 of Paper I.

The total grand potential density $\Omega$ for the crust must include
contributions for the Landau-quantized electron gas, the magnetic pressure,
the nuclear lattice and the dripped neutron superfluid. As we will not be
overly concerned with the latter two contributions, we can write $\Omega$ as

$\Omega(\mu_{\text{e}},n_{\text{b}},Y_{\text{N}},w,A,Z,B,E^{s}_{ij})=-P_{\text{e}}(\mu_{\text{e}},B)+\frac{B^{2}}{8\pi}+\Omega_{\text{b}}(n_{\text{b}},Y_{\text{N}},w,A,Z,E^{s}_{ij}),$
(1)

where $\mu_{\text{e}}$ is the electron chemical potential, $n_{\text{b}}$ is
the total number density of baryons including those in the lattice nuclei and
the dripped neutrons, $Y_{\text{N}}$ is the fraction of baryons within the
nuclei, $w$ is the volume fraction of each unit cell occupied by the nucleus,
$A$ and $Z$ are the atomic and mass numbers of the nuclei, $E^{s}_{ij}$ is the
(deviatoric) shear strain tensor and $B$ is the magnetic field magnitude.
$P_{\text{e}}$ is the electron pressure given in e.g., Paper I, Eq. (40), and
$\Omega_{\text{b}}$ is the grand potential associated with the baryons i.e.,
the nuclei and the dripped neutrons, including lattice energy, surface energy
terms, elastic deformation, etc., which typically appear in treatments of the
crust (e.g., Baym et al. (1971a, b)). The dripped neutron contributions are
excluded in the outer crust. Though the natural independent variables for
$\Omega_{\text{b}}$ would be the chemical potentials associated with the
lattice nuclei and dripped neutrons, we assume that we can perform a change of
variable and write $\Omega_{\text{b}}$ in terms of the number densities as is
often done in equation of state (EOS) calculations. The vacuum
Euler–Heisenberg Lagrangian, included in $\Omega$ in Paper I, is not included
in this paper as its effects on the stability are small; its inclusion would
add an additional term to the vacuum field pressure term $B^{2}/(8\pi)$. We
have neglected temperature in this expression, as its main role at
temperatures relevant to most of the lifetime of a magnetar is to regulate
divergences in thermodynamic functions at filled Landau levels. If finite
temperature is included, Eq. (1) and the following equations can be easily
generalized to include it, though it adds unnecessary complication for our
purposes.

Charge neutrality requires that the electron number density $n_{\text{e}}$
equals $Yn_{\text{b}}$ where the proton fraction is $Y=Y_{\text{N}}Z/A$, but
this must be imposed after taking required partial derivatives of $\Omega$ or
the internal energy density $u$. The latter quantity is given by (noting that
temperature is neglected and hence $u$ and the Helmholtz free energy density
are equal)

$u(n_{\text{e}},n_{\text{b}},Y_{\text{N}},w,A,Z,B,E^{s}_{ij})=u_{\text{e}}(n_{\text{e}},B)+\frac{B^{2}}{8\pi}+u_{\text{b}}(n_{\text{b}},Y_{\text{N}},w,A,Z,E^{s}_{ij}),$
(2)

where $u_{\text{e}}$ and $u_{\text{b}}$ are the electron and baryonic (nuclei
plus dripped neutron) internal energy density respectively. $u$ and $\Omega$
are related by the Legendre transformation

$u=\Omega+n_{\text{e}}\mu_{\text{e}}+n_{\text{b}}\mu_{\text{b}},$ (3)

where $\mu_{\text{b}}=\partial u/\partial
n_{\text{b}}|_{Y_{\text{N}},w,A,Z,E_{ij}^{s}}$ is the chemical potential of
the baryons.

The first law of thermodynamics for both $u$ and $\Omega$ are used to compute
thermodynamic derivatives appearing in the MHD equations. Using arguments from
Paper I Section 4, we have $H=4\pi\partial u/\partial
B|_{n_{\text{e}}}=4\pi\partial\Omega/\partial B|_{\mu_{\text{e}}}$, and can
define $M\equiv-\partial u_{\text{e}}/\partial B|_{n_{e}}=-u_{B}$ where $M$ is
the usual definition of the magnetization and $u_{B}$ is the notation used in
Paper I. We then have

$\text{d}M=\left.\frac{\partial^{2}P_{\text{e}}}{\partial
B^{2}}\right|_{\mu_{\text{e}}}\text{d}B+\frac{\partial^{2}P_{\text{e}}}{\partial
B\partial\mu_{\text{e}}}\text{d}\mu_{\text{e}}=\left[\left.\frac{\partial^{2}P_{\text{e}}}{\partial
B^{2}}\right|_{\mu_{\text{e}}}-\left(\left.\frac{\partial^{2}P_{\text{e}}}{\partial\mu_{\text{e}}^{2}}\right|_{B}\right)^{-1}\left(\frac{\partial^{2}P_{\text{e}}}{\partial
B\partial\mu_{\text{e}}}\right)^{2}\right]\text{d}B+\left(\left.\frac{\partial^{2}P_{\text{e}}}{\partial\mu_{\text{e}}^{2}}\right|_{B}\right)^{-1}\frac{\partial^{2}P_{\text{e}}}{\partial
B\partial\mu_{\text{e}}}\text{d}n_{\text{e}}\equiv\chi_{n}\text{d}B+\mathcal{M}_{n}\text{d}n_{\text{e}},$
(4)

$\chi_{n}=-u_{BB}$ and $\mathcal{M}_{n}=-u_{Bn}$ in the notation of Paper I,
with the former being the differential magnetic susceptibility at fixed
$n_{\text{e}}$. Eq. (4) is used to conveniently compute $\chi_{n}$ and
$\mathcal{M}_{n}$ in terms of partial derivatives of $P_{\text{e}}$. We will
also need partial derivatives of $P_{\text{e}}$, in particular

$\left.\frac{\partial P_{\text{e}}}{\partial
B}\right|_{n_{\text{e}}}=\left.\frac{\partial P_{\text{e}}}{\partial
B}\right|_{\mu_{\text{e}}}-n_{\text{e}}\left(\left.\frac{\partial^{2}P_{\text{e}}}{\partial\mu_{\text{e}}^{2}}\right|_{B}\right)^{-1}\frac{\partial^{2}P_{\text{e}}}{\partial
B\partial\mu_{\text{e}}}=M-n_{\text{e}}\mathcal{M}_{n}.$ (5)

The expressions for the partial derivatives of $P_{\text{e}}$ with respect to
$B$ and $\mu_{\text{e}}$ can be found in the online supplemental material for
Paper I. The functional behaviour of $u_{BB}=-\chi_{n}$ is identical to
$u_{BB}$ defined there but here it only includes the electrons, and
$u_{Bn}=-\mathcal{M}_{n}$ is defined analogously to $u_{\rho B}$, though they
are not identical. We treat the electrons as a free Fermi gas, ignoring the
effects of the crystal lattice and electron band structure on their Landau
quantization (see Shoenberg (1984) for how band structure modifies this
picture).

##  3 Perturbative Hall MHD modes in strong magnetic fields

The MHD stability of the crust must be studied in a different framework than
in the core. This is because in the crust the electric current is due to the
electrons moving relative to the (approximately) stationary, neutralizing
nuclear lattice. For sufficiently weak fields, the ions are effectively
decoupled from the magnetic field, which moves with the electrons (Cumming et
al., 2004). Ohm’s Law thus must be modified to include the Hall term (e.g.,
Goldreich & Reisenegger (1992)):

$E_{i}=-\frac{1}{c}\epsilon_{ijk}v_{j}B_{k}+\frac{1}{n_{\text{e}}ec}\epsilon_{ijk}J^{j}_{e}B^{k}+\frac{1}{\sigma}J_{e,i}=-\frac{1}{c}\epsilon_{ijk}v^{\text{e}}_{j}B_{k}+\frac{1}{\sigma}J_{e,i}.$
(6)

$E_{i}$ is the electric field, $\epsilon_{ijk}$ is the Levi-Civita tensor,
$B^{i}$ is the magnetic field, and $v_{\text{e}}^{i}$ and $v^{i}$ are the
velocity of the electrons and nuclear lattice respectively. $\sigma$ is the
conductivity, which is assumed isotropic for simplicity here, but in strong
magnetic fields this is not true in general. MHD when this is used as Ohm’s
Law is known as Hall MHD. Additional terms accouting for e.g., the electron
inertia and the electron pressure gradient (Biskamp, 2000; Cramer, 2001) are
neglected. The (free) current density follows from Ampère’s Law

$J^{i}_{e}=e(Zn_{\text{N}}v^{i}-n_{\text{e}}v_{\text{e}}^{i})=\frac{c}{4\pi}\epsilon^{ijk}\nabla_{j}H_{k},$
(7)

where $Z$ the atomic number of the nuclei and $n_{\text{N}}$ the number
density of nuclei. We assume local charge neutrality
$Zn_{\text{N}}=n_{\text{e}}$.

First consider the case where the nuclear lattice is fixed $v^{i}=0$. The
magnetic induction equation
$\partial_{t}B^{i}=-c\epsilon^{ijk}\nabla_{j}E_{k}$ becomes

$\partial_{t}B^{i}=-\frac{c}{4\pi
e}\epsilon^{ijk}\nabla_{j}\left(\frac{1}{n_{\text{e}}}\epsilon_{k\ell
m}B^{m}\epsilon^{\ell
np}\nabla_{n}H_{p}\right)-\frac{c^{2}}{4\pi}\epsilon^{ijk}\nabla_{j}\left(\frac{1}{\sigma}\epsilon_{k\ell
m}\nabla^{\ell}H^{m}\right).$ (8)

Since the Ohmic dissipation timescale is much longer than the Hall timescale
for magnetars (Goldreich & Reisenegger, 1992), we ignore Ohmic dissipation in
our analytic calculation. The relevant magnetic induction equation is thus

$\partial_{t}B^{i}=-\frac{c}{4\pi
e}\epsilon^{ijk}\nabla_{j}\left(\frac{1}{n_{\text{e}}}\epsilon_{k\ell
m}B^{m}\epsilon^{\ell np}\nabla_{n}H_{p}\right)=\epsilon^{ijk}\epsilon_{k\ell
m}\nabla_{j}\left(v_{\text{e}}^{\ell}B^{m}\right).$ (9)

This says that the magnetic field is frozen to the electron fluid. Magnetic
fields that are static solutions to this equation i.e., $\partial_{t}B^{i}=0$,
are Hall equilibria.

As mentioned in the introduction, there is no canonical energy-based stability
criterion for Hall MHD (Lyutikov, 2013). Instead, we examine the local
stability of the low-frequency oscillation modes that are perturbations of Eq.
(9): the Hall modes/waves, also referred to as whistler modes (Cramer, 2001).
These are a class of low-frequency, circularly-polarized MHD modes. We first
consider the case where the nuclear lattice is decoupled from the electrons,
which is the original case considered by Goldreich & Reisenegger (1992). We
then consider the modifications to the Hall mode dispersion due to the
coupling to the nuclear lattice via the Lorentz force acting on the nuclei
i.e. the relaxation of the $v^{i}=0$ approximation– this was first considered
by Cumming et al. (2004).

###  3.1 No coupling to lattice

We take the Eulerian perturbation of Eq. (9), assuming uniform background
quantities $B^{i}$, $v_{\text{e}}$, $n_{\text{e}}$, $M$, etc. The perturbation
of $H_{k}$ is

$\delta H_{k}=(1-4\pi M/B)\delta
B_{k}-4\pi(\chi_{n}-M/B)\hat{B}_{k}(\hat{B}^{m}\delta
B_{m})-4\pi\mathcal{M}_{n}\hat{B}_{k}\delta n_{\text{e}},$ (10)

where we employ $\mathcal{M}_{n}\equiv\partial^{2}u_{\text{e}}/\partial
B\partial n_{\text{e}}$. We also assume incompressible perturbations and
uniform $n_{\text{e}}$, so $\delta n_{\text{e}}=0$. Assuming harmonic spatial-
temporal dependence $\delta B^{k}\propto\exp(ik_{j}x^{j}-i\omega t)$, we
obtain

$-i\omega\delta B^{i}=\frac{c(k_{j}H^{j})}{4\pi
en_{\text{e}}}\epsilon^{ijk}k_{j}\delta
B_{k}-\frac{c(k_{j}B^{j})}{en_{\text{e}}}\left(\chi_{n}-\frac{M}{B}\right)(\hat{B}_{m}\delta
B^{m})\epsilon^{ijk}k_{j}\hat{B}_{k},$ (11)

where $\hat{B}^{i}=B^{i}/B$. Contract this with $\epsilon_{\ell
ni}\hat{k}^{n}$, where $\hat{k}_{n}=k_{n}/k$, then insert the result back into
the original equation to eliminate $\epsilon^{ijk}k_{j}\hat{B}_{k}$.
Contracting the resulting expression with $\hat{B}_{i}$ gives

$\omega=\pm\frac{ck|k_{j}B^{j}|}{4\pi
en_{\text{e}}}\sqrt{\frac{H^{2}}{B^{2}}-4\pi\frac{H}{B}\left(\chi_{n}-\frac{M}{B}\right)\sin^{2}\theta_{B}}\approx\pm\omega_{\text{H}}\sqrt{1-4\pi\chi_{n}\sin^{2}\theta_{B}},$
(12)

$\hat{B}^{i}=B^{i}/B$ $\omega_{\text{H}}$ is the normal ($B=H$) Hall
frequency, $\cos\theta_{B}\equiv\hat{k}_{j}\hat{B}^{j}$, and we used that
$H\approx B$ to within a few percent even for fields above the quantum
critical field (e.g., Figure 4(a) of Paper I), implying $|M|/B\ll 1$. This
mode is clearly unstable if the term under the square root is negative, and as
shown in Paper I, this can clearly be true at locations in $B$–$n_{\text{e}}$
parameter space where additional Landau levels are filled. If
$\hat{B}_{m}\delta B^{m}=0$, Eq. (12) reduces to
$\omega=\omega_{\text{H}}\frac{H}{B}$: in this case, the modes are stable for
a spatially uniform background.

In the unstable regions of parameter space, the unstable Hall modes have an
approximate growth time (using Eq. (12) with $B\approx H$)

$\tau=\frac{1}{\text{Im}(\omega)}\sim 6\times
10^{3}|\sec\theta_{B}|\left(\frac{10^{15}\text{
G}}{B}\right)\left(\frac{n_{\text{e}}}{10^{-4}\text{
fm}^{-3}}\right)\left(\frac{R_{c}^{-1}}{k}\right)^{2}\left(\frac{1}{|1-4\pi\chi_{n}\sin^{2}\theta_{B}|}\right)^{1/2}\text{
yrs},$ (13)

where $R_{c}=10^{5}$ cm is the approximate thickness of the crust. This
instability should be contrasted with the unstable Hall modes first discussed
in Rheinhardt & Geppert (2002) and later by Rheinhardt et al. (2004), Cumming
et al. (2004) and Pons & Geppert (2010). That is sometimes called the Hall
instability or the Hall drift instability, and depends on a strong shearing
(background) electron velocity as the source of instability. It is a long-
wavelength instability, in contrast to the instability here which grows faster
for shorter wavelengths. The instability here will occur in any sufficiently
strongly magnetized crust at up to hundreds of spatially-restricted regions
where $1-4\pi\chi_{n}\sin^{2}\theta_{B}<0$. However, it is much slower than
the MPR-type instability in Paper I due to the Hall phase velocity
$\omega_{\text{H}}/k$ being much slower than the Alfvén velocity, which set
the characteristic instability timescales in this paper and Paper I
respectively. We note that regions where $\chi_{n}>1/(4\pi)$ are unstable to
magnetic domain formation; this is discussed in further detail in Section 5.2.

###  3.2 Coupling to lattice

Motion of the nuclear lattice is described by the Euler equation with strong-
field magnetization terms (Eq. (3) of Paper I). In the solid crust we
additionally include shear stresses, since for magnetar-strength magnetic
fields in the crust, the elasticity must be included to balance the Lorentz
force exerted by the magnetic field. Using Eq. (1) and $P\equiv
P_{\text{e}}+\Omega_{\text{b}}$, we can write the Euler equation as

$\rho(\partial_{t}v_{i}+v^{j}\nabla_{j}v_{i})+\nabla_{i}P+\rho\nabla_{i}\Phi=\nabla^{j}\left[\frac{B^{2}}{8\pi}g_{ij}-\frac{1}{4\pi}H_{k}B^{k}g_{ij}+\frac{1}{4\pi}H_{i}B_{j}\right]+\nabla^{j}\sigma_{ij}\equiv\nabla^{j}T^{B}_{ij}+\nabla^{j}\sigma_{ij}.$
(14)

where $g_{ij}$ is the metric tensor and the (deviatoric) shear stress tensor
is

$\sigma_{ij}=\check{\mu}\left(\nabla_{i}X_{j}+\nabla_{j}X_{i}-\frac{2}{3}\nabla_{k}X^{k}g_{ij}\right)=2\check{\mu}E_{ij}^{s}.$
(15)

$X^{i}$ is the local displacement of the crust from its relaxed state,
$\check{\mu}$ is the shear modulus, and $E_{ij}^{s}$ is the (deviatoric) shear
strain tensor. The shear modulus of a Coulombic nuclear lattice is (Strohmayer
et al., 1991)

$\check{\mu}\approx 0.1925e^{2}Z^{2/3}n_{\text{e}}^{4/3}=2.4\times
10^{28}\left(\frac{Z}{40}\right)^{2/3}\left(\frac{n_{\text{e}}}{10^{35}\text{
cm}^{-3}}\right)^{4/3}\text{ g s}^{-2}\text{ cm}^{-1},$ (16)

where we ignore thermal corrections. For comparison, the electron pressure is
$P_{\text{e}}\approx(3\pi^{2})^{1/3}\hbar cn_{\text{e}}^{4/3}/4\approx
53(40/Z)^{2/3}{\check{\mu}}$ if many Landau levels are occupied.

We have three quantities describing the perturbations: the Eulerian
perturbation of the magnetic field $\delta B^{i}$, and the Lagrangian
displacement fields for the electrons and lattice $\xi^{i}_{\text{e}}$ and
$\xi^{i}$ respectively. $\delta v^{i}=\partial_{t}\xi^{i}$ and $\delta
v^{i}_{\text{e}}=\partial_{t}\xi^{i}_{\text{e}}$ in the limit of zero
background velocities. The perturbed forms of three equations, Eq. (7,9,14),
determine these perturbations, so two of the perturbed quantities can be
eliminated.

Taking the Eulerian perturbation of Eq. (7) and using that
$v^{i}=0=v^{i}_{\text{e}}$ in the background,

$en_{\text{e}}(\delta v^{i}-\delta
v^{i}_{\text{e}})=\frac{c}{4\pi}\epsilon^{ijk}\nabla_{j}\delta H_{k}.$ (17)

Assuming $\xi^{i},\xi^{i}_{\text{e}}\propto\exp(ik_{j}x^{j}-i\omega t)$ and
using Eq. (10) and $\delta
n_{\text{e}}=-n_{\text{e}}\nabla_{j}\xi^{j}_{\text{e}}$ for uniform
$n_{\text{e}}$, we obtain

$-i\omega
en_{\text{e}}(\xi^{i}-\xi^{i}_{\text{e}})=\frac{ic}{4\pi}\left[\frac{H}{B}\epsilon^{ijk}k_{j}\delta
B_{k}-4\pi\left(\chi_{n}-\frac{M}{B}\right)\epsilon^{ijk}k_{j}\hat{B}_{k}(\hat{B}_{m}\delta
B^{m})+4\pi
n_{\text{e}}\mathcal{M}_{n}(k_{m}\xi^{m}_{\text{e}})\epsilon^{ijk}k_{j}\hat{B}_{k}\right].$
(18)

Contracting this equation with $k_{i}$ implies
$k_{j}\xi^{j}=k_{j}\xi^{j}_{\text{e}}$; using this condition to eliminate
$k_{j}\xi^{j}_{\text{e}}$ gives an expression for $\xi^{i}_{\text{e}}$ in
terms of $\xi^{i}$ and $\delta B^{i}$ alone:

$\xi^{i}_{\text{e}}=\xi^{i}+\frac{c\mathcal{M}_{n}}{e\omega}(k_{m}\xi^{m})\epsilon^{ijk}k_{j}\hat{B}_{k}+\frac{c}{4\pi
en_{\text{e}}\omega}\left[\frac{H}{B}\epsilon^{ijk}k_{j}\delta
B_{k}-4\pi\left(\chi_{n}-\frac{M}{B}\right)\epsilon^{ijk}k_{j}\hat{B}_{k}(\hat{B}_{m}\delta
B^{m})\right].$ (19)

For uniform background field the perturbed induction equation, using as
$E^{i}$ the second form of Eq. (6), we have

$\partial_{t}\delta B^{i}=B^{k}\nabla_{k}\delta
v_{\text{e}}^{i}-B^{i}\nabla_{k}\delta
v_{\text{e}}^{k}=B^{k}\nabla_{k}\partial_{t}\xi^{i}_{\text{e}}-B^{i}\nabla_{k}\partial_{t}\xi^{k}_{\text{e}}.$
(20)

Eliminating $\xi_{\text{e}}$ using Eq. (19) gives

$-i\omega\delta
B^{i}=(k_{j}B^{j})\omega\xi^{i}-B^{i}(k_{j}\xi^{j})\omega+\frac{c(k_{j}B^{j})}{e}\mathcal{M}_{n}(k_{m}\xi^{m})\epsilon^{ijk}k_{j}\hat{B}_{k}+\frac{c(k_{j}H^{j})}{4\pi
en_{\text{e}}}\epsilon^{ijk}k_{j}\delta
B_{k}+\frac{c(k_{j}B^{j})}{en_{\text{e}}}\left(\chi_{n}-\frac{M}{B}\right)(\hat{B}_{m}\delta
B^{m})\epsilon^{ijk}k_{j}\hat{B}_{k}$ (21)

Taking the Eulerian perturbation of Eq. (14) and using the magnetohydrostatic
form of Eq. (14), ignoring gravity and gradients of background quantities, we
find

$\partial_{t}^{2}\xi_{i}+\frac{1}{\rho}\nabla_{i}\delta
P=\frac{1}{\rho}\nabla^{j}\delta
T^{B}_{ij}+\frac{1}{\rho}\nabla^{j}\delta\sigma_{ij},$ (22)

where

$\displaystyle\delta T^{B}_{ij}=\frac{1}{4\pi}\left(B_{k}\delta
B^{k}-B_{k}\delta H^{k}-H_{k}\delta
B^{k}\right)g_{ij}+\frac{1}{4\pi}\left(B_{j}\delta H_{i}+H_{i}\delta
B_{j}\right),$
$\displaystyle\delta\sigma_{ij}=\check{\mu}\left(\nabla_{i}\xi_{j}+\nabla_{j}\xi_{i}-\frac{2}{3}\nabla_{k}\xi^{k}\delta_{ij}\right),$
(23)

assuming spatially uniform $\check{\mu}$. Since we are only interested in
incompressible modes, as the Hall modes are incompressible, we can ignore the
$\nabla_{i}\delta P/\rho$ term. Eq. (22) can then be written as

$\displaystyle(k^{2}\check{\mu}-\rho\omega^{2})\xi_{i}+\left(i\frac{\delta
P}{\rho}+\frac{\check{\mu}}{3}(k_{j}\xi^{j})\right)k_{i}={}$ $\displaystyle
i\Bigg{[}\frac{(k_{j}H^{j})}{4\pi}\delta
B_{i}-(k_{j}B^{j})\left\\{\left(\chi_{n}-\frac{M}{B}\right)(\hat{B}_{m}\delta
B^{m})+in_{\text{e}}\mathcal{M}_{n}(k_{j}\xi^{j})\right\\}\hat{B}_{i}$
$\displaystyle\quad-\left\\{in_{\text{e}}\mathcal{M}_{n}B(k_{j}\xi^{j})+\frac{H}{4\pi}(\hat{B}_{j}\delta
B^{j})-B\chi_{n}(\hat{B}_{j}\delta B^{j})\right\\}k_{i}\Bigg{]}.$ (24)

Using Eq. (21,24), we can find the dispersion relation for Hall MHD modes when
the electron fluid and the nuclear lattice are coupled.

Since the Hall modes are incompressible, we set $k_{j}\xi^{j}=0$. The
incompressible mode approximation breaks down in the limit that the magnetic
field pressure approaches the matter pressure, so strictly speaking the
following results will only be valid at sufficiently high densities within the
crust such that this limit is not reached. Define the Alfvén frequency
$\omega_{\text{A}}$ and shear mode frequency $\omega_{s}$ through

$\displaystyle\omega^{2}_{\text{A}}\equiv\frac{(k_{j}B^{j})(k_{j}H^{j})}{4\pi\rho},$
$\displaystyle\omega^{2}_{s}\equiv\frac{\check{\mu}k^{2}}{\rho}.$ (25)

Eq. (21) and (24) now become

$\displaystyle-i\omega\delta B^{i}$
$\displaystyle{}=\frac{c(k_{\ell}H^{\ell})}{4\pi
en_{\text{e}}}\epsilon^{ijk}k_{j}\delta
B_{k}-\frac{c(k_{j}B^{j})}{en_{\text{e}}}\left(\chi_{n}-\frac{M}{B}\right)\epsilon^{ijk}k_{j}\hat{B}_{k}(\hat{B}_{m}\delta
B^{m})+\omega(k^{\ell}B_{\ell})\xi^{i},$ (26)
$\displaystyle(\omega_{s}^{2}-\omega^{2})\xi_{i}{}$
$\displaystyle=i\frac{(k_{j}H^{j})}{4\pi}\delta
B_{i}-i(k_{j}B^{j})\left(\chi_{n}-\frac{M}{B}\right)(\hat{B}_{m}\delta
B^{m})\left(\hat{B}_{i}-\hat{k}_{i}(\hat{k}_{j}\hat{B}^{j})\right).$ (27)

since $\xi^{i}$ cannot have a component along $k^{i}$, and where
$\hat{B}^{i}-\hat{k}^{i}(\hat{k}_{j}\hat{B}^{j})$ was used to replace
$\hat{B}^{i}$ as it is the projection of $\hat{B}^{i}$ perpendicular to
$k^{i}$.

Using Eq. (27) to eliminate $\xi^{i}$ from Eq. (26) , then contracting the
result with $B^{i}$ and $\epsilon^{ijk}k_{j}\hat{B}_{k}$ gives

$\displaystyle-i\omega\left(\frac{\omega^{2}-\omega_{s}^{2}-\omega^{2}_{\text{A}}}{\omega^{2}-\omega_{s}^{2}}\right)\epsilon_{ijk}\hat{B}^{i}\hat{k}^{j}\delta
B^{k}{}$ $\displaystyle=-\frac{ck(k_{\ell}H^{\ell})}{4\pi
en_{\text{e}}}(\hat{B}^{j}\delta
B_{j})+\frac{ck(k_{\ell}B^{\ell})}{en_{\text{e}}}\left(\chi_{n}-\frac{M}{B}\right)(\hat{B}^{j}\delta
B_{j})\left(1-(\hat{k}^{j}\hat{B}_{j})^{2}\right),$ (28)
$\displaystyle-i\omega\left[\frac{\omega^{2}-\omega_{s}^{2}-\omega^{2}_{\text{A}}(1-\alpha)}{\omega^{2}-\omega_{s}^{2}}\right](\hat{B}^{j}\delta
B_{j}){}$ $\displaystyle=\frac{ck(k_{\ell}H^{\ell})}{4\pi
en_{\text{e}}}\epsilon_{ijk}\hat{B}^{i}\hat{k}^{j}\delta B^{k}.$ (29)

where we’ve defined $\hat{k}^{j}\hat{B}_{j}=\cos\theta_{B}$ and

$\alpha\equiv\frac{4\pi
B}{H}\left(\chi_{n}-\frac{M}{B}\right)\sin^{2}\theta_{B}\approx
4\pi\chi_{n}\sin^{2}\theta_{B},$ (30)

where the approximate form is valid because $H\approx B$ and $|M|/B\ll 1$.
Combining these equations we obtain the dispersion relation

$\omega^{2}=\frac{\omega_{\text{H}}^{2}(\omega^{2}-\omega_{s}^{2})^{2}(1-\alpha)}{\left(\omega^{2}-\omega_{s}^{2}-\omega^{2}_{\text{A}}\right)\left(\omega^{2}-\omega_{s}^{2}-\omega_{\text{A}}^{2}(1-\alpha)\right)}.$
(31)

In the $\omega_{s}\rightarrow\infty$ limit this recovers Eq. (12) as expected,
but in the $H\rightarrow B$ limit it does not recover Eq. (47) of Cumming et
al. (2004). In the low frequency limit
$|\omega^{2}|\ll\omega_{s}^{2},\omega_{\text{A}}^{2}$, it becomes

$\omega^{2}\approx\frac{\omega_{\text{H}}^{2}(1-\alpha)}{\left(1+\omega^{2}_{\text{A}}/\omega_{s}^{2}\right)\left(1+\omega_{\text{A}}^{2}(1-\alpha)/\omega_{s}^{2}\right)}.$
(32)

which reproduces Cumming et al. (2004) Eq. (48) (their Eq. (47) contains a
typographical error).

Eq. (32) supports an unstable Hall mode if $\alpha\approx
4\pi\chi_{n}\sin^{2}\theta_{B}>1$, which is the same source of instability as
the uncoupled case in Eq. (12). For typical magnetar field strengths and crust
densities, $\omega_{\text{A}}$ can be much greater than $\omega_{s}$, so for a
crust with finite rigidity the mode frequency is reduced and the instability
growth time is increased compared to the estimate made in Eq. (13) for an
infinitely rigid crust.

The low-frequency approximation Eq. (32) for the strong-field, crust-coupled
Hall MHD modes is plotted in Figure 1 for uniform field $B=10^{16}$ G and a
realistic crust profile. For the crust EOS we employ the BSk24 EOS (Pearson et
al., 2018). The crust used is that for a $1.4M_{\odot}$ star with radius
$R=12.59$ km, crust thickness $R_{c}=1.04$ km and neutron drip line $0.48$ km
from the surface. The significant magnetic forces that would act on the crust
have been ignored in obtaining the background model, so these results
illustrate typical numerical values associated with these Hall modes rather
than providing exact values. Panel (a) shows $\omega^{2}>0$ i.e., where the
Hall modes are oscillatory, whereas panel (b) shows $\omega^{2}<0$ i.e., where
they are unstable. The instability is possible in regions in
$B$–$n_{\text{e}}$ parameter space where new Landau levels are populated.
However, there is no instability if $\cos\theta_{B}$ is not sufficiently
small, explaining the absence of the $\cos\theta_{B}=1$ curve from panel (b).

Figure 1:  Approximate (low frequency) dispersion relation for strong-field
Hall MHD, normalized by $\omega_{\text{H}}^{2}$, as a function of position
outward from the crust-core boundary $z$ and computing using Eq. (32). Uniform
field $B=10^{16}$ G and varying angles between the wave vector $k$ and
magnetic field $B$ are shown. Note that the “kinks” not associated with Landau
quantization originate from density discontinuities at transitions from
nuclear layers. (a): Positive part. The standard crust-coupled Hall MHD result
$\omega_{\text{H}}^{2}/(1+\omega_{\text{A}}^{2}/\omega_{s}^{2})$ (Cumming et
al., 2004) is overlaid for comparison. (b): Negative of the negative part.
Unstable regions are those with nonzero negative part.

The growth time of the unstable Hall mode at a particular location within the
star will change as the field evolves. Figure 2 shows this evolution for
different locations within the model crust, plotting Eq. (13) with $\omega$
from Eq. (32) used instead of Eq. (12) in panel (a) and Eq. (13) directly in
panel (b). Note that the particular locations within the crust are
distinguished by their electron chemical potentials and not $n_{\text{e}}$,
which is a function of $B$ and changes up to a few percent as the field is
varied. Typical growth times at the specified wave number $k=R_{c}/10=10^{-4}$
cm-1 are $\sim 10^{2}$–$10^{4}$ yrs, which are comparable to magnetar
lifetimes.

Importantly, the fluid becomes stable as the strongly-quantized threshold (all
electrons in a single Landau level) is approached and $n\lesssim 10$ Landau
levels are occupied. This only occurs within the plotted range of fields for
$\mu_{\text{e}}<24.33$ MeV, corresponding to regions of the crust with
$z\gtrsim 0.45$ km. Increasing temperature (not shown) has a stabilizing
effect as in Paper I, as it smooths the sharp transitions where the
instability occurs and either increases $\tau$ or stabilizes the fluid. In
contrast, decreasing the temperature allows the instability to persist at
lower field strengths, lowering the apparent low-field cutoffs in the $z=0.1$
km and $z=0.5$ km growth times. The growth time is very sensitive to $B$, and
can change by orders of magnitude and even stabilize with a single percent
change in $B$. However stable Hall mode propagation can change the local field
value by small amounts and destabilize a region of the crust that was
previously stable.

Figure 2:  Growth time $\tau$ of the instability associated with the strong-
field Hall modes as a function of $B$. $\tau$ is computed using the low
frequency approximate dispersion relation with coupling to the lattice, Eq.
(32) (panel (a)) and without coupling to the lattice, Eq. (12) (panel (b)).
The growth times for three different locations within the crust for a single
representative wave number $k=(R_{c}/10)^{-1}=10^{-4}$ cm-1 is shown. The
curves are evaluated at $\mu_{\text{e}}=71$, $53$, $38$ MeV, correspondingly
approximately to the electron density at $z=0.1$ (black), $0.25$ (blue) and
$0.4$ km (red) respectively inside the crust. $\chi_{n}$ is evaluated at
$T=4\times 10^{7}$ K.

In the high frequency limit (the incompressible ideal MHD limit), noting that
$\omega_{\text{A}}^{2},|\omega^{2}|\gg\omega_{s}^{2},\omega^{2}_{\text{H}}$,
we find solutions

$\omega^{2}=\omega^{2}_{\text{A}},~{}\omega^{2}_{\text{A}}(1-\alpha),$ (33)

which have $\hat{B}^{i}\delta B_{i}=0$ and $\hat{B}^{i}\delta B_{i}\neq 0$
respectively. The first solution is clearly stable Alfvén modes, whereas the
second solution can be both stable and unstable depending on $B$ and
$\mu_{\text{e}}$ (or $n_{\text{e}}$). The unstable Alfvén modes are the analog
to the unstable modes discussed in Paper I; the assumption of
incompressibility means that the term proportional to the cross derivative of
$u$ with respect to $B$ and $n_{\text{e}}$ (in Paper I, the mass density
$\rho$) is absent here. These modes thus have much faster growth times than
the unstable Hall modes. Since we have ignored the displacement current term
in Ampère’s Law, the resulting Alfvén modes can be superluminal for
sufficiently large $B^{2}/\rho$ as is possible in magnetar crusts; including
the displacement field term, we instead obtain

$\omega^{2}=\frac{\omega^{2}_{\text{A}}}{1+B^{2}/(4\pi\rho
c^{2})},~{}\frac{\omega_{\text{A}}^{2}(1-\alpha)}{1+B^{2}\cos^{2}\theta_{B}/(4\pi\rho
c^{2})}.$ (34)

At sufficiently low density, $B^{2}/4\pi$ exceeds $\rho$ times the sound speed
squared $c_{s}^{2}$ and the incompressible approximation fails. Since
electrons dominate the pressure at low enough density, magnetic fields
dominate once $n_{\text{e}}\approx 1.4\times 10^{34}B_{15}^{3/2}{\rm
cm^{-3}}$. Note that approximately $\sqrt{3\pi/4e^{2}}=13.7$ Landau levels are
filled at this density.

###  3.3 Effect of finite conductivity

It is interesting to ask whether Ohmic decay can damp the unstable modes
discussed in the previous section. It can be incorporated into our analysis by
retaining the conductivity term in Eq. (6). Following analogous manipulations
to those in Section 3.2, we obtain the dispersion relation including finite
conductivity for incompressible modes:

$\omega^{2}=\frac{\omega_{\text{H}}^{2}(\omega^{2}-\omega_{s}^{2})^{2}(1-\alpha)}{\left[\omega^{2}-\omega_{s}^{2}-\omega^{2}_{\text{A}}(1-\alpha)+i\gamma_{\text{O}}(\omega^{2}-\omega_{s}^{2})\left(1-4\pi\frac{M}{B}-\frac{H}{B}\alpha\right)\right]\left[\omega^{2}-\omega_{s}^{2}-\omega^{2}_{\text{A}}+i\gamma_{\text{O}}(\omega^{2}-\omega_{s}^{2})\left(1-4\pi\frac{M}{B}\right)\right]},$
(35)

where $\gamma_{\text{O}}\equiv k^{2}/(4\pi\sigma\omega)$ and $\alpha$ is as
defined in Eq. (30). Taking $\omega=\omega_{0}-i\Gamma$ and using
$\omega^{2}\approx\omega_{0}^{2}-2i\omega_{0}\Gamma$, $|M|\ll B$, we find in
the low frequency limit $|\omega^{2}|\ll\omega_{s}^{2},\omega_{\text{A}}^{2}$

$\displaystyle\omega_{0}^{2}\approx\frac{\omega^{2}_{\text{H}}(1-\alpha)}{\left(1+\omega_{\text{A}}^{2}/\omega_{s}^{2}\right)\left(1+\omega_{\text{A}}^{2}(1-\alpha)/\omega_{s}^{2}\right)},$
$\displaystyle\Gamma\approx\frac{k^{2}}{8\pi\sigma}\frac{\left(2-\alpha+2\omega_{\text{A}}^{2}(1-\alpha)/\omega_{s}^{2}\right)}{\left(1+\omega_{\text{A}}^{2}/\omega_{s}^{2}\right)\left(1+\omega_{\text{A}}^{2}(1-\alpha)/\omega_{s}^{2}\right)}.$
(36)

The instability conditions are the same as for the case with infinite
conductivity i.e., instability if $\alpha>1$. However, the modes may be
Ohmically damped even if they are unstable provided that
$1<\alpha<(1+\omega_{\text{A}}^{2}/\omega_{s}^{2})/(1+2\omega_{\text{A}}^{2}/\omega_{s}^{2})$,
or that $1<\alpha<2$ in the perfectly rigid limit
$\omega_{s}\rightarrow\infty$. The Ohmic decay rate can even be comparable to
the instability growth rate: for example, in the perfectly rigid limit, we
have

$\frac{\Gamma}{\omega_{0}}\approx\frac{n_{\text{e}}e}{2\sigma
H\cos\theta_{B}}\frac{(2-\alpha)}{\sqrt{1-\alpha}},$ (37)

and so the ratio $n_{\text{e}}/\sigma$ determines whether Ohmic decay could
suppress the instability. The Alfvén modes are also damped: in the high
frequency limit Eq. (35) gives

$\omega^{2}=\omega_{\text{A}}^{2}(1-i\gamma_{\text{O}}),~{}\omega_{\text{A}}^{2}(1-\alpha)(1-i\gamma_{\text{O}}(1-\alpha)),$
(38)

and so again taking $\omega=\omega_{0}-i\Gamma$ we find

$\displaystyle\omega_{0}^{2}\approx\omega_{\text{A}}^{2},~{}\omega_{\text{A}}^{2}(1-\alpha),$
$\displaystyle\Gamma\approx\frac{k^{2}}{8\pi\sigma},~{}\frac{k^{2}(1-\alpha)}{8\pi\sigma}.$
(39)

###  3.4 Effect of plastic flow

The magnetic fields attainable in magnetar crusts are sufficiently strong to
exceed the yield strength of the crust, resulting in its plastic deformation.
The crust can be modeled (e.g. Beloborodov & Levin (2014)) as an elasto-
viscoplastic solid (Irgens, 2008). Below the yield criterion, the crust
behaves as an elastic material with shear stress tensor given by Eq. (15). The
yield stress $\sigma_{y}$ of a neutron star crust has been estimated through
numerical simulations to be

$\sigma_{y}\approx
0.03143e^{2}Z^{2/3}n_{\text{e}}^{4/3}\approx\check{\mu}E_{y}.$ (40)

Here $E_{y}\sim 0.1$ is the yield strain (Horowitz & Kadau, 2009; Chugunov &
Horowitz, 2010) and we have ignored thermal corrections. The exact value of
$E_{y}$ is not well known, and may be lower than the value used here (e.g.,
Baiko & Chugunov (2018)). The transition to the plastic regime begins when
$\sigma_{\text{M}}>\sigma_{y}$ where $\sigma_{\text{M}}$ is the von Mises
stress

$\sigma_{\text{M}}=\sqrt{\frac{3}{2}\sigma_{ij}\sigma^{ij}}.$ (41)

Beyond this point, the shear stress is determined using the Bingham–Maxwell
model

$\frac{1}{2\eta}\left[1-\frac{\sigma_{y}}{\sigma_{\text{M}}}\right]\theta(\sigma_{\text{M}}-\sigma_{y})\sigma_{ij}+\frac{1}{2\check{\mu}}\partial_{t}\sigma_{ij}=\partial_{t}E^{s}_{ij},$
(42)

where $\theta(x)$ is the step function and $\eta$ is the effective dynamic
viscosity in the plastic flow regime. $\eta$ is poorly understood
theoretically for neutron star crusts: estimates for its value based on
studies of magnetar magneto-plastic field evolution are of order
$10^{36}$–$10^{38}$ g cm-1 s-1 (Lander, 2016; Lander & Gourgouliatos, 2019). A
theoretical calculation of its (temperature-dependent) value found $\eta\sim
10^{31}$ g cm-1 s-1 (Kwang-Hua, 2018), comparable to the effective value of
$\approx 7\times 10^{31}$ g cm-1 s-1 employed by Li et al. (2016), though the
previous references argue that this value is too low to support the persistent
twisting of magnetic field lines anchored to the crust which is a leading
mechanism for generating magnetar outbursts.

Assuming harmonic time dependence $\exp(-i\omega t)$ for the displacement, Eq.
(42) implies a linear relation between $\sigma_{ij}$ and $E^{s}_{ij}$ of

$\sigma_{ij}=2\left[\frac{\check{\mu}}{1+i\check{\mu}/(\eta\omega)\left[1-\sigma_{y}/\sigma_{\text{M}}\right]\theta(\sigma_{\text{M}}-\sigma_{y})}\right]\equiv
2\tilde{\mu}(\omega)E^{s}_{ij},$ (43)

so to determine the effect of crust plasticity on the Hall mode dispersion
relation amounts to replacing $\check{\mu}$ in $\omega_{s}^{2}$ with the
dynamic modulus $\tilde{\mu}(\omega)$. Under the approximation that the
pressure gradient and gravity will balance each other, the von Mises stress in
the crust can be estimated by equating $\sigma_{ij}=-T^{B}_{ij}$, hence

$\sigma_{\text{M}}\approx\sqrt{\frac{9}{128\pi^{2}}}B^{2}\approx\frac{B^{2}}{4\pi},$
(44)

so the crust will be in the plastic flow regime for densities below the
threshold

$\rho_{\text{crit}}\approx\frac{m_{\text{n}}n_{\text{e}}}{Y}\approx 5\times
10^{12}\left(\frac{40}{Z}\right)^{1/2}\left(\frac{0.5}{Y}\right)\left(\frac{B}{10^{15}\text{
G}}\right)^{3/2}\left(\frac{0.1}{E_{y}}\right)^{3/4}\text{ g cm}^{-3},$ (45)

for neutron mass $m_{\text{n}}$. We note that $\rho_{\text{crit}}$ is greater
than the density at which $B^{2}/(8\pi)\approx P_{\text{e}}$. For
$\sigma_{\text{M}}\geq\sigma_{y}$, there are two limiting cases. Define
$\overline{\eta}=\eta(1-\sigma_{y}/\sigma_{\text{M}})$ and an effective
relaxation time $\tau=\overline{\eta}/\check{\mu}$ (Landau & Lifshitz, 1970).
In the solid limit $|\tau\omega|\gg 1$ we have
$\tilde{\mu}(\omega)\approx\check{\mu}(1-i\check{\mu}/(\eta\omega))$, while in
the fluid limit $|\tau\omega|\ll 1$ and
$\tilde{\mu}(\omega)\approx-i\overline{\eta}\omega$.

For the $|\tau\omega|\gg 1$ case, the modification of the oscillation modes in
the low frequency limit due to crustal plasticity can be estimated using the
Newton–Raphson method with Eq. (32) as the initial guess, giving

$\omega\approx\frac{\omega_{\text{H}}\sqrt{1-\alpha}}{\sqrt{\left(1+\omega^{2}_{\text{A}}/\omega_{s}^{2}\right)\left(1+\omega_{\text{A}}^{2}(1-\alpha)/\omega_{s}^{2}\right)}}-i\frac{\check{\mu}\left[2\left(1+\omega_{\text{A}}^{2}/\omega_{s}^{2}\right)-\alpha\left(1+2\omega_{\text{A}}^{2}/\omega_{s}^{2}\right)\right]}{2\overline{\eta}\left(1+\omega_{\text{A}}^{2}(1-\alpha)/\omega_{s}^{2}\right)\left(1+\omega_{\text{A}}^{2}/\omega_{s}^{2}\right)},$
(46)

where $\omega_{s}^{2}$ is defined as in Eq. (25) in terms of $\check{\mu}$.
This retains the instability for $\alpha>1$, but it can be completely damped
out by the plastic viscosity, with lower $\check{\mu}/\overline{\eta}$
providing greater damping. If $\omega_{\text{A}}^{2}/\omega_{s}^{2}\approx
0.1\sigma_{\text{M}}/\sigma_{y}\gg 1$, which is the elasto-viscoplastic
regime, then this approximates to

$\omega\approx\frac{\omega_{\text{H}}\omega_{s}^{2}}{\omega_{\text{A}}^{2}}-i\frac{\check{\mu}}{\overline{\eta}},$
(47)

where the instability vanishes and the modes are only damped.

In the $|\tau\omega|\ll 1$ case, we replace $\omega_{s}^{2}\rightarrow-i\omega
k^{2}\overline{\eta}/\rho=-i\omega k^{2}\nu$ in Eq. (31). In the high
frequency limit $|\omega|^{2}\gg\omega_{\text{H}}^{2}$, this gives two modes
with dispersions

$\displaystyle\omega_{1}=-i\frac{\nu
k^{2}}{2}\pm\omega_{\text{A}}\sqrt{1-\frac{\nu^{2}k^{4}}{4\omega^{2}_{\text{A}}}}\approx-i\frac{\nu
k^{2}}{2}\pm\omega_{\text{A}},$ $\displaystyle\omega_{2}=-i\frac{\nu
k^{2}}{2}\pm\omega_{\text{A}}\sqrt{1-\alpha-\frac{\nu^{2}k^{4}}{4\omega^{2}_{\text{A}}}}\approx-i\frac{\nu
k^{2}}{2}\pm\omega_{\text{A}}\sqrt{1-\alpha}.$ (48)

Since $\nu k^{2}\ll\check{\mu}k^{2}/\rho$ by assumption and
$\check{\mu}k^{2}/\rho\ll\omega_{\text{A}}^{2}$ in the elasto-viscoplastic
regime, the imaginary damping part of these modes is much smaller in magnitude
than the real part and the modes are only slightly affected by the plastic
viscosity. In the low frequency limit
$|\omega|^{2},\omega_{\text{H}}^{2}\ll\omega_{\text{A}}^{2}$, we obtain

$\omega\approx-i\frac{\nu
k^{2}\omega_{\text{A}}^{2}(2-\alpha)}{2(\omega_{\text{H}}^{2}(1-\alpha)+\nu^{2}k^{4})}\pm\frac{\omega_{\text{A}}^{2}\omega_{\text{H}}(1-\alpha)}{\omega_{\text{H}}^{2}(1-\alpha)+\nu^{2}k^{4}}\sqrt{1-\frac{\alpha^{2}\nu^{2}k^{4}}{4\omega_{\text{H}}^{2}(1-\alpha)^{2}}}.$
(49)

This can be simplified if $\omega_{\text{H}}\ll\nu k^{2}$ to give

$\omega\approx-i\frac{\nu
k^{2}(2-\alpha)}{2(1-\alpha)}\frac{\omega_{\text{A}}^{2}}{\omega_{\text{H}}^{2}}\pm\frac{\omega_{\text{A}}^{2}}{\omega_{\text{H}}},$
(50)

which is strongly damped since $\omega_{\text{A}}\gg\omega_{\text{H}}$. But
these are not the usual Hall modes, which therefore do not exist in the
$|\tau\omega|\ll 1$ elasto-viscoplastic regime.

##  4 Nonzero temperatures and the instability parameter space

It was shown in Paper I that finite temperatures can suppress the instability
associated with Landau quantized fermions in ideal MHD using the full
Fermi–Dirac integrals for the thermodynamic quantities relevant for the
instabilities i.e. the analogs to $M$, $\chi_{n}$ and $\mathcal{M}_{n}$ as
defined in this paper. In this section we first describe approximations for
the finite temperature corrections to $M$ and $\chi_{n}$, then use these
approximations to determine the critical temperature for onset of the strong-
field Hall MHD instability as discussed in the previous section. In doing so,
we also derive an approximate criterion for the maximum number of occupied
Landau levels at zero temperature $n_{\text{max}}$ required for the
instability. This allows us to constrain the regions in $B$–$n_{\text{e}}$–$T$
parameter space in which the instability could be active. In this section we
work in units $\hbar=c=k_{B}=1$.

###  4.1 Low temperature corrections to thermodynamic quantities

The standard method to determine the temperature-dependence of a Fermi gas,
for low to moderate temperatures, is to employ the Sommerfeld expansion (e.g.,
Kittel (2005); Ashcroft & Mermin (1976)). However, this expansion breaks down
for the thermodynamic potentials for Landau-quantized fermions. Consider $M$:
the finite temperature integral is given by

$M=\left(\frac{\partial P_{\text{e}}}{\partial
B}\right)_{\mu_{\text{e}},T}=-\frac{e}{2\pi^{2}}\sum_{n=0}^{\infty}g_{n}\int_{m_{n}}^{\infty}\frac{dE(E^{2}-m_{n}^{2}-neB)}{\sqrt{E^{2}-m_{n}^{2}}\left\\{\exp[\beta(E-\mu_{\text{e}})]+1\right\\}},$
(51)

where $m_{n}\equiv\sqrt{m_{\text{e}}^{2}+2eBn}$, $\beta=T^{-1}$ and
$g_{n}=2-\delta_{n,0}$. This can be expanded at finite temperature as

$M(T)=M(T=0)+\frac{e\mu_{\text{e}}T^{2}}{12}\sum_{n=0}^{n_{\text{max}}-1}g_{n}\frac{\mu_{\text{e}}^{2}-m_{n}^{2}+eBn}{(\mu_{\text{e}}^{2}-m_{n}^{2})^{3/2}}+M_{n\geq
n_{\text{max}}}(T\neq 0),$ (52)

where the terms on the right are the $T=0$ solution, the standard lowest-order
Sommerfeld expansion for $n<n_{\text{max}}$, and additional corrections for
$n\geq n_{\text{max}}$. The Sommerfeld expansion is valid for
$\mu_{\text{e}}-m_{n}\gg T$, which is possible for $n<n_{\text{max}}$, but it
fails for $|\mu_{\text{e}}-m_{n}|\lesssim T$ or $m_{n}-\mu_{\text{e}}\gg T$
which occurs for $n\geq n_{\text{max}}$. Moreover, while $M(0)$ is finite, its
derivatives with respect to $\mu_{\text{e}}$ and $B$ are not, with divergences
arising from the $n=n_{\text{max}}$ term. Hence a different approximation must
be used to include finite temperature corrections for $n\gtrsim
n_{\text{max}}$ i.e., for small $|\mu_{\text{e}}-m_{n}|/T$.

After a change of integration variable to $y=\sqrt{\beta(E-m_{n})}$, Eq. (51)
can be rewritten in the form

$M=\frac{e\sqrt{T}}{2\pi^{2}}\sum_{n=0}^{\infty}g_{n}\sqrt{2m_{n}}\left[2T\int_{0}^{\infty}\frac{dyy^{2}\sqrt{1+y^{2}T/(2m_{n})}}{\exp(y^{2}-Y_{n})+1}-\frac{neB}{{m_{n}}}\int_{0}^{\infty}\frac{dy}{\sqrt{1+y^{2}T/(2m_{n})}[\exp(y^{2}-Y_{n})+1]}\right],$
(53)

where $Y_{n}\equiv\beta(\mu_{\text{e}}-m_{n})$. For low temperatures
$\mu_{\text{e}}\gg T$ and $n<n_{\text{max}}$, $Y_{n}\gg 0$. This is not
necessarily true for $n\geq n_{\text{max}}$. For $n=n_{\text{max}}$, the
contributions to both integrals for low temperatures are exponentially
suppressed for $y^{2}>Y_{n_{\text{max}}}$, and since $\mu_{\text{e}}\sim
m_{n_{\text{max}}}$, $Y_{n_{\text{max}}}T/(2m_{n_{\text{max}}})\ll 1$ when
$y^{2}\lesssim Y_{n_{\text{max}}}$. So Taylor expanding to zeroth order the
integrands in Eq. (53) about $Y_{n}T/(2m_{n})\ll 1$ gives the temperature-
dependent contribution to the sum over Landau levels from $n=n_{\text{max}}$

$\left(M\right)_{n=n_{\text{max}}}\approx\frac{e}{\pi^{2}}\sqrt{\frac{m_{n_{\text{max}}}T}{2}}g_{n_{\text{max}}}\left[2T\int_{0}^{\infty}\frac{dyy^{2}}{\exp(y^{2}-Y_{n_{\text{max}}})+1}-\frac{n_{\text{max}}eB}{m_{n_{\text{max}}}}\int_{0}^{\infty}\frac{dy}{\exp(y^{2}-Y_{n_{\text{max}}})+1}\right].$
(54)

Since both integrals are $\sim 1$ for small value of $Y_{n}$, in this limit
the second term dominates if $T\ll n_{\text{max}}eB/(2m_{n_{\text{max}}})$.
Note that the two integrals only depend on a single parameter
$Y_{n_{\text{max}}}$. In general, this approximation could be made more
accurate by including higher-order terms in the Taylor expansions of
$(1+y^{2}T/(2m_{n}))^{\pm 1/2}$.

Recall that $n_{\text{max}}$ is the integer part of $p_{F}^{2}/2eB$ for Fermi
momentum $p_{F}=\sqrt{\mu^{2}_{\text{e}}-m^{2}_{\text{e}}}$; therefore
$\mu_{\text{e}}^{2}=p_{F}^{2}+m_{\text{e}}^{2}>2eBn_{\text{max}}+m_{\text{e}}^{2}=m_{n_{\text{max}}}^{2}$.
At $T=0$ there are no terms in the sums with $n>n_{\text{max}}$, but at
nonzero $T$ these terms are present. The upshot is that there are two cases to
worry about.

Now let us consider what happens when $\mu_{\text{e}}$ is near but just above
$m_{n_{\text{max}}}$. In this first case, all terms with $n<n_{\text{max}}$
can be approximated by their $T=0$ values, assuming that

$\mu_{\text{e}}-m_{n_{\text{max}}-1}\simeq
m_{n_{\text{max}}}-m_{n_{\text{max}}-1}=m_{n_{\text{max}}}\sqrt{m_{n_{\text{max}}}^{2}-2eB}\simeq\frac{eB}{m_{n_{\text{max}}}}\simeq\frac{eB}{\mu_{\text{e}}}\gg
T,$ (55)

which is equivalent to the requirement that the Landau level spacing at the
Fermi surface is larger than the temperature. We are therefore left with
considering just the terms $n\geq n_{\text{max}}$. The dominant terms are

$-\frac{ne^{2}B}{\pi^{2}}\sqrt{\frac{2T}{m_{n}}}\int_{0}^{\infty}\frac{\text{d}y}{\sqrt{(1+y^{2}T/2m_{n})}[\exp(y^{2}-Y_{n})+1]}\simeq-\frac{ne^{2}B}{\pi^{2}}\sqrt{\frac{2T}{m_{n}}}\int_{0}^{\infty}\frac{\text{d}y}{\exp(y^{2}-Y_{n})+1}\equiv-\frac{ne^{2}BI_{0}(Y_{n})}{\pi^{2}}\sqrt{\frac{2T}{m_{n}}}~{}.$
(56)

Since $Y_{n_{\text{max}}}>0$ in this case, the integrand is largest near
$y=\sqrt{Y_{n_{\text{max}}}}$, and cuts off exponentially within $\delta y\sim
1$ of its peak value. As a result $I_{0}(Y_{n_{\text{max}}})\sim 1$ for small
values of $Y_{n_{\text{max}}}$; the contribution from $n=n_{\text{max}}$ to
the magnetization is of order $(ep_{F}^{2}/\pi^{2})\sqrt{T/\mu_{\text{e}}}$,
which goes to zero as $T\to 0$. The $T=0$ limit is approached once
$Y_{n_{\text{max}}}$ is larger than $\sim$ a few i.e., for
$\mu_{\text{e}}-m_{n_{\text{max}}}\gtrsim({\rm a\,few})\times T$, so
$\mu_{\text{e}}$ is still very near $m_{n_{\text{max}}}\gg T$. The
contribution from $n=n_{\text{max}}+1$ is suppressed because
$Y_{n_{\text{max}}+1}=(\mu_{\text{e}}-m_{n_{\text{max}}+1})/T\simeq(m_{n_{\text{max}}}-m_{n_{\text{max}}+1})/T\simeq-
eB/\mu_{\text{e}}T$, which is very negative for $T$ small compared to the
Landau level spacing, and therefore
$I_{0}(Y_{n_{\text{max}}+1}\sim\exp(-|Y_{n_{\text{max}}+1}|)\ll 1$.

Next consider what happens when $\mu$ is near but just below
$m_{n_{\text{max}}}$. In this case, the contributions from $n<n_{\text{max}}$
are once again close to their $T=0$ values. The contribution from
$\nu=n_{\text{max}}$ can still be found from Eq. (56) but now
$Y_{n_{\text{max}}}<0$. Although this contribution is exponentially small for
large $|Y_{n_{\text{max}}}|$, for
$0>\mu_{\text{e}}-m_{n_{\text{max}}}\gtrsim-T$ it is of the same order as the
contribution for $0<\mu_{\text{e}}-m_{n_{\text{max}}}\lesssim T$. Thus, the
magnetization is approximately

$M\simeq M_{0\leq
n<n_{\text{max}}}(B,\mu_{\text{e}},0)-\frac{ep_{F}^{2}I_{0}(Y_{n_{\text{max}}})}{2\pi^{2}}\sqrt{\frac{2T}{\mu_{\text{e}}}}~{}$
(57)

where the first term on the right-hand side is the contribution from all
Landau levels below $n=n_{\text{max}}$ calculated at $T=0$, and where
$Y_{n_{\text{max}}}$ can be positive or negative. Completely analogous
reasoning for the electron density

$n_{\text{e}}=\frac{eB}{2\pi^{2}}\sum_{n=0}^{\infty}g_{n}\int_{m_{n}}^{\infty}\frac{\text{d}E\,E}{\sqrt{E^{2}-m_{n}^{2}}\,\\{\exp[(E-\mu_{\text{e}})/T]+1\\}}=\frac{eB\sqrt{T}}{\pi^{2}\sqrt{2}}\sum_{n=0}^{\infty}{g_{n}\sqrt{m_{n}}}\int_{0}^{\infty}\frac{dy(1+y^{2}T/m_{n})}{\sqrt{1+y^{2}T/(2m_{n})}[\exp(y^{2}-Y_{n})+1]}$
(58)

leads to the approximation

$n_{\text{e}}\simeq n_{\text{e},0\leq
n<n_{\text{max}}}(B,\mu_{\text{e}},0)+\frac{eB\sqrt{2T\mu_{\text{e}}}I_{0}(Y_{n_{\text{max}}})}{\pi^{2}}.$
(59)

In the $n=n_{\text{max}}$ terms in Eqs. (57) and (59) only the leading terms
in $T$ have been retained and $\mu_{\text{e}}$ has been substituted for
$m_{n_{\text{max}}}$ since the most temperature-sensitive dependence in the
derivatives $\partial(M,n_{e})/\partial B$ will arise from differentiating
$I_{0}(Y)$.

Since the instability criterion as found in Eq. (12) depends on $\chi_{n}$, we
also want to determine the lowest-order temperature-dependent corrections to
this quantity. From Eq. (4) we have

$\left(\frac{\partial H}{\partial
B}\right)_{n_{\text{e}},T}=\left(\frac{\partial H}{\partial
B}\right)_{\mu_{\text{e}},T}-\left(\frac{\partial
H}{\partial\mu_{\text{e}}}\right)_{B,T}\left(\frac{\partial
n_{\text{e}}}{\partial\mu_{\text{e}}}\right)_{B,T}^{-1}\left(\frac{\partial
n_{\text{e}}}{\partial B}\right)_{\mu_{\text{e}},T}.$ (60)

When $\beta(\mu_{\text{e}}-m_{n})\lesssim 1$ i.e., when
$(\mu_{\text{e}}^{2}-m_{\text{e}}^{2})/(2eB)$ is nearly an integer and where
the instability occurs, we can drop the contribution from the
$n=n_{\text{max}}+1$ Landau level in Eq. (57–59). Here we hold $n_{\text{e}}$
fixed because it, as opposed to $\mu_{\text{e}}$, is a natural variable for
the internal energy density using which the equations of motion were derived.
Taking the required partial derivatives for computing Eq. (60) (following
e.g., Paper I Appendix B (online supplement) for the required partial
derivatives of $H$ at $T=0$) gives

$\displaystyle\left(\frac{\partial H(B,\mu_{\text{e}},T)}{\partial
B}\right)_{\mu_{\text{e}},T}$ $\displaystyle\simeq\left(\frac{\partial
H_{0\leq n<n_{\text{max}}}(B,\mu_{\text{e}},0)}{\partial
B}\right)_{\mu_{\text{e}},T}-\frac{2e^{2}n_{\text{max}}}{\pi}\sqrt{\frac{2\mu_{\text{e}}}{T}}I^{\prime}(Y_{n_{\text{max}}})$
$\displaystyle\simeq
1-\frac{2e^{2}n_{\text{max}}}{\pi}\Bigg{[}2\left\\{v_{F}\Sigma_{1/2}-(1-v_{F}^{2})\Sigma_{S}-2v_{F}^{2}\Sigma_{S,n}\right\\}-\left(\frac{2+v_{F}^{2}}{v_{F}}\right)\Sigma_{1/2}+v_{F}\Sigma_{-1/2}-v_{F}(1-v_{F}^{2})\Sigma_{X}$
$\displaystyle\qquad\qquad\qquad\qquad+2\left(\frac{1}{v_{F}^{2}}-1\right)\Sigma_{S}+\sqrt{\frac{2\mu_{\text{e}}}{T}}I^{\prime}(Y_{n_{\text{max}}})\Bigg{]},$
(61a) $\displaystyle\left(\frac{\partial
H(B,\mu_{\text{e}},T)}{\partial\mu_{\text{e}}}\right)_{B,T}$
$\displaystyle\simeq\left(\frac{\partial H_{0\leq
n<n_{\text{max}}}(B,\mu_{\text{e}},0)}{\partial\mu_{\text{e}}}\right)_{B,T}+\frac{2e\mu_{\text{e}}}{\pi}\sqrt{\frac{2\mu_{\text{e}}}{T}}I^{\prime}(Y_{n_{\text{max}}})=-\frac{2e\mu_{\text{e}}}{\pi}\left[3\Sigma_{1/2}-\Sigma_{-1/2}-\sqrt{\frac{2\mu_{\text{e}}}{T}}I^{\prime}(Y_{n_{\text{max}}})\right],$
(61b) $\displaystyle=-4\pi\left(\frac{\partial
n_{\text{e}}(B,\mu_{\text{e}},T)}{\partial B}\right)_{\mu_{\text{e}},T}$
$\displaystyle\left(\frac{\partial
n_{\text{e}}(B,\mu_{\text{e}},T)}{\partial\mu_{\text{e}}}\right)_{B,T}$
$\displaystyle\simeq\left(\frac{\partial n_{\text{e},0\leq
n<n_{\text{max}}}(B,\mu_{\text{e}},0)}{\partial\mu_{\text{e}}}\right)_{B,T}+\frac{\mu_{\text{e}}^{2}}{2\pi^{2}n_{\text{max}}}\sqrt{\frac{2\mu_{\text{e}}}{T}}I^{\prime}(Y_{n_{\text{max}}})=\frac{\mu_{\text{e}}^{2}}{2\pi^{2}n_{\text{max}}}\left[\Sigma_{-1/2}+\sqrt{\frac{2\mu_{\text{e}}}{T}}I^{\prime}(Y_{n_{\text{max}}})\right],$
(61c)

where $v_{F}^{2}\equiv(1-m_{\text{e}}^{2}/\mu_{\text{e}}^{2})$ and we use
$\mu_{\text{e}}^{2}\approx 2eBn_{\text{max}}$ so

$\displaystyle\Sigma_{1/2}={}$
$\displaystyle\sum_{n=0}^{n_{\text{max}}-1}\frac{g_{n}}{2}\sqrt{1-\frac{n}{n_{\text{max}}}},$
$\displaystyle\Sigma_{S}={}$
$\displaystyle\sum_{n=0}^{n_{\text{max}}-1}\frac{g_{n}}{2}\text{arsinh}\sqrt{\frac{v_{F}^{2}(n_{\text{max}}/n-1)}{v_{F}^{2}+(1-v_{F}^{2})n_{\text{max}}/n}}$
$\displaystyle\Sigma_{S,n}={}$
$\displaystyle\sum_{n=1}^{n_{\text{max}}-1}\frac{n}{n_{\text{max}}}\text{arsinh}\sqrt{\frac{v_{F}^{2}(n_{\text{max}}/n-1)}{v_{F}^{2}+(1-v_{F}^{2})n_{\text{max}}/n}},$
$\displaystyle\Sigma_{-1/2}={}$
$\displaystyle\sum_{n=0}^{n_{\text{max}}-1}\frac{g_{n}}{2\sqrt{1-n/n_{\text{max}}}},$
(62a) $\displaystyle\Sigma_{X}={}$
$\displaystyle\sum_{n=0}^{n_{\text{max}}-1}\frac{g_{n}}{2}\frac{n\sqrt{1-n/n_{\text{max}}}}{n_{\text{max}}(1+v_{F}^{2}(n/n_{\text{max}}-1))},$
$\displaystyle\Sigma_{-3/2}={}$
$\displaystyle\sum_{n=0}^{n_{\text{max}}-1}\frac{g_{n}}{2(1-n/n_{\text{max}})^{3/2}}.$

$I^{\prime}(Y)$ is

$I^{\prime}(Y)=\frac{\text{d}I(Y)}{\text{d}Y}=\int^{\infty}_{0}\frac{\text{d}y\exp(y^{2}-Y)}{[\exp(y^{2}-Y)+1]^{2}},$
(63)

and is plotted in Figure 3 panel (a). It has asymptotic forms
$I^{\prime}(Y)\sim 1/(2\sqrt{Y})$ for $Y\gg 1$ and $I^{\prime}(Y)\sim\exp(Y)$
for $Y\ll 0$, and a peak value of $I^{\prime}(Y\approx 1.1)\approx 0.3955$.

(a)  (b)

Figure 3: (a) $I^{\prime}(Y)=\text{d}I/\text{d}Y$ as a function of $Y$ (solid
line). The asymptotic form $I^{\prime}(Y)\sim 1/(2\sqrt{Y})$ for $Y\gg 1$ is
superposed (dashed line). (b) Scaled asymptotic forms of different sums as
defined in Eq. (62a) as a function of $n_{\text{max}}$, justifying the
scalings given in Eq. (65). $v_{F}^{2}=1$ is assumed in computing
$\Sigma_{S,n}$.

Combining Eq. (61a–61c), working in the ultrarelativistic regime $v_{F}\approx
1$, and only retaining terms $\mathcal{O}(T^{-1/2})$ and $\mathcal{O}(T^{0})$
gives the desired approximation for $(\partial H/\partial B)_{n_{\text{e}}}$

$\left(\frac{\partial H}{\partial B}\right)_{n_{\text{e}}}\simeq
1-\frac{4e^{2}n_{\text{max}}}{\pi}\left(\Sigma_{1/2}-2\Sigma_{S,n}\right)+\frac{6e^{2}n_{\text{max}}\Sigma_{1/2}}{\pi}\left[-1+\frac{3\Sigma_{1/2}}{\Sigma_{-1/2}+\sqrt{2\mu_{\text{e}}/T}I^{\prime}(Y_{n_{\text{max}}})}\right].$
(64)

###  4.2 Unstable region of $B$–$n_{\text{e}}$–$T$ parameter space

Since the instability criterion is roughly $1-4\pi\chi_{n}<0$ or equivalently
$\partial H/\partial B|_{n_{\text{e}}}<0$, we can use Eq. (64) to estimate the
region of parameter space where the strong field Hall MHD instability may be
active. We can first estimate the values of the sums in Eq. (62a). These
expressions are shown in Figure 3 panel (b) as a function of $n_{\text{max}}$,
which shows that they have the following asymptotic forms for large
$n_{\text{max}}$ in the $m_{\text{e}}\rightarrow 0$ limit:

$\displaystyle\Sigma_{1/2}\simeq 2n_{\text{max}}/3,$
$\displaystyle\Sigma_{-1/2}\simeq 2n_{\text{max}},$ $\displaystyle
3\Sigma_{1/2}-\Sigma_{-1/2}\simeq 1.5\sqrt{n_{\text{max}}},$
$\displaystyle\Sigma_{-3/2}\simeq 2.6n_{\text{max}}^{3/2},$
$\displaystyle\Sigma_{1/2}-2\Sigma_{S,n}\sim 1/n_{\text{max}}.$ (65)

In the zero-temperature limit and working to leading order in
$n_{\text{max}}$, Eq. (64) reduces to

$\left(\frac{\partial H}{\partial B}\right)_{n_{\text{e}}}\simeq
1-\frac{4e^{2}n_{\text{max}}}{\pi}\left(\Sigma_{1/2}-2\Sigma_{S,n}\right)-\frac{6e^{2}n_{\text{max}}\Sigma_{1/2}}{\pi}\approx
1-\frac{4e^{2}n_{\text{max}}^{2}}{\pi}.$ (66)

This condition implies that $n_{\text{max}}\gtrsim\sqrt{\pi/(4e^{2})}=10.37$
is necessary for the instability (assuming $n_{\text{e}}$ is held fixed in the
derivatives), and since $n_{\text{max}}$ is determined by $B$ and
$\mu_{\text{e}}$, it sets the possible region in parameter space where the
instability could be active under the given approximations. Figure 4 shows the
potentially unstable region in $B$–$\mu_{\text{e}}$ and $B$–$n_{\text{e}}$
parameter spaces at $T=0$ and in the large $n_{\text{max}}$ approximation.
Since $n_{\text{e}}$ is a function of $B$, in the latter we use the
$n_{\text{e}}(B=0)$ value as an independent variable.

(a)  (b)

Figure 4: The contour of constant $n_{\text{max}}=\sqrt{\pi/(4e^{2})}$
(approximating $n_{\text{max}}$ as a real number and not an integer) as a
function of $B$ and $\mu_{\text{e}}$ (a) and $B$ and $n_{\text{e}}(T=0)$ (b),
with the potentially stable and unstable regions of parameter space split by
these contours labelled.

Including the lowest-order finite temperature corrections, we find that the
condition for instability is

$\frac{3\Sigma_{1/2}}{\Sigma_{-1/2}+\sqrt{2\mu_{\text{e}}/T}I^{\prime}(Y_{n_{\text{max}}})}<1-\frac{\pi}{6e^{2}n_{\text{max}}\Sigma_{1/2}}-\frac{2(2\Sigma_{S,n}-\Sigma_{1/2})}{3\Sigma_{1/2}}=\frac{5}{3}-\frac{4\Sigma_{S,n}}{3\Sigma_{1/2}}-\frac{\pi}{6e^{2}n_{\text{max}}\Sigma_{1/2}}.$
(67)

For values of $n_{\text{max}}\lesssim 10$, the right-hand side of this
equation cannot be positive, and since left-hand side of the equation is
always positive the instability is not active for such low values of
$n_{\text{max}}$. Assuming that $n_{\text{max}}$ is large enough to make the
right-hand side positive, and noting that $I^{\prime}(Y)\lesssim 0.3955$, the
necessary condition for stability is

$\sqrt{\frac{2\mu_{\text{e}}}{T}}>2.528\left(\frac{3\Sigma_{1/2}}{5/3-4\Sigma_{S,n}/(3\Sigma_{1/2})-\pi/(6e^{2}n_{\text{max}}\Sigma_{1/2})}-\Sigma_{-1/2}\right)\equiv\beta_{\text{crit}}.$
(68)

and hence a critical temperature above which the instability is inactive
$T_{\text{crit}}=2\mu_{\text{e}}/\beta_{\text{crit}}^{2}$. For large
$n_{\text{max}}$ Eq. (68) reduces to $\beta_{\text{crit}}\approx
2.528(3\Sigma_{1/2}-\Sigma_{-1/2})\approx 3.8\sqrt{n_{\text{max}}}$. This
implies a critical temperature of

$T_{\text{crit}}\approx 1.6\times 10^{8}\left(\frac{\mu_{\text{e}}}{10\text{
MeV}}\right)\left(\frac{n_{\text{max}}}{100}\right)^{-1}\text{ K}.$ (69)

Figure 5 panel (a) shows Eq. (64) as a function of $Y$ for fixed
$n_{\text{max}}$ and different values of $\sqrt{2\mu_{\text{e}}/T}$, labelled
by their scaling compared to
$\beta_{\text{crit}}=\sqrt{2\mu_{\text{e}}/T_{\text{crit}}}$ as defined by Eq.
(68). It demonstrates the suppression of the instability as the temperature is
increased, with the range of unstable values of $Y$ decreasing to a single
point at $T_{\text{crit}}$.

(a)  (b)

Figure 5: (a) Eq. (64) for fixed $n_{\text{max}}=50$ as a function of $Y$.
Points above $(\partial H/\partial B)_{n_{\text{e}}}=0$ are stable.
$\beta_{\text{crit}}=37.00384$, defined by Eq. (68), is the value marking the
transition between stability and instability, and Eq. (64) is plotted for this
value and for different multiples of it. $v_{F}^{2}=1$ is assumed. (b)
$T_{\text{crit}}/\mu_{\text{e}}$ as a function of $n_{\text{max}}$ computed
using Eq. (68) for $n_{\text{e}}$ fixed and Eq. (61a) for $\mu_{\text{e}}$
fixed. $v_{F}^{2}=1$ is assumed for both curves. Note the critical temperature
is zero for $n_{\text{max}}\leq 10$ only in the fixed $n_{\text{e}}$ case.

The contours of Figure 6 show the critical temperature for the instability in
$B$–$\mu_{\text{e}}$ and $B$–$n_{\text{e}}(T=0)$ parameter space using the
large $n_{\text{max}}$ approximation Eq. (69). The critical temperature
increases as the density decreases and the magnetic field increases i.e. as
$n_{\text{max}}$ increases. But the approximation is not valid in the upper
left-hand regions of each plot.

Figure 7 shows the contours of fixed critical temperature for the instability
using the full form of $\beta_{\text{crit}}$ in Eq. (68). The critical
temperature generally increases as the density decreases and the magnetic
field increases i.e., as $n_{\text{max}}$ increases. Moving from high
$\mu_{\text{e}}$ to low $\mu_{\text{e}}$ along nearly constant $B$, the
contours eventually turn around and follow curves of nearly constant
$\mu_{\text{e}}$ for increasing $B$. This behaviour can be understood by
plotting $T_{\text{crit}}$ as a function of $n_{\text{max}}$, as is done in
Figure 5 (b). This shows that the critical temperature is zero for
$n_{\text{max}}\leq 10$, increases to a peak value for $n_{\text{max}}=48$,
where $T_{\text{crit}}/\mu_{\text{e}}\approx 0.00146$, and then decreases for
larger $n_{\text{max}}$. This explains the “turn-around” in the critical
temperature contours: on a curve of fixed $\mu_{\text{e}}$ in the left panel
of Figure 7, as $B$ is increased, the curve can intersect the same constant
$T_{\text{crit}}$ contour twice. As $B$ is increased, $n_{\text{max}}$
decreases, and so going from right to left across Figure 5 panel (b),
$T_{\text{crit}}$ will have the same value for two different values of
$n_{\text{max}}$ or two different values of $B$ for fixed $\mu_{\text{e}}$.
$T_{\text{crit}}$ at fixed $\mu_{\text{e}}$ is also plotted for comparison:
this is found by using Eq. (61a) while setting $\partial H/\partial
B|_{\mu_{\text{e}},T}=0$ and taking $I^{\prime}(Y)=0.3955$. Note that there is
no critical value of $n_{\text{max}}$ below which the instability is absent,
unlike in the fixed $n_{\text{e}}$ case. For large $n_{\text{max}}$ the two
curves are identical.

Figure 5 panel (b) also informs us of how the system evolves under the
instability. Holding $\mu_{\text{e}}$ fixed, for a system at $T_{\text{crit}}$
and if no heating occurs, then the evolution of $n_{\text{max}}$ will occur
along a constant $T_{\text{crit}}/\mu_{\text{e}}$ line from higher to lower
$n_{\text{max}}$ i.e., the field will increase until reaching the branch of
$T_{\text{crit}}/\mu_{\text{e}}$ left of the peak at $n_{\text{max}}=48$. An
approximation for the amount the field could grow by the instability is found
by taking the two values of $n_{\text{max}}$ at which the
$T_{\text{crit}}/\mu_{\text{e}}$ line intersects the curve and using
$B\approx(\mu^{2}_{\text{e}}-m^{2}_{\text{e}})/(2en_{\text{max}})$. So for
fixed $\mu_{\text{e}}=25$ MeV and $T_{\text{crit}}/\mu_{\text{e}}=0.001$ i.e.,
$T_{\text{crit}}=2.9\times 10^{8}$ K, as the
$T_{\text{crit}}/\mu_{\text{e}}=0.001$ line intersects the curve at
$n_{\text{max}}=24$ and $n_{\text{max}}=119$, as the field evolves
isothermally it increases from $4.4\times 10^{14}$ G to $2.2\times 10^{15}$ G,
a factor of five growth. This growth will not occur continuously, as the fluid
stabilizes itself when not near a filled Landau level region of parameter
space, so this is an estimate of the maximum growth that could occur for a
particular region. This is of course a simplification, since (1) the
instability criterion is more complicated than assumed in this section, and
also depends on the relative strengths of components of the magnetic field;
(2) the evolution will not happen isothermally, and in fact the field growth
by the instability could lead to heating as we discuss in the next section.

(a)  (b)

Figure 6: Critical temperature $T_{\text{crit}}$ contours in K as a function
of $B$ and $\mu_{\text{e}}$ (left) and $B$ and $n_{\text{e}}(T=0)$ (right),
computed using the large $n_{\text{max}}$ approximation for $T_{\text{crit}}$
(Eq. (69)).

(a)  (b)

Figure 7: Critical temperature $T_{\text{crit}}$ contours in K as a function
of $B$ and $\mu_{\text{e}}$ (a) and $B$ and $n_{\text{e}}(T=0)$ (b), computed
using $T_{\text{crit}}/\mu_{\text{e}}=2/\beta_{\text{crit}}^{2}$ and the full
expression for $\beta_{\text{crit}}$, Eq. (68). The grey region is that for
which the right-hand side of Eq. (67) is negative, which is always stable.

The assumption that the electrons will be Landau quantized breaks down when
the energy spacing between successive Landau levels is near the thermal
energy, at which point the Landau levels are smeared out. The temperature at
which this occurs can be estimated using (Harding & Lai, 2006)

$T_{B}=\left(\sqrt{m_{\text{e}}^{2}+2eB(n_{\text{max}}+1)}-\sqrt{m_{\text{e}}^{2}+2eBn_{\text{max}}}\right).$
(70)

The contours of constant $T_{B}$ are plotted in Figure 8 as a function of $B$
and $\mu_{\text{e}}$ or $n_{\text{e}}(T=0)$. In general, at given location in
parameter space, $T_{B}$ will be about a factor of five larger than
$T_{\text{crit}}$, and hence the instability will always be inactive before
significant thermal smearing of the Landau levels occurs. As temperature is
increased we need to include higher Landau levels than $n=n_{\text{max}}+1$ in
our computation of $T_{\text{crit}}$, but the temperatures at which this is
necessary are above $T_{B}$, so since $T_{B}>T_{\text{crit}}$ the system will
already be unstable at these temperatures.

(a)  (b)

Figure 8: Contours of constant $T_{B}$, the temperature at which the spacing
between adjacent Landau levels equals the thermal energy and an indicator of
thermal smearing of the Landau levels. Values are in units of K. Plotted as a
function of $B$ and $\mu_{\text{e}}$ (a) and $B$ and $n_{\text{e}}(T=0)$ (b).

We note that the results here are limited by the approximations used, since
Section 4.1 assumed low temperatures $T\leq T_{B}$. Overall the most important
number for determining whether the instability can be active in a certain
region of $B$–$\mu_{\text{e}}$ or $B$–$n_{\text{e}}(T=0)$ parameter space is
$n_{\text{max}}$, and it must be greater than around $10$ for the instability
assuming the required partial derivatives are computed at fixed
$n_{\text{e}}$. Thus strongly quantizing fields $n_{\text{max}}\approx 0$ will
not be unstable. The most interesting range of $n_{\text{max}}$ for the
instability is between 10 and $\sim 100$, as $T_{\text{crit}}$ drops below the
expected temperatures for a magnetar crust as $n_{\text{max}}$ increases too
far beyond this point. Our analysis also applies to the core and Paper I,
since the instability condition there also depends on the negativity of
$1-4\pi\chi_{n}$. The range in parameter space that was examined here was
chosen for its relevance to the crust.

##  5 Physical implications

###  5.1 Enhanced Ohmic dissipation

Since the electric current density is related by Ampère’s Law to the curl of
$H$, not the curl of $B$, the Ohmic dissipation will be modified by the same
terms which lead to an instability. Estimating the current density as

$J^{i}_{e}=\frac{c}{4\pi}\epsilon^{ijk}\nabla_{j}H_{k}=\frac{c}{4\pi}\left[\frac{H}{B}\epsilon^{ijk}\nabla_{j}B_{k}-4\pi\left(\chi_{n}-\frac{M}{B}\right)\epsilon^{ijk}\nabla_{j}B\hat{B}_{k}\right],$
(71)

the Ohmic heating rate per unit volume is given by

$\dot{q}_{\text{O}}=\frac{1}{\sigma}J^{2}_{e}\approx\frac{c^{2}}{16\pi^{2}\sigma}\left[\epsilon^{ijk}\nabla_{j}B_{k}-4\pi\chi_{n}\epsilon^{ijk}\nabla_{j}B\hat{B}_{k}\right|^{2}$
(72)

where we used the approximations $|\chi_{n}|\gg|M|/B$ and $1\gg 4\pi|M|/B$ in
the second line. To estimate the relative size of this heating rate compared
to the $H=B$ result, assuming a straight field we obtain

$\dot{q}_{\text{O}}=\frac{c^{2}}{16\pi^{2}\sigma}\left(1-4\pi\chi_{n}\right)^{2}|\epsilon_{ijk}\nabla_{j}B_{k}|^{2}.$
(73)

So the Ohmic dissipation is enhanced by a factor of $(1-4\pi\chi_{n})^{2}$
compared to the result ignoring Landau quantization in regions where
$p_{F}^{2}/(2eB)$ is close to $n_{\text{max}}$. Even at high temperatures
$|4\pi\chi_{n}|$ can still be of order one as Figure 5 panel (a) shows, so the
heating can be increased by a factor $\gtrsim 4$ or more. The increased Ohmic
heating is possible for both signs of $\chi_{n}$, unlike the instability which
requires it to be negative. However, since this enhancement is only present in
limited spatial regions, the overall increase in Ohmic dissipation may be
modest.

This estimate of increased Ohmic dissipation has assumed that the conductivity
is unchanged by the Landau quantization of electrons, which is false. For
strongly-quantizing fields $\sigma$ can change by a factor $\sim 10$ compared
to its classical value due to Shubnikov–de Haas oscillations (Potekhin, 1999).
However, for $n_{\text{max}}$ in the unstable regime $\gtrsim 10$, $\sigma$
does not differ from its classical value by more than a factor of $\sim 2$.
Thus the enhancement of the Ohmic dissipation due to $\chi_{n}$ may not be
completely suppressed by an increase in the conductivity in those regions of
$B$–$n_{\text{e}}$ parameter space where we expect an instability. Unlike the
Shubnikov–de Haas oscillations and the instability, the increase in the Ohmic
dissipation is likely to persist at higher temperatures to some extent.

###  5.2 Domain formation and domain-formation heating

As discussed in previous studies of magnetized neutron star crusts (Blandford
& Hernquist, 1982; Suh & Mathews, 2010; Wang et al., 2013, 2016), regions of
parameter space where $\chi_{n}>1/(4\pi)$ are thermodynamically unstable to
the formation of magnetic domains. This Shoenberg effect (Shoenberg, 1984) and
the resulting Condon domains (Condon, 1966; Egorov, 2010) are well studied in
laboratory settings. The resulting magnetic domain-forming regions are not
composed of alternating domains of opposite magnetization as in ferromagnets.
This fact has sometimes been obscured in the literature by the tendency to
split the grand potential $\Omega$ and hence the magnetization derived from it
into “non-oscillatory” and “oscillatory” parts, the latter of which is split
into domains of opposite magnetization, but when added to the former the
overall magnetization of neighbouring domains will differ but point in the
same direction.

In the unstable regions, the equilibrium configuration is found using the
Maxwell construction for $H(B)$ at fixed $\mu_{\text{e}}$: $\mu_{\text{e}}$
must be held constant so that electrons can be exchanged between neighbouring
domains without energy cost. Working at $T=0$ for simplicity, the equilibrium
value of $H$ in the unstable region is found by finding the value of $H$ at
which the thermodynamic potential

$\overline{\Omega}=\Omega-\frac{H^{i}B_{i}}{4\pi},$ (74)

intersects itself. If the equilibrium value of $H$ in a domain-forming region
is $H_{\text{eq}}$, and the two values of $B$ such that $H(B)=H_{\text{eq}}$
and $(\partial H/\partial B)_{n_{\text{e}}}>0$ bounding the domain-forming
region are $B_{1}$ and $B_{2}$ where $B_{1}<B_{2}$, then the grand potential
density in the domain-forming region of parameter space at a stable
equilibrium is

$\Omega_{\text{eq}}=\Omega(B_{1})+\frac{1}{4\pi}H_{\text{eq}}(B-B_{1}).$ (75)

$\Omega_{\text{eq}}$ joins continuously to $\Omega$ outside this range of $B$.
The slope $(\partial\Omega/\partial B)_{\mu_{\text{e}}}=H/(4\pi)$ is also
continuous, but $(\partial H/\partial B)_{\mu_{\text{e}}}$ is discontinuous.
The Maxwell construction applied to $H(B)$ is shown in Figure 9, and the
energy difference $\Delta\Omega=\Omega-\Omega_{\text{eq}}$ between the
unstable and equilibrium configuration in a domain-forming region is shown in
Figure 10. Note that the sections of the original $H(B)$ curve included within
the domain-forming regions where $1-4\pi\chi_{n}>0$ or equivalently $(\partial
H/\partial B)_{n_{\text{e}}}>0$ are metastable, and domains do not necessarily
form in this region of parameter space absent nucleation sites for domain
walls.

Figure 9: $H(B)$ for fixed $\mu_{\text{e}}$ and $T=0$, showing how the Maxwell
construction is used to determine the equilibrium value of the field $H$ and
the span of $B$ for which $H$ is uniform. The construction has been applied to
two different domain-forming regions of parameter space distinguished by
superscripts $a$ and $b$. Shaded areas $a_{1}$ and $a_{2}$, and areas $b_{1}$
and $b_{2}$, are equal.

(a)  (b)

Figure 10: The difference in the grand potential density between the unstable
(no domains) and stable equilibrium (with domains) state at two distinct
domain-forming ranges of $B$ for fixed $\mu_{\text{e}}$ and $T=0$. Note the
difference fixed values of $\mu_{\text{e}}$ for panels (a) and (b).

Blandford & Hernquist (1982) discuss how the electrodynamic boundary
conditions at the no-domain to domain-forming region interface can modify the
picture discussed above. In the absence of surface currents, the
electrodynamic boundary conditions are

$\hat{n}_{i}B^{i}_{1}=\hat{n}_{i}B^{i}_{2},\quad\epsilon_{ijk}\hat{n}^{j}H_{1}^{k}=\epsilon_{ijk}\hat{n}^{j}H_{2}^{k},$
(76)

for $\hat{n}^{i}$ normal to the no-domain to domain-forming region interface
and $1$ and $2$ denote these distinct regions. Under the assumption that the
angle $\theta$ between $B$ and the interface varies only by a small amount
$\delta\theta\ll\theta$ between the two sides of the interface and that the
magnitudes of $H^{i}$ and $B^{i}$ are similar on each side, one can show that

$\frac{H_{2}-H_{1}}{B_{2}-B_{1}}\approx-\tan^{2}\theta,$ (77)

and hence the Maxwell construction as used above is only correct if $\theta=0$
and $H_{1}=H_{2}$. The correct Maxwell construction in the general case is

$\int^{B_{2}}_{B_{1}}H(B)\textrm{d}B=\frac{1}{2}(H_{1}+H_{2})(B_{2}-B_{1}),$
(78)

where $H_{1}=H(B_{1})$ and $H_{2}=H(B_{2})$ are the values of $H$ at the two
ends of the linear “equilibrium” solution for $H(B)$. $H_{2}\leq H_{1}$ by Eq.
(77), so if $H_{1}\neq H_{2}$, then the two points on the $H(B)$ curve which
are connected by the “equilibrium” solution are metastable, and the
“equilibrium” solution is still thermodynamically unstable since $(\partial
H/\partial B)_{n_{\text{e}}}<0$ in the range $B_{1}<B<B_{2}$. The case where
$H_{1}=H_{2}$ is energetically favourable and is the true equilibrium unless
additional forces not considered here (e.g., force balance between the Lorentz
force and crustal elasticity) prevent it from being realized.

Unlike the previous studies on neutron star crusts, we are less interested in
the magnetostriction associated with the magnetic domains and its ability to
e.g., “crack” the crust, than we are with a potential mechanism to release
heat during the domain formation process. However, a newly-born magnetar will
not initially contain magnetic domains due to their suppression at high
temperatures. As the star cools below $T_{\text{crit}}$ (which depends on $B$
and $n_{\text{e}}$), domain formation is possible if $B$ and $n_{\text{e}}$ in
a particular region of the star are such that a Landau level is nearly
empty/almost filled. The timescale to form these domains is of order the Ohmic
decay timescale (Blandford & Hernquist, 1982) as the domain formation process
is dominated by eddy currents.

The instability of both the Hall modes and the Alfvén modes as discussed in
Section 3 requires $\chi_{n}>1/(4\pi)$ to be true, so it is only possible in a
region of parameter space which is already thermodynamically unstable. Since
the timescale of this instability is shorter than the Ohmic decay timescale,
which is also the domain formation timescale, this suggests that unstable
field growth through this MHD instability could occur before domain formation
locally stabilizes the system. Since the growth time of the unstable Alfvén
modes is much shorter than for the unstable Hall modes, the unstable field
growth would be dominated by the former. Because this field growth would occur
over length scales of order the domain size, which is of order the spacing
between consecutive regions where new Landau levels become occupied $\sim
10-100$ cm, it could lead to enhanced Ohmic dissipation, since these length
scales are relatively short in the context of neutron star crusts.
Additionally, metastable phases of non-uniform magnetization can exist in the
absence of domain wall nucleation site. If disturbed, domain formation and a
release of energy is possible in these regions.

We consider the possibility of the dynamic destabilization of the crust due to
thermal or magnetic field evolution, and the resulting heating by Ohmic
dissipation of magnetic fields that grow unstably due to the Hall MHD
instability described in Section 3. We can make a “back-of-the-envelope”
estimate of the maximum heating that the domain formation process could
provide. The maximum energy per unit volume available to be released during
the domain formation process is $\Delta\Omega$. Not all of this will be
available to be converted into field growth and then decay to heat the crust,
as forming domain walls requires energy, though this will be of order a
percent of the total energy available (Wang et al., 2016).

The heating occurs in the unstable regions which are spaced throughout the
crust, which we approximate as spherical shells. Their radial extent $\Delta
z$ varies from order tens to hundreds of cm, becoming more compressed where
the field is weaker at fixed density or where the density is increased at
fixed field. To determine how many of these regions there is, we assume a
uniform field strength $B_{0}$ and note that the crust extends from
$\mu_{\text{e}}=m_{\text{e}}$ to $\mu_{\text{e},\text{cc}}\approx 70$ MeV, so
starting from the surface and moving inward successive Landau levels from
$n=0$ to

$n_{\text{max},\text{cc}}=\left\lfloor\frac{\mu_{\text{e},\text{cc}}^{2}-m_{\text{e}}^{2}}{2eB}\right\rfloor=\left\lfloor\frac{\mu_{\text{e},\text{cc}}^{2}-m_{\text{e}}^{2}}{11.83B_{15}\text{
MeV}^{2}}\right\rfloor,$ (79)

will be occupied. For $B_{15}=B/(10^{15}\text{ G})=1$,
$n_{\text{max},\text{cc}}\approx 400$, so there are 400 regions which are
potentially unstable to domain formation. We know from Section 4 that nonzero
temperature will suppress the instability in many of these regions, so we use
a reduced value $N_{\text{in}}\approx 100$. We also make the very optimistic
assumption that these regions are all destabilized simultaneously and that all
have similar values of $\Delta\Omega$.

The relevant timescale for the heating through this mechanism is the Ohmic
dissipation timescale

$\tau_{\text{O}}=\frac{4\pi L^{2}\sigma}{c^{2}},$ (80)

for length scale $L$ and conductivity $\sigma$, not the much shorter
instability timescale. We take $L\approx\Delta z$. The conductivity of a
neutron star crust ranges from $\sigma\approx 10^{17}$–$10^{27}$ s-1 (Potekhin
et al., 1999); for temperatures $\sim 5\times 10^{8}$ K, in the bulk of the
crust $\sigma\sim 10^{23}$ s-1. The heat flux from this mechanism is thus

$F_{\text{O}}\approx N_{\text{in}}\frac{\Delta\Omega\Delta
z}{\tau_{\text{O}}}\approx 3.3\times
10^{21}\left(\frac{N_{\text{in}}}{100}\right)\left(\frac{10^{24}\text{
s}}{\sigma}\right)\left(\frac{\Delta\Omega}{10^{23}\text{
erg/cm}^{3}}\right)\left(\frac{100\text{ cm}}{\Delta
z}\right)\frac{\text{erg}}{\text{cm}^{2}\text{ s}},$ (81)

assuming that the heating is continuous, though this is not true since it
ceases once the domains form. Beloborodov & Li (2016) estimated that for
magnetars to sustain their strong surface flux, the heating mechanism in their
crust must provide a heat flux of order $F\sim 10^{24}(\epsilon/0.01)^{-1}$
erg cm-2 s-1, where $\epsilon$ is an efficiency factor accounting for heat
loss to the crust and to neutrino emission. Thus, even based on our very
optimistic analysis, magnetic domain-associated Ohmic dissipation is far too
small to be of any significance as a crustal heating mechanism.

###  5.3 Location of domains

In the ultrarelativistic limit of massless electrons
$p_{F}\approx\mu_{\text{e}}$, which is an excellent approximation for most of
a neutron star crust, we can determine the location in $B$–$\mu_{\text{e}}$
parameter space of the magnetic domain-forming regions relatively easily.
Under this approximation we can write $P_{\text{e}}$ and $M$ as

$\displaystyle P_{\text{e}}={}$
$\displaystyle\frac{eBp_{F}^{2}}{2\pi^{2}}\sum_{n=0}^{n_{\text{max}}}\frac{g_{n}}{2}\left[\sqrt{1-bn}-bn\,\text{arsinh}\left(\sqrt{\frac{1}{bn}-1}\right)\right]=\frac{p_{F}^{4}}{4\pi^{2}}b\left[\Sigma_{1/2}(b)-\Sigma_{S,n}(b)\right],$
(82) $\displaystyle M={}$
$\displaystyle\frac{ep_{F}^{2}}{2\pi^{2}}\sum_{n=0}^{n_{\text{max}}}\frac{g_{n}}{2}\left[\sqrt{1-bn}-2bn\,\text{arsinh}\left(\sqrt{\frac{1}{bn}-1}\right)\right]=\frac{ep_{F}^{2}}{2\pi^{2}}\left[\Sigma_{1/2}(b)-2\Sigma_{S,n}(b)\right],$
(83)

where $b\equiv 2eB/p_{F}^{2}$, and where $\Sigma_{1/2}(b)$ and
$\Sigma_{S,n}(b)$ are given by $\Sigma_{1/2}$ and $\Sigma_{S,n}$ in Eq. (62a)
with $n_{\text{max}}\rightarrow b^{-1}$ and including the $n=n_{\text{max}}$
term in the sum. Note that these equations are functions of only a single
parameter $b$. For fixed $\mu_{\text{e}}$, Eq. (78) implies that

$\frac{P_{\text{e}}(B_{2},\mu_{\text{e}})-P_{\text{e}}(B_{1},\mu_{\text{e}})}{B_{2}-B_{1}}=\frac{1}{2}\left(M(B_{1},\mu_{\text{e}})+M(B_{2},\mu_{\text{e}})\right).$
(84)

Note that this is inconsistent with Eq. (17) of Blandford & Hernquist (1982):
we find that their Eq. (17) is only possible if the two possible
magnetizations in the domain-forming regions are of equal magnitude but
opposite direction, but this is not true in general.

Using Eq. (82–83), the locations of the domains as a function of $b$ can be
readily obtained. The upper and lower bounding values of $b$ of the domain-
forming regions are plotted as a function $n_{\text{max}}$ in Figure 11. The
constant values of $h=2eH/p_{F}^{2}$ in the domain-forming regions is also
shown. Since $b$ is a ratio between $B$ and
$p_{F}^{2}\approx\mu_{\text{e}}^{2}$, these values only need to be computed
once; then for a given $B$, they can be translated into the ranges of
$\mu_{\text{e}}$ where domain formation will occur or vice versa.

Figure 11: Scaled, shifted values of $b=2eB/p_{F}^{2}\approx
2eB/\mu_{\text{e}}^{2}$ as a function of $n_{\text{max}}$ at the boundaries of
the domain-forming regions, with subscript $+$ and $-$ denoting the upper and
lower boundary values respectively. The constant values of $h=2eH/p_{F}^{2}$
in the domain-forming regions, scaled and shifted by the same factors as the
$b$ values, are also displayed.

##  6 Discussion and Conclusion

Magnetohydrodynamic instabilities within magnetars are commonly invoked in
mechanisms for energetic magnetar outbursts, including many proposed fast
radio burst models. These instabilities could cause crust yielding and hence
magnetospheric twisting and bursting activity. The thermal and magnetic
evolution of magnetars are intricately linked, and the decay of strong surface
fields typical of magnetars possibly powers a heating mechanism responsible
for the higher luminosities observed for magnetars compared to rotation-
powered pulsars. We thus studied Hall MHD applied to a neutron star crust
threaded by magnetic fields in the range characteristic of magnetars, with the
goal of looking for instabilities that could affect magnetic field evolution
or crustal heating.

Like in the core, we find instabilities associated with the population of new
Landau levels originating from the differential magnetic susceptibility. The
novel instability in the crust is associated with the circularly-polarized
Hall modes, though like in the core we also find unstable Alfvén modes, which
grow much faster. Using linear perturbation theory we estimated the Hall mode
instability’s growth time, finding a slow-growing mode with timescale set
approximately by the inverse Hall mode frequency. In these same regions,
unstable Alfvén modes grow on timescales of order the local Alfvén mode
period. We showed that the instability is self-limiting; unstable field growth
will eventually stabilize the magnetized fluid locally and magnetic domain
formation will occur. It is also suppressed above a field-and density-
dependent critical temperature. This suggests the instability will only become
active after the magnetar cools sufficiently, to temperatures $\sim
10^{7}$–$10^{9}$ K depending on the location within the crust and the strength
of the magnetic field. In particular, we showed (see Figure 5 panel (b)) that
the peak value of the critical temperature for instability is about
$0.0015p_{F}$ where electrons are relativistic, and it happens at
$n_{\text{max}}=p_{F}^{2}/2eB=45$. This also means that at a given $B$ the
peak temperature for stability is where $p_{F}\approx\sqrt{90eB}$ and
therefore has a value of about $0.0015\sqrt{90eB}=8.4\times
10^{7}\sqrt{B/B_{q}}$ K where $B_{q}=m_{\text{e}}^{2}/e$.

We studied the possibility that the domain formation transition could be a
heat source in magnetar crusts. We estimated the Ohmic dissipation of magnetic
fields that grow by the Hall MHD instability described in this paper in
unstable regions of the crust. Even using an optimistic estimate, this
mechanism is far too intermittent, and the energy available to be converted
into magnetic field and then heat via Ohmic dissipation is many orders of
magnitude too small, to dominate magnetar crustal heating. We also showed that
the magnetization alone can increase the Ohmic dissipation by a factor of a
few in a strongly magnetized neutron star due to the oscillations of the
differential magnetic susceptibility. This effect is independent of the
instability and persists to some extent at temperatures at which the
instability is suppressed, and also to regions of field-density parameter
space where the instability is not active (i.e., the strongly-quantized
regime). However, this enhancement will only occur in limited spatial regions,
so its overall effect on magneto-thermal evolution of neutron stars may be
modest. Numerical simulations may be necessary to determine whether this
effect can have observable implications. An additional possible heating source
is viscosity in the plastic flow, though the uncertain strength of the
viscosity as discussed in Section 3.4 means that further study is necessary to
determine its overall significance.

##  7 Acknowledgements

P. B. R. was supported in part by the INT’s U.S. Department of Energy grant
No. DE-FG02-00ER41132.

## Data Availability

There are no new data associated with this article.

## References

  * Akgün et al. (2018) Akgün T., Cerdá-Durán P., Miralles J. A., Pons J. A., 2018, Mon. Not. R. Astron. Soc., 481, 5331
  * Ashcroft & Mermin (1976) Ashcroft N. W., Mermin N. D., 1976, Solid State Physics. Saunders College Publishing, Philadelphia
  * Baiko & Chugunov (2018) Baiko D. A., Chugunov A. I., 2018, Mon. Not. R. Astron. Soc., 480, 5511
  * Baym et al. (1971a) Baym G., Pethick C. J., Sutherland P., 1971a, Astrophys. J., 170, 299
  * Baym et al. (1971b) Baym G., Bethe H. A., Pethick C. J., 1971b, Nucl. Phys., A175, 225
  * Beloborodov (2017) Beloborodov A. M., 2017, Astrophys. J., 843, L26
  * Beloborodov & Levin (2014) Beloborodov A. M., Levin Y., 2014, Astrophys. J., 794, L24
  * Beloborodov & Li (2016) Beloborodov A. M., Li X., 2016, Astrophys. J., 833, 261
  * Biskamp (2000) Biskamp D., 2000, Magnetic Reconnection in Plasmas. Cambridge University Press, Cambridge, U.K.
  * Blandford & Hernquist (1982) Blandford R. D., Hernquist L., 1982, J. Phys. C Solid State Phys., 15, 6233
  * Chamel (2017) Chamel N., 2017, J. Low Temp. Phys., 189, 328
  * Chamel et al. (2021) Chamel N., Fantina A. F., Suleiman L., Zdunik J.-l., Haensel P., 2021, Universe, 7, 193
  * Chugunov & Horowitz (2010) Chugunov A. I., Horowitz C. J., 2010, Mon. Not. R. Astron. Soc. Lett., 407, 54
  * Condon (1966) Condon J. H., 1966, Phys. Rev., 145, 526
  * Cramer (2001) Cramer N. F., 2001, The Physics of Alfvén Waves. Wiley, Berlin
  * Cumming et al. (2004) Cumming A., Arras P., Zweibel E., 2004, Astrophys. J., 609, 999
  * Egorov (2010) Egorov V. S., 2010, Phys.-Usp., 53, 755
  * Friedman & Schutz (1978) Friedman J. L., Schutz B. F., 1978, Astrophys. J., 222, 281
  * Glampedakis & Andersson (2007) Glampedakis K., Andersson N., 2007, Mon. Not. R. Astron. Soc., 377, 630
  * Goldreich & Reisenegger (1992) Goldreich P., Reisenegger A., 1992, Astrophys. J., 395, 250
  * Gourgouliatos & Cumming (2015) Gourgouliatos K. N., Cumming A., 2015, Mon. Not. R. Astron. Soc., 446, 1121
  * Gourgouliatos et al. (2013) Gourgouliatos K. N., Cumming A., Reisenegger A., Armaza C., Lyutikov M., Valdivia J. A., 2013, Mon. Not. R. Astron. Soc., 434, 2480
  * Harding & Lai (2006) Harding A. K., Lai D., 2006, Reports Prog. Phys., 69, 2631
  * Hernquist (1984) Hernquist L., 1984, Astrophys. J. Suppl. Ser., 56, 325
  * Heyl & Hernquist (2005) Heyl J. S., Hernquist L., 2005, Astrophys. J., 618, 463
  * Horowitz & Kadau (2009) Horowitz C. J., Kadau K., 2009, Phys. Rev. Lett., 102, 191102
  * Irgens (2008) Irgens F., 2008, Continuum Mechanics. Springer, Berlin
  * Jones (1988) Jones P. B., 1988, Mon. Not. R. Astron. Soc., 233, 875
  * Jones (2003) Jones P. B., 2003, Astrophys. J., 595, 342
  * Kaminker et al. (2006) Kaminker A. D., Yakovlev D. G., Potekhin A. Y., Shibazaki N., Shternin P. S., Gnedin O. Y., 2006, Mon. Not. R. Astron. Soc., 371, 477
  * Kaspi & Beloborodov (2017) Kaspi V. M., Beloborodov A. M., 2017, Annu. Rev. Astron. Astrophys., 55, 261
  * Kittel (2005) Kittel C., 2005, Introduction to Solid State Physics, 8 edn. John Wiley & Sons, Hoboken, NJ
  * Kumar et al. (2017) Kumar P., Lu W., Bhattacharya M., 2017, Mon. Not. R. Astron. Soc., 468, 2726
  * Kwang-Hua (2018) Kwang-Hua C. W., 2018, Astrophys. Space Sci., 363, 184
  * Landau & Lifshitz (1970) Landau L. D., Lifshitz E. M., 1970, Theory of Elasticity, 2 edn. Pergamon Press, Oxford
  * Lander (2016) Lander S. K., 2016, Astrophys. J., 824, L21
  * Lander & Gourgouliatos (2019) Lander S. K., Gourgouliatos K. N., 2019, Mon. Not. R. Astron. Soc., 486, 4130
  * Levin & Lyutikov (2012) Levin Y., Lyutikov M., 2012, Mon. Not. R. Astron. Soc., 427, 1574
  * Li et al. (2016) Li X., Levin Y., Beloborodov A. M., 2016, Astrophys. J., 833, 189
  * Lu & Kumar (2018) Lu W., Kumar P., 2018, Mon. Not. R. Astron. Soc., 477, 2470
  * Lu et al. (2020) Lu W., Kumar P., Zhang B., 2020, Mon. Not. R. Astron. Soc., 498, 1397
  * Lyubarsky (2014) Lyubarsky Y., 2014, Mon. Not. R. Astron. Soc. Lett., 442, L9
  * Lyutikov (2003) Lyutikov M., 2003, Mon. Not. R. Astron. Soc., 346, 540
  * Lyutikov (2006) Lyutikov M., 2006, Mon. Not. R. Astron. Soc., 367, 1594
  * Lyutikov (2013) Lyutikov M., 2013, Phys. Rev. E, 88, 053103
  * Lyutikov (2015) Lyutikov M., 2015, Mon. Not. R. Astron. Soc., 447, 1407
  * Margalit & Metzger (2018) Margalit B., Metzger B. D., 2018, Astrophys. J., 868, L4
  * Mereghetti et al. (2015) Mereghetti S., Pons J. A., Melatos A., 2015, Space Sci. Rev., 191, 315
  * Metzger et al. (2019) Metzger B. D., Margalit B., Sironi L., 2019, Mon. Not. R. Astron. Soc., 485, 4091
  * Pearson et al. (2018) Pearson J. M., Chamel N., Potekhin A. Y., Fantina A. F., Ducoin C., Dutta A. K., Goriely S., 2018, Mon. Not. R. Astron. Soc., 481, 2994
  * Pons & Geppert (2010) Pons J. A., Geppert U., 2010, Astron. Astrophys., 513, L12
  * Pons et al. (2009) Pons J. A., Miralles J. A., Geppert U., 2009, Astron. Astrophys., 496, 207
  * Potekhin (1999) Potekhin A. Y., 1999, Astron. Astrophys., 351, 787
  * Potekhin & Chabrier (2018) Potekhin A. Y., Chabrier G., 2018, Astron. Astrophys., 609, A74
  * Potekhin & Yakovlev (2001) Potekhin A. Y., Yakovlev D. G., 2001, Astron. Astrophys., 374, 213
  * Potekhin et al. (1999) Potekhin A. Y., Baiko D. A., Haensel P., Yakovlev D. G., 1999, Astron. Astrophys., 346, 345
  * Rau & Wasserman (2021) Rau P. B., Wasserman I., 2021, Mon. Not. R. Astron. Soc., 506, 4632
  * Rheinhardt & Geppert (2002) Rheinhardt M., Geppert U., 2002, Phys. Rev. Lett., 88, 4
  * Rheinhardt et al. (2004) Rheinhardt M., Konenkov D., Geppert U., 2004, Astron. Astrophys., 420, 631
  * Shoenberg (1984) Shoenberg D., 1984, Magnetic Oscillations of Metals. Cambridge University Press, Cambridge, U.K.
  * Strohmayer et al. (1991) Strohmayer T., van Horn H. M., Ogata S., Iyetomi H., Ichimaru S., 1991, Astrophys. J., 375, 679
  * Suh & Mathews (2010) Suh I. S., Mathews G. J., 2010, Astrophys. J., 717, 843
  * Thompson & Duncan (1995) Thompson C., Duncan R. C., 1995, Mon. Not. R. Astron. Soc., 275, 255
  * Turolla et al. (2015) Turolla R., Zane S., Watts A. L., 2015, Reports Prog. Phys., 78, 116901
  * Viganò et al. (2013) Viganò D., Rea N., Pons J. A., Perna R., Aguilera D. N., Miralles J. A., 2013, Mon. Not. R. Astron. Soc., 434, 123
  * Wang et al. (2013) Wang Z., Lü G., Zhu C., Huo W., 2013, Astrophys. J., 773, 2
  * Wang et al. (2016) Wang Z., Lü G., Zhu C., Wu B., 2016, Publ. Astron. Soc. Pacific, 128, 1
  * Wang et al. (2019) Wang W., Zhang B., Chen X., Xu R., 2019, Astrophys. J., 875, 84
  * Waxman (2017) Waxman E., 2017, Astrophys. J., 842, 34