Magnetic-field evolution in a plastically-failing neutron-star crust

TL;DR

This paper models the plastic failure and magnetic-field evolution in neutron-star crusts, revealing surface motions and magnetic field weakening, with implications for understanding magnetar activity.

Contribution

It introduces a coupled model of plastic flow and magnetic evolution in neutron-star crusts, highlighting the role of plasticity in magnetic stress release and surface motion.

Findings

Surface motions of 10-100 cm/year occur after crust yields.

Plastic flow viscosity must be around 10^{36}-10^{37} g cm^{-1}s^{-1}.

Localized magnetic fields are weaker with plastic flow than without.

Abstract

Under normal conditions in a neutron-star crust, ions are locked in place in the crustal lattice and only electrons are mobile, and magnetic-field evolution is thus directly related to the electron velocity. The evolution, however, builds magnetic stresses that can become sufficiently large for the crust to exceed its elastic limit, and to flow plastically. We consider the nature of this plastic flow and the back-reaction on the crustal magnetic field evolution. We formulate a plane-parallel model for the local failure, showing that surface motions are inevitable once the crust yields, in the absence of extra dissipative mechanisms. We perform numerical evolutions of the crustal magnetic field under the joint effect of Hall drift and Ohmic decay, tracking the build-up of magnetic stresses, and diagnosing crustal failure with the von Mises criterion. Beyond this point we solve for the coupled evolution of the plastic velocity and magnetic field. Our results suggest that to have a coexistence of a magnetar corona with small-scale magnetic features, the viscosity of the plastic flow must be roughly $10^{36}-10^{37}\ \textrm{g cm}^{-1}\textrm{s}^{-1}$. We find significant motion at the surface at a rate of $10-100$ centimetres per year, and that the localised magnetic field is weaker than in evolutions without plastic flow. We discuss astrophysical implications, and how our local simulations could be used to build a global model of field evolution in the neutron-star crust.

Tables

Table 1. Table 1: Summary of the simulations, all of which have a total run time of 5 kyr. The first column is the run name, and subsequent columns are: the maximum value of magnetic field occurring in the crust; the reference state described by equations ( 46 )-( 53 ) (where S ​ C , D ​ C 𝑆 𝐶 𝐷 𝐶 SC,DC refer to Single Closed loop and Double Closed loop, and S ​ O , D ​ O 𝑆 𝑂 𝐷 𝑂 SO,DO the corresponding states with open field lines at the surface); the maximum value of the plastic viscosity; the maximum value of conductivity used; whether plastic flow is permitted (Y) or not (N); whether we allow for plastic flow as soon as the failure criterion is satisfied anywhere – effectively meaning a shallow failure (S) – or whether we wait for it to be satisfied in the deeper part of the crust (D); and whether the reference state is rescaled (Y) or not (N).

Model	$B_{\max }$	RS	${\nu }_{\max }$	$\sigma$	Plastic	Failure	Rescaling
	(${10}^{15}$G)		($1.6{\text{g cm}}^{-1}{\text{s}}^{-1}$)	(s-1)
H0	1	DC	-	${10}^{24}$	N	-	-
N0	1	DC	${10}^{38}$	${10}^{24}$	Y	S	N
P0	1	DC	${10}^{38}$	${10}^{24}$	Y	S	Y
H1	1.3	DC	-	${10}^{24}$	N	-	-
P1	1.3	DC	${10}^{38}$	${10}^{24}$	Y	S	Y
N1	1.3	DC	${10}^{38}$	${10}^{24}$	Y	S	N
H2	2	DC	-	${10}^{24}$	N	-	-
P2	2	DC	${10}^{38}$	${10}^{24}$	Y	S	Y
N2	2	DC	${10}^{38}$	${10}^{24}$	Y	S	N
H3	2	DO	-	${10}^{24}$	N	-	-
P3	2	DO	${10}^{38}$	${10}^{24}$	Y	S	Y
N3	2	DO	${10}^{38}$	${10}^{24}$	Y	S	N
H4	4	SC	-	${10}^{24}$	N	-	-
P4	4	SC	${10}^{38}$	${10}^{24}$	Y	S	Y
N4	4	SC	${10}^{38}$	${10}^{24}$	Y	S	N
H5	4	SO	-	${10}^{24}$	N	-	-
P5	4	SO	${10}^{38}$	${10}^{24}$	Y	S	Y
N5	4	SO	${10}^{38}$	${10}^{24}$	Y	S	N
H6	4	DC	-	$2.5\times {10}^{23}$	N	-	-
P6	4	DC	${10}^{38}$	$2.5\times {10}^{23}$	Y	S	Y
N6	4	DC	${10}^{38}$	$2.5\times {10}^{23}$	Y	S	N
P7	2	DC	${10}^{37}$	${10}^{24}$	Y	S	Y
N7	2	DC	${10}^{37}$	${10}^{24}$	Y	S	N
P8	2	DC	${10}^{36}$	${10}^{24}$	Y	S	Y
N8	2	DC	${10}^{36}$	${10}^{24}$	Y	S	N
P9	2	DC	${10}^{38}$	${10}^{24}$	Y	D	Y
N9	2	DC	${10}^{38}$	${10}^{24}$	Y	D	N

Equations

$$0=-\nabla P_0-{\rho }_0\nabla {\Phi }_0+\nabla \cdot {\mathcal{M}}_0,$$

$${\mathcal{M}}_{ij}=\frac{1}{4\pi }\left(B_i{B}_j-\frac{1}{2}{B}^2{\delta }_{ij}\right)$$

$$\nabla \cdot \mathcal{M}=\frac{1}{4\pi }(\mathbf{B}\cdot \nabla )\mathbf{B}-\frac{1}{8\pi }\nabla (B^2),$$

$$0=-\nabla P-\rho \nabla \Phi +\nabla \cdot \mathcal{M}+\nabla \cdot \text{mn@boldsymbol}\tau ,$$

$${\color[rgb]{0,0,0}{\rho\frac{\partial{\boldsymbol{v}_{\mathrm{pl}}}}{\partial{t}}+\rho(\boldsymbol{v}_{\mathrm{pl}}\cdot\nabla)\boldsymbol{v}_{\mathrm{pl}}}}=-\nabla P-\rho\nabla\Phi+\nabla\cdot(\boldsymbol{\mathcal{M}}+\mn@boldsymbol{\tau}_{\textrm{el}})+\nu\nabla^{2}\boldsymbol{v}_{\mathrm{pl}},$$

$$P\approx P_0,\Phi \approx {\Phi }_0.$$

$$\text{mn@boldsymbol}\tau ={\mathcal{M}}_0-\mathcal{M}.$$

$$\nu {\nabla }^2{\mathbf{v}}_{\mathrm{pl}}=-\nabla \cdot (\mathcal{M}-{\mathcal{M}}_{\text{el}}),$$

$${\tau }_{ij}={\tilde{\tau }}_{ij}+\eta {\delta }_{ij},$$

$${\tau }_{\mathrm{el}}\le \sqrt{\frac{1}{2}{\tilde{\tau }}_{ij}{\tilde{\tau }}_{ij}}=\sqrt{\frac{1}{2}({\tilde{\mathcal{M}}}_{ij}^0-{\tilde{\mathcal{M}}}_{ij})({\tilde{\mathcal{M}}}_{ij}^0-{\tilde{\mathcal{M}}}_{ij})}$$

$${\tilde{\mathcal{M}}}_{ij}={\mathcal{M}}_{ij}-\frac{1}{3}{\mathcal{M}}_{kk}{\delta }_{ij}=\frac{1}{4\pi }\left(B_i{B}_j-\frac{1}{3}{B}^2{\delta }_{ij}\right),$$

$$\nu {\nabla }^2{\mathbf{v}}_{\mathrm{pl}}=-\nabla \cdot (\tilde{\mathcal{M}}-{\tilde{\mathcal{M}}}_{\mathrm{el}}).$$

$$4\pi \nu {\nabla }^2{\mathbf{v}}_{\mathrm{pl}}$$

$$=({\mathbf{B}}_0\cdot \nabla ){\mathbf{B}}_0-(\mathbf{B}\cdot \nabla )\mathbf{B}-\frac{1}{3}\nabla (B_0^2-B^2).$$

$${\tau }_{\mathrm{el}}\le \frac{1}{4\pi }\sqrt{\frac{1}{3}{B}_0^4+\frac{1}{3}{B}^4+\frac{1}{3}{B}_0^2{B}^2-(\mathbf{B}\cdot {\mathbf{B}}_0{)}^2}.$$

$$\frac{\partial \mathbf{B}}{\partial t}=-c\nabla \times \mathbf{E},$$

$${\mathbf{v}}_e=-\frac{\mathbf{j}}{e{n}_e}=-\frac{c}{4\pi e{n}_e}\nabla \times \mathbf{B},$$

$$\mathbf{E}=-\frac{{\mathbf{v}}_e\times \mathbf{B}}{c}+\frac{\mathbf{j}}{\sigma },$$

$$\frac{\partial \mathbf{B}}{\partial t}=-\nabla \times \left(\frac{c}{4\pi e{n}_e}\left(\nabla \times \mathbf{B}\right)\times \mathbf{B}+\frac{c^2}{4\pi \sigma }\nabla \times \mathbf{B}\right).$$

$$\mathbf{j}=-e{n}_e({\mathbf{v}}_e-{\mathbf{v}}_{\mathrm{pl}}).$$

$$\mathbf{E}=\frac{1}{c}\left(\frac{\mathbf{j}\times \mathbf{B}}{e{n}_e}-{\mathbf{v}}_{\mathrm{pl}}\times \mathbf{B}\right)+\frac{\mathbf{j}}{\sigma }.$$

$$\displaystyle\frac{\partial{\boldsymbol{B}}}{\partial t}=-\nabla\times\bigg{(}$$

$$\displaystyle+\frac{c^{2}}{4\pi\sigma}\nabla\times{\boldsymbol{B}}\bigg{)}\,.$$

$$n_e=\frac{Z}{A}(1-x_{fn})\frac{\rho }{m_b},$$

$$\mathcal{R}\equiv \frac{(r-R_{cc})}{(R_{nd}-R_{cc})},$$

$$\tilde{\rho }\equiv \frac{\rho }{{\rho }_{cc}}=410{\left(1+\frac{R_{cc}}{R_{nd}}\right)}^2(1-\mathcal{R}{)}^2+0.004,$$

$$n_e=10^{36}(1.5{\tilde{\rho }}^{2/3}+1.9{\tilde{\rho }}^2){\text{ cm}}^{-3}.$$

$$\Gamma =\frac{Z^2{e}^2}{a_I{k}_{B}T}$$

$${\tau }_{el}=\left(0.0195-\frac{1.27}{\Gamma -71}\right)\frac{Z^2{e}^2{n}_I}{a_I}.$$

$${\tau }_{el}=5.1\times 10^{29}(0.4\tilde{\rho }+0.5{\tilde{\rho }}^3){\text{ g cm}}^{-1}{\text{s}}^{-2}.$$

$$\nu =1.6\times 10^{38}(0.4\tilde{\rho }+0.5{\tilde{\rho }}^3){\text{ g cm}}^{-1}{\text{s}}^{-1}.$$

Full text

††pagerange: Magnetic-field evolution in a plastically-failing neutron-star crust–References††pubyear: 0000

1 Introduction

Many of the diverse phenomena observed from neutron stars (NSs) are believed to be related to magnetic-field evolution in their solid crusts. Most dramatically, the activity of magnetars – bursts and flares in the hard X-ray/soft gamma-ray band, and their irregular enhanced spindown – are most naturally explained in a paradigm where the solid crust has to fail. The idea is that evolution of the star’s magnetic field induces a secular build-up of stresses in the star’s elastic crust, which eventually reaches its yield point; at this point there is a dynamical release of excess energy that heats the crust and produces coronal loops outside the star, in turn triggering emission processes that lead to the observed magnetar phenomena (Thompson & Duncan, 1995; Beloborodov & Thompson, 2007; Beloborodov, 2009).

Turning this picture into a quantitative model is, however, extremely challenging. It relies on detailed information about the star’s internal magnetic-field strength and geometry, in both the core and the crust, and the processes by which the field evolves. It demands an understanding of how a neutron-star crust fails, both at a microscopic and macroscopic level, and the role played by the magnetic field in this failure. Finally, it requires knowledge of where the released elastic energy is transferred after crustal failure, and how this induces emission processes. No one could reasonably claim to have all of these ingredients under control yet, but there have been many recent promising developments.

Firstly, the emission part of the problem is arguably understood to a relatively satisfactory extent. Thermal blackbody emission is associated with the hot crust; persistent X-ray emission is likely related to coronal loops (Beloborodov, 2009); timing changes and enhanced spindown to outer coronal regions (Beloborodov & Thompson, 2007).

As for field evolution, if one restricts oneself to considering the magnetic field of a (never-yielding) crust in isolation from the core, the problem reduces to the simpler framework of electron MHD. This limiting case is so named because the ions are assumed fixed, and so the electric current (and consequently the magnetic field) depend only on one variable: the electron velocity through the crust. Magnetic-field evolution proceeds under the joint effect of the conservative Hall drift and the dissipative Ohmic decay. Electron MHD has been used extensively and successfully over the last decade or so to make quantitative studies of the magnetic-field evolution which takes place in a NS crust, and this field of study is now quite mature and well-explored. Building on the classic study of evolution mechanisms by Goldreich & Reisenegger (1992) and some early analytic work (Cumming et al., 2004), numerical studies – in particular by Pons and collaborators – have made steady progress towards full magneto-thermal evolutions in 2D (Pons & Geppert, 2007; Aguilera et al., 2008; Pons et al., 2009; Viganò et al., 2013). Recently, full 3D studies of the problem have appeared, demonstrating the radical differences between poloidal- and toroidal-dominated fields (Wood & Hollerbach, 2015; Gourgouliatos et al., 2016; Gourgouliatos & Hollerbach, 2018), along with 3D studies of the generalised neutron-star induction equation Viganò et al. (2018). Akgün et al. (2018) have made a first effort to connect the crustal evolution to the magnetosphere (and therefore to observable emission).

Core-field evolution is a much more vexed issue, mainly due to the uncertain role played by superfluid components (which are present even in magnetars younger than any observed to date; Ho et al. (2012)). There are debates about whether any core process is fast enough to be relevant to magnetars; providing arbitration of the dispute is beyond the scope of this paper, but for a flavour of the recent literature see, e.g., Jones (2006); Glampedakis et al. (2011b); Graber et al. (2015); Gügercinoğlu & Alpar (2016); Dommes & Gusakov (2017); Passamonti et al. (2017); Castillo et al. (2017); Ofengeim & Gusakov (2018). There is also little understanding about how core and crust are linked, and the physics at the boundary that separates them. Even if core-field evolution is negligible, however, the magnetic field residing there is likely still important in providing the inner boundary condition for crustal-field evolution. In this case, one needs to know whether the majority of the core consists of type-I or type-II superconducting protons. In the latter case the microscopic field is quantised into thin fluxtubes (Baym et al., 1969), which may be averaged over to produce a macroscopic magnetic force (Easson & Pethick, 1977; Mendell, 1991; Glampedakis et al., 2011a), and the equilibrium equations may then be solved to find candidate macroscopic magnetic field configurations (Akgün & Wasserman, 2008; Lander, 2013; Henriksson & Wasserman, 2013; Lander, 2014). In the former case, the local field does not have a definite geometry, and it is less clear what form the macroscopic field would take (Sedrakian, 2005).

Regardless of where the field evolution takes place, however, understanding the response of the crust to high stress is pivotal: partially because the build-up of stresses is inevitable, and exceeding the yield limit seems highly likely (Thompson & Duncan, 1995; Perna & Pons, 2011; Lander et al., 2015; Lander, 2016; Li et al., 2016; Bransgrove et al., 2018), but also because models to explain magnetar emission rely on crustal motions (Beloborodov, 2009). The wide range of magnetar activity with a characteristic $0.1$-second timescale (e.g. Göğüş et al. (2001)), comparable with that for an elastic shear wave in the crust, also buttresses the argument that these phenomena are driven by a mechanism involving the transfer of crustal energy to the magnetosphere. The aim of this paper is to make the first quantitative study of the crust’s coupled magneto-plastic dynamics, studying the plastic flow induced by an excess of magnetically-induced crustal stress, and the back-reaction of this flow on the magnetic field evolution.

2 Neutron-star field evolution with plastic flow

2.1 Magnetic-induced plastic flow

To understand how magnetic stresses build in a neutron-star crust, and the crust’s non-elastic response beyond its yield limit, we take an approach based on generalising the crust’s momentum equation. This equation, whose form is effectively redundant in electron MHD, plays a key role in understanding the crust’s behaviour at high stress. Our approach, outlined in this subsection, is a revised presentation of the approach in Lander (2016).

In its early life, a neutron star is entirely fluid, with a magnetic field developed during a variety of complex processes at birth. This field should complete its rearrangement into a stable hydromagnetic equilibrium state over a few Alfvén timescales – a matter of seconds or less. This field is likely to be a dominantly large-scale structure dominated by low multipoles like the dipole component. The crust begins to freeze very gradually over the following years, and finishes the process after a far longer time: roughly 1000 yr. Because this timescale is so much longer than the Alfvén one, the crust must freeze in a relaxed equilibrium state with the magnetic field; initially the crust harbours no stresses. The force balance is therefore given by

[TABLE]

where

[TABLE]

is the Maxwell stress tensor, $P,\rho$ and $\Phi$ are, respectively, the pressure, mass density and gravitational potential; and the subscript zero indicates the value of quantities at the point when the crust freezes, and therefore when it is relaxed. The divergence of $\boldsymbol{\mathcal{M}}$ gives the familiar Lorentz force:

[TABLE]

where the first of the right-hand-side terms is the magnetic tension, and the latter the magnetic pressure. Because we are studying a problem involving crustal stresses, however, we will often find it clearer to use $\boldsymbol{\mathcal{M}}$ instead of ${\boldsymbol{B}}$.

Over time the crustal magnetic field changes due to the action of Hall drift and Ohmic decay, and therefore evolves away from its relaxed state; the force balance now becomes

[TABLE]

where $\mn@boldsymbol{\tau}$ is the elastic stress tensor. This is still a stationary equilibrium state, but only because the crustal lattice inhibits motion, at the cost of building up stresses.

Let $\mn@boldsymbol{\tau}_{\mathrm{el}}$ be the maximum stress which the crust is able to sustain before failure. This is a tensorial quantity and need not be unique; various different combinations of the individual tensor components could lead to failure. Once $\mn@boldsymbol{\tau}_{\mathrm{el}}$ is exceeded, we will assume the crust deforms through an irreversible plastic flow; we discuss the physics of crustal failure and a quantitative failure criterion in the following subsection. We will assume that any crustal failure is a simple shearing process; in this case the equation of motion for the crust gains an extra viscous-flow term corresponding to the piece of the elastic stress beyond the maximum value which the crust can sustain, and closely resembles the standard Navier-Stokes equation (Lander, 2016):

[TABLE]

Within this paper, therefore, the ‘plastic flow’ is physically the same as a standard viscous flow. Note, however, that if the failure geometry is less simple than the case we consider, the dynamics may depend explicitly on the (tensorial) rate of strain, through an extra term in the above equation (Prager, 1961).

The crust is stably stratified, so that any compressible motions are inhibited by buoyancy forces. We therefore demand that the flow satisfies $\nabla\cdot\boldsymbol{v}_{\mathrm{pl}}=0$ (note that this is less restrictive than assuming constant density). We anticipate plastic flow in a neutron-star crust to be slow and steady; slow meaning that we may neglect the advective term $(\boldsymbol{v}_{\mathrm{pl}}\cdot\nabla)\boldsymbol{v}_{\mathrm{pl}}$ for being quadratic in the velocity, steady meaning that we may neglect the acceleration term, $\dot{\boldsymbol{v}}_{\rm pl}=0$. Under these assumptions the plastic flow reduces to a standard Stokes flow, with just the viscous term being important.

We do not expect major changes in the hydrostatic equilibrium state over time, and so assume that

[TABLE]

Strictly speaking, these neglected changes are formally of the same order of magnitude as the plastic flow itself – quadratic in differences in the magnetic field. In proceeding, therefore, we should be aware that some small re-adjustment of the hydrostatic equilibrium could relieve some stresses and thus ultimately reduce the amount of energy transferred into the plastic flow compared with that predicted by our model. Note that if we had allowed for rotation, the corresponding centrifugal force term would be expected to vary considerably over time, since spindown at magnetar-field strengths can be very rapid.

With the above assumptions, we first see from equations (1) and (4) that the elastic stress (below the yield point) is sourced by Maxwell stresses alone:

[TABLE]

Next we compare the equations of motion above and below the yield stress, (5) and (4); subtracting one from the other leads us to the (unsurprising) result that changes in the Lorentz force greater than those sustainable by the crust induce a viscous flow:

[TABLE]

where $\boldsymbol{\mathcal{M}}_{\textrm{el}}$ is the value of the Maxwell tensor immediately before the yield stress is attained – recalling that this maximum tensorial stress is not uniquely defined. So far we have made concrete the notion of how magnetic stresses can eventually result in plastic motion of the crust, but have not described the mechanism of crustal failure in any detail. We will discuss this next, and in doing so will slightly revise the above equation.

2.2 Failure criterion

Arguably the biggest challenge in modelling any seismic phenomenon in a neutron star is that there are major uncertainties about how exactly the crystalline lattice of the crust yields at large stress. The terms ‘breaking’ or ‘cracking’ from older literature are misleading; at the high pressures relevant to a NS crust the material cannot break in the canonical fashion for terrestrial materials, in which a void propagates through the medium (Jones, 2003; Levin & Lyutikov, 2012). Microscopic simulations of crustal failure have provided valuable insight (Horowitz & Kadau, 2009; Hoffman & Heyl, 2012): they show that the crystal yields in a collective failure, undergoing an irreversible plastic deformation, but only once extreme levels of stress have built up (at least an order of magnitude higher than predicted by earlier work). Despite this result, we note that these simulations – performed by a gradual shearing applied at the boundaries of the local domain – may not unveil every failure mode of the system. It is also not clear whether crustal failure would proceed in the same way on larger scales. A failure event might appear collective on microscopic scales, whilst still being highly local with respect to the whole star. The plastically-flowing failed region may well be a wide patch or band – but a narrow failure band between two large non-failing pieces of crust could mimic ‘cracking’, notwithstanding the plastic nature of the localised failure.

The notion of collective failure and the resultant plastic flow is difficult to capture in a large-scale model of the crust, since the boundaries of the region which fail may be irregular, set by impurities or the crust’s seismic history. It is also not obvious what the ‘correct’ boundary conditions would be for the failed region(s). In this paper we will sidestep this highly challenging issue by considering a local, plane-parallel, portion of the crust: macroscopic, but not large enough that we need to account for the curvature of the crust, nor for coexisting yielding and non-yielding regions within our numerical domain. We will start with a prescription that once any gridpoint reaches failure criterion, the whole numerical domain – i.e. the entire segment of the crust we model – is assumed to fail. The plastic-flow velocity (8) is then calculated everywhere (even though most of the numerical domain is below its yield stress), and the resultant additional term is evolved in the induction equation. This notion of collective failure is faithful to the physical picture of failure from previous work. However, it is clearly a rather extreme case, since in practice the crust is liable to fail first in the outermost, lowest-density region. As a second prescription, we will instead only begin to apply the additional, plastic-flow parts of the evolution once a gridpoint further into the crust reaches its local yield stress. Fortunately, we will see later that our numerical results are rather insensitive to which presciption we use.

We will assume that failure occurs in the form of shearing motions, for a number of reasons. They satisfy our incompressibility condition; they result in the simplification of potentially complex viscoplastic dynamics to standard viscous dynamics (as discussed in the text around equation (5)); and they also allow us to adopt results from the aforementioned molecular-dynamics simulations with some degree of self-consistency (since they used shearing to induce failure). Furthermore, this class of motions are of most observational interest too, since they result in twisting of the magnetosphere and associated emission changes.

Terrestrial elasticity theory provides a variety of different relations used to predict yielding. A classic, simple yield criterion consistent with the physical picture we have for crust failing is that of von Mises, which is based on the concept that a material fails due to excessive shearing distortion rather than from volumetric changes. Baiko & Chugunov (2018) argue that this criterion is not reliable for a perfect Coulomb crystal, but as long as our local domain is large enough that we may regard the material as polycrystalline, it should be safe to proceed. In order to employ the von Mises criterion, we must first separate the stress tensor into two pieces:

[TABLE]

where $\tilde{\tau}_{ij}$ and $\eta\delta_{ij}=\frac{1}{3}\tau_{kk}\delta_{ij}$ are, respectively, the deviatoric and volumetric pieces of the stress tensor. The former, trace-free, component is the relevant one in the von Mises criterion, which states that a medium yields when

[TABLE]

for some scalar yield stress $\tau_{el}$ (this has no quantitative relation to the tensorial quantity $\mn@boldsymbol{\tau}_{\mathrm{el}}$ used in an indicative fashion in the earlier discussion).

In our model, elastic stresses build up due to changes in the Maxwell tensor (2) alone, and so we need to find the deviatoric piece of this. We do so from its definition as the traceless piece of the full tensor:

[TABLE]

which differs from the full tensor only in the prefactor of the diagonal piece. Now, with the understanding that yielding of the crust is sourced by the deviatoric piece of the Maxwell stresses alone, we see that we must revise equation (8) to read:

[TABLE]

This is now consistent with the failure criterion we adopt, but the numerical implementation of the above equation is problematic. In principle one can monitor elastic failure locally, by checking (at each timestep, and for each gridpoint) whether inequality (10) is satisfied. On the first timestep for which the inequality is satisfied at some gridpoint, one may record the components of the Maxwell stress deviator tensor to get the local value of $\tilde{\boldsymbol{\mathcal{M}}}_{\textrm{el}}$. Over time more gridpoints will reach the yield stress, and the results recorded to build up an array for $\tilde{\boldsymbol{\mathcal{M}}}_{\textrm{el}}$. However, in practical terms we cannot allow individual gridpoints to fail, since all the surrounding gridpoints provide a boundary condition to any motion and it would be physically improbable and numerically impossible to deal with; instead we expect larger regions to fail collectively. Nor can we complete building up the $\tilde{\boldsymbol{\mathcal{M}}}_{\textrm{el}}$ array before we allow any failure to occur, because the outer region would have vastly exceeded the breaking stress before anything occurred near the inner boundary. Instead, for the sake of definiteness, we will use the relaxed-state $\tilde{\boldsymbol{\mathcal{M}}}_{0}$ in lieu of $\tilde{\boldsymbol{\mathcal{M}}}_{\mathrm{el}}$, keeping in mind that doing so will likely result in an overestimated plastic flow. Having done so, our final form of the equation sourcing the plastic flow is

[TABLE]

One might have guessed that the plastic flow would depend only on differences in magnetic tension; the above relation shows that the changing magnetic pressure is also important. We may now also write the von Mises criterion in a more explicit form:

[TABLE]

The prefactors of these terms differ from those in (Lander et al., 2015), which calculated a failure criterion using the full Maxwell stress tensor (thus implicitly allowing crustal yielding to be sourced in part by the expansion/contraction of the crust). In the context of that paper the calculation was reasonable, since it was not concerned with what happens to the crust after yielding. Here, however, we have had to be more explicit about the failure mechanism, so that our calculation of crustal failure is consistent with the resulting plastic flow.

2.3 Electron MHD with plastic flow

To understand the evolution of the crustal magnetic field in the presence of a plastic flow, we need to revise the usual derivation leading to electron MHD (Goldreich & Reisenegger, 1992). We begin with Faraday’s equation:

[TABLE]

where ${\boldsymbol{E}}$ is the electric field and $c$ the speed of light. Now, electron MHD assumes that the advection of the magnetic field lines is due to the motion of free electrons, which may be thought of as the only carriers of the electric currents in the crust (since the ions, locked into the crustal lattice, are static). Thus, the electron velocity can be related to the magnetic field via Ampère’s law:

[TABLE]

where $n_{e}$ is the electron number density and $e$ the elementary electron charge. We can then evaluate the electric field through Ohm’s law:

[TABLE]

where $\sigma$ is the electric conductivity of the crust. Finally, we substitute into Faraday’s induction equation to obtain the usual electron MHD equation:

[TABLE]

The electron MHD equation holds provided the crust is strong enough to balance the force terms. Because of this, the momentum equation that appears in standard MHD (see equation (4)), is trivially satisfied. Should part of crust fail and consequently initiate a plastic flow, we need to take into account its effect on the overall evolution. We will associate the plastic flow with the motion of the ions alone, and make the approximation that the electrons are not themselves advected by the flow; for this, we require $|\boldsymbol{v}_{\mathrm{pl}}|\ll|{\boldsymbol{v}}_{e}|$. The new effect is that the electric current is no longer directly linked to ${\boldsymbol{v}}_{e}$ alone, but takes its more usual form as being due to a relative flow of oppositely-charged fluids:

[TABLE]

Ohm’s law – with or without the slow plastic flow – follows from the momentum equation of the electron momentum equation (Goldreich & Reisenegger, 1992). What changes in the case of plastic flow, is that eliminating ${\boldsymbol{v}}_{e}$ in favour of ${\boldsymbol{j}}$ now picks up a new $\boldsymbol{v}_{\mathrm{pl}}$ term not present before. As a result Ohm’s law becomes:

[TABLE]

Consequently the electron MHD equation with plastic flow becomes:

[TABLE]

2.4 Material properties of the crust

To solve the equations above, we need a model for the structure of a neutron-star crust and information about a number of its transport properties. One key quantity is the electron number density $n_{e}$, which is given by

[TABLE]

where $Z$ is the number of protons per ion, $A$ the number of baryons per ion, $x_{fn}$ the fraction of ‘free’ neutrons outside ions and $m_{b}$ the baryon mass (more specifically, we use the value of the atomic mass unit here). $Z,A$ and $x_{fn}$ are all functions dependent on the density alone, so first we need $\rho$. We obtain this by solving the TOV equations for a hydrostatic equilibrium configuration (Misner et al., 1973), using the equation of state of Douchin & Haensel (2001) based on the SLy4 interaction. We choose a configuration with a mass of $1.4$ solar masses, for which the stellar radius is $11.7$ km, the neutron-drip point is at $11.3$ km (i.e. the radius $R_{nd}$ of the density contour $\rho=\rho_{nd}\equiv 4\times 10^{11}\textrm{g cm}^{-3}$) and the crust-core boundary at $10.8$ km (i.e. the radius $R_{cc}$ of the density contour $\rho=\rho_{cc}\equiv 1.4\times 10^{14}\textrm{g cm}^{-3}$).

In the bulk of the literature on crustal dynamics, it is naively assumed that the relevant region of the star is the entire region for which $\rho<\rho_{cc}$. However, the outer crust ($\rho<\rho_{nd}$) only crystallises for temperatures below around $10^{9}$ K (Haensel & Pichon, 1994; Krüger et al., 2015) and so all or part of this region will be liquid for young neutron stars, or for those with an internal heating mechanism like magnetars (Kaminker et al., 2009). Furthermore, even if parts of the outer crust are formally just cold enough to solidify, the crystalline lattice would presumably be so weak that it would not be able to harbour any significant stress. For these reasons, in this paper we consider magnetic-field evolution in the inner crust only, with the outermost layer assumed to respond passively to this evolution, like a magnetosphere. It is therefore convenient for us to define a rescaled radius variable

[TABLE]

which is zero at the crust-core boundary and unity at the neutron-drip point.

Our numerical solution for $\rho$ is itself accurate, but as we will see later, some transport properties we require for this work are far more poorly understood. Therefore, rather than adopting the exact output from our calculation of $\rho$, we will use the following simple approximation:

[TABLE]

which deviates little from the exact result in the region $\rho_{nd}<\rho\leq\rho_{cc}$.

Douchin & Haensel (2001) provides tables of $Z,A,x_{fn}$ and $\rho$ throughout the crust. We use these results to construct polynomial fits to the former three quantities in terms of $\rho$, and from these quantities and equation (22) we calculate the electron number density $n_{e}$. The calculated result may be approximated to good accuracy with the following simple expression, which we use in place of the full numerical result:

[TABLE]

Whilst the crust remains below its yield stress, the equations of electron MHD apply, for which we only need to know one transport coefficient: the electrical conductivity $\sigma$. In this work we will assume $\sigma\propto n_{e}^{2/3}$, and set the conductivity at the base of the crust to $\sigma(R_{cc})=10^{24}$s*-1*. Next, the Coulomb parameter is given by

[TABLE]

where $a_{I}=(4\pi n_{I}/3)^{-1/3}$ is the ion sphere radius. This quantity is of fundamental importance: it tells us when it is energetically favoured for the outer envelope of the neutron star to crystallise. In particular, a lattice forms when $\Gamma>175$.

For the breaking stress we use the fit of Chugunov & Horowitz (2010) to their molecular-dynamics simulations of NS crust-yielding:

[TABLE]

It is again convenient to use a fit to these results instead of the full output. We adopt the following approximation:

[TABLE]

Finally, we need an expression for the viscosity $\nu$ of the crust in its plastic phase. This is essentially completely unknown from first principles, and so we will begin by using the estimate from Lander (2016), which comes from demanding that a magnetar’s corona can in principle be a persistent rather than a transient phenomenon. In other words, the viscosity must be low enough that crustal motions can twist the magnetosphere on a shorter timescale than the $10$-yr timescale on which exterior electric current is thought to dissipate (Beloborodov & Thompson, 2007). This reasoning leads to a plastic viscosity of the same functional form as that for $\tau_{el}$, with only a different magnitude:

[TABLE]

Note that a calculation of this quantity from theoretical grounds has been attempted by Kwang-Hua (2018), and is orders of magnitude lower than the above. If this were correct, plastic evolution would be so fast that the crust would be immediately returned to its elastic regime, and it seems doubtful that any large-scale gradual magnetospheric twisting could occur.

2.5 Boundary conditions

We have argued above that the natural outer boundary for our work is the boundary between the inner and outer parts of the crust, defined by the neutron-drip isopycnic contour where $\rho\approx 4\times 10^{11}$ g cm*-3*. We will assume that this boundary is free for both the magnetic field and the plastic velocity, and implement this by demanding that the radial derivatives of both quantities vanish there (rather than the quantities themselves). The magnetic field beyond this point is matched to a vacuum field, and any dynamics of the outer crust are assumed to be slave to the evolution occurring at greater depths.

The inner boundary is the crust-core boundary, at which we assume the radial magnetic magnetic field and plastic velocity to vanish. For the velocity this is quite natural; the stronger lattice at high densities is less susceptible to fail anyway, and if it does the motion will be much slower, due to the higher plastic viscosity there. The magnetic-field condition is standard and given little attention, but is harder to justify; the pragmatic reason to adopt it is to ensure that is no coupling between the crust and the (poorly-understood) magnetic-field evolution at the crust-core boundary and into the core. Since timescales for core evolution tend to be considerably longer than those of the crust, however, there is some reason to believe that the crustal evolution would be somewhat decoupled from core processes. Setting the radial field component to zero is one way to ensure this, albeit by brute force. For the tangential field component, one could demand that the tangential electric field vanishes, in accordance with the ‘Meissner boundary condition’ described in Hollerbach & Rüdiger (2002). However, this leads to the inversion of a non-linear equation, drastically complicating the calculation and only leading to minimal difference from the straightforward zero-boundary condition that we implement here.

2.6 Energy conservation and surface motions

Three different effects drive the magnetic-field evolution we consider here, and may lead to changes in the star’s total magnetic energy

[TABLE]

The properties of two of these effects are well known: Ohmic decay results in the dissipation of magnetic energy, whilst the Hall drift effect is conservative (Goldreich & Reisenegger, 1992). The behaviour of the final term, due to plastic flow, is less clear. We know, from equation (13), that the plastic flow is sourced from the unbalanced Lorentz force in a process involving the plastic-phase viscosity $\nu$ – but despite its connotations of dissipation, here $\nu$ acts more to regulate the transfer of energy from magnetic to (plastic) kinetic, rather than to dissipate it.

To gain some intuition for the effect of plastic flow on the total crustal energy, let us first consider a different but related effect: that of the dynamical term in the usual induction equation for MHD (i.e. in a fluid, not an elastic solid). With this term alone, the induction equation reads

[TABLE]

and under the action of this term alone the magnetic energy of the system evolves as follows:

[TABLE]

where we have used the product rule. This term does not appear to conserve magnetic energy, but this is not surprising: it describes a dynamical coupling between ${\boldsymbol{v}}$ and ${\boldsymbol{B}}$ which leads to (e.g.) Alfvén waves. Let us, then, instead check whether the term preserves total energy. In the simplest possible set-up, assuming the velocity evolves in a way independent of the stellar pressure and gravity, we have

[TABLE]

and changes in the total energy should be given by the sum of changes in the kinetic and magnetic pieces:

[TABLE]

where the first two terms in the final equality come from the magnetic-energy piece, upon using a couple of vector identities, and the final piece comes from substituting in equation (33) for $\dot{\boldsymbol{v}}$. Clearly the second and third terms cancel, leaving us (after using the divergence theorem) with a surface integral:

[TABLE]

where $\hat{\boldsymbol{n}}$ is the unit outward-pointing normal vector to the surface. So, if there are no losses from the surface due to electromagnetic waves, then the dynamical term of equation (31) is conservative.

Returning to the plastic-flow problem, as we have formulated it, we see that there is no such cancellation of terms. We are assuming that the flow is steady on the timescales of interest, $\dot{\boldsymbol{v}}_{\rm pl}=0$, meaning that the kinetic energy is roughly constant and only the magnetic energy is time-varying:

[TABLE]

Comparing with equation (34), we see that this may be written as

[TABLE]

where the second term gives the power of the plastic-flow process: the work done by the Lorentz force at a rate set by the plastic-flow velocity. From the above, we see that plastic flow is neither conservative or dissipative, a priori. For it to be conservative, we require (using another vector identity) that:

[TABLE]

We see that the power of the process – involving the interior plastic flow calculated from equation (13) – must source a surface plastic flow that satisfies the above relation.

Clearly it is not possible to satisfy the above relation unless $(\boldsymbol{v}_{\mathrm{pl}}\times{\boldsymbol{B}})$ has some component parallel to $({\boldsymbol{B}}\times\hat{\boldsymbol{n}})$ at the surface. By the definition of $\hat{\boldsymbol{n}}$ we know that $({\boldsymbol{B}}\times\hat{\boldsymbol{n}})$ is tangent to the surface. For $(\boldsymbol{v}_{\mathrm{pl}}\times{\boldsymbol{B}})$ to have a tangential component, either $\boldsymbol{v}_{\mathrm{pl}}$ or ${\boldsymbol{B}}$ must have a component parallel to $\hat{\boldsymbol{n}}$. In other words, necessary – but not sufficient – conditions for conservative magnetic-field evolution under plastic flow are:

[TABLE]

and:

[TABLE]

The most important implication of this is that once a neutron star crust yields under magnetic stresses, it must generate surface motions in order to conserve energy. If there is no surface motion, all of the power must be expended in some other way – presumably through viscoplastic heating in the crust.

3 Numerical implementation

We numerically integrate equation Eq. (21) (and, above the yield stress, equation (13)) in Cartesian geometry with orthogonal coordinates $(x,y,z)$. We consider a square block in $y-z$ of size $y\in[0,R_{nd}-R_{cc}]$ and $z\in[0,R_{nd}-R_{cc}]$, where $z$ is the depth of the crust and $z=0$ corresponds to the crust-core boundary. The system has plane-parallel geometry with the $x-$derivatives of the physical quantities vanishing. For this purpose we are using a modified version of a code previously developed and used for electron MHD evolution in Gourgouliatos et al. (2015); Gourgouliatos & Hollerbach (2016).

Some classes of magnetic field do not evolve under the Hall effect at all; these are known as Hall equilibria (Gourgouliatos et al., 2013). In order to ensure our initial conditions cause immediate evolution, we take care to ensure that they are deliberately out of Hall equilibrium. This initial field corresponds to the crust-freezing state, so that although it has non-vanishing Maxwell stresses $\boldsymbol{\mathcal{M}}_{0}$, these are in hydromagnetic equilibrium, and so the crust is unstressed with them. As the magnetic field evolves, we monitor the von Mises criterion (Eq. (14)) to see whether it is satisfied at some location within the integration domain. Since we are considering a block of the crust, we assume that if the crust fails at some point this will initiate a global plastic flow within the block. To find the plastic flow velocity we then solve Eq. (13).

3.1 Plane-parallel form of equations

We express the magnetic field in terms of two scalars

[TABLE]

This expression ensures that the magnetic field is divergence-free and incorporates its plane-parallel symmetry. The plastic flow is incompressible; $\nabla\cdot{\boldsymbol{v}}_{pl}=0$. Moreover, radial motions are energetically disfavoured, and thus we take the $z-$component of the velocity to be identically zero. Finally, the plane-parallel assumption we have made leads to $\partial_{x}{\boldsymbol{v}}_{pl}={\bf 0}$. Subject to these constraints the plastic velocity can be written in the following form

[TABLE]

Next we rewrite the electron MHD equation with plastic evolution (Eq. 21), which reduces to two scalar partial differential equations:

[TABLE]

[TABLE]

We apply periodic boundary conditions at $y=0$ and $y=R_{nd}-R_{cc}$, which are implemented by a series of ghost-points on either side of the boundary. At the crust-core boundary we apply the zero boundary condition, therefore $\Psi(y,0)=B_{x}(y,0)=0$. Finally, beyond $z=R_{nd}$ we consider the magnetic field to be current-free. Thus, $B_{y}(y,z>R_{nd})=0$, and further demanding that the Laplacian of $\Psi$ vanishes: $\Psi_{yy}(y,z>R_{nd})+\Psi_{zz}(y,z>R_{nd})=0$, which is implemented through two rows of ghost points.

A plastic flow will appear only if the crust fails. We scan through the crust checking whether the von Mises criterion is satisfied; if this is the case we integrate Eq. (13). This velocity is then included in the numerical integration of Eq. (44). We note that the $y-$component has only a $z-$dependence because of its zero-divergence, whereas the right-hand-side of Eq. (13) in principle depends on both $(y,z)$. To reconcile this issue we take the $y-$average of Eq. (13), which leads to an equation where both sides depend on $z$ only.

Numerically solving for $\boldsymbol{v}_{\mathrm{pl}}$ at every timestep would drastically increase the computational time required to evolve the system. For this purpose, we repeat the survey of crust failure and solution for $\boldsymbol{v}_{\mathrm{pl}}$ not at every integration time-step, but less frequently. For the typical resolution employed here ($100^{2}$), we found that updating $\boldsymbol{v}_{\mathrm{pl}}$ every $10^{3}$ steps gives the same results as solving for $\boldsymbol{v}_{\mathrm{pl}}$ at every time-step, and hence we have adopted this for our calculation.

3.2 Numerical tests

We have run the code experimenting with various initial conditions and resolutions. For all initial conditions, we have run simulations where the field evolves only via electron MHD. These simulations exhibit the basic characteristics of such studies: the initial evolution is drastic, essentially responding to the non-equilibrium starting state. Following this initial phase, the system later relaxes to a Hall-attractor. The pure electron MHD state serves as our benchmark simulation for a crust that does not yield and will be used for comparison with the runs with plastic flow. We have repeated some runs at higher resolution ($200^{2}$) and we found that the overall evolution is in good agreement with the integration times we present in this work. In the calculation we have included a finite resistivity $\sigma$ which leads to long-term decay of the magnetic field and an overall decrease in the magnetic energy of the field. The relative importance of Hall and Ohmic terms in the evolution may be determined from the magnetic Reynolds number (sometimes known as the Hall parameter):

[TABLE]

While $R_{M}$ is not constant throughout the crust, in the range of parameters employed we have a maximum value in our simulations in the order of few times $10^{2}$. We remark that in the setups with the strongest magnetic fields the maximum value was close to $10^{3}$ and the simulations were becoming numerically unstable. There we adopted a resistivity four times our standard value, to ensure numerical stability and to give a manageable $R_{M}\sim 10^{2}$.

3.3 Magnetic-field reference state

The magnetic field must reach an equilibrium state whilst the star is still fluid, so that when the crust first freezes it does so with the field in a particular state. It is born relaxed, naturally. Later it will deviate from this and so stresses build up. In studies of electron MHD without crustal failure the late-time evolution of the crustal field was somewhat independent of the initial state, with these early features having been washed out over time. One novel feature about the problem we study here is that the initial-state magnetic field remains important long after the crustal field has evolved, because it continues to set the crust’s relaxed state. In other words, in principle one could remove all elastic stresses and the crust would return to this state (or vice-versa). We do not anticipate that plastic flow will allow significant adjustment of the reference state, since it does not act to relieve all crustal stresses – only those exceeding the yield limit.

The initial conditions we implement are motivated by the presence of a strong crustal field, which appears in the form of loops resembling the structure of a poloidal loop field. We primarily consider an initial configuration with a double loop of closed field lines:

[TABLE]

[TABLE]

Here $B_{0}$ is a parameter that sets the magnetic field strength, which is chosen so that the maximum value of the magnetic field is $\sim 10^{15}$G.

To verify that the conclusions we draw from our evolutions are general, and not peculiar to this choice of initial magnetic field, we will also try other initial conditions. One of these is a single-loop configuration:

[TABLE]

[TABLE]

although this state, like the previous double-loop configuration, features a magnetic field that vanishes on the surface. We will, therefore, also relax this assumption by allowing for some of the initial field lines to extend beyond $R_{nd}$, both in the double loop

[TABLE]

[TABLE]

and the single loop configuration:

[TABLE]

[TABLE]

In the non-dissipative problem, these initial conditions would set the reference state. However, in some models the dissipation is higher than what it would be for a realistic neutron star. This is primarily because we simulate a slab instead of the whole crust, thus large-scale configurations which should be long-lived are confined to the size of the simulation box – and this sets the length-scale of the system and leads to faster dissipation. Secondly, in some case a higher resistivity is used than the realistic approximation. This could cause considerable energy loss, i.e. a dissipation of half of the magnetic energy within 1 kyr. Thus the initial conditions do not represent a state which can ever be attained again by the field, since there is simply not enough magnetic flux left. Numerically, we have found that a permanent stress – whose location and magnitude do not evolve significantly after a while – develops independently of the details of the evolution. For this reason it is necessary to update the reference state, so that it represents a configuration with the same energy as the present-day one, with just the geometry being different. We do this by rescaling the $\boldsymbol{\mathcal{M}}_{0}$, multiplying it by $E(t)/E(0)$ where $E(t)$ is the total magnetic energy of the simulated domain at time $t$ and $E(0)$ is its initial value. This is better than pumping up the evolving magnetic field because of possible numerical instabilities. The only thing to bear in mind is that whilst the whole magneto-plastic evolution is made self-consistent in this way, the characteristic timescale grows due to the reduction of $B$. So, the evolutions we see slow down at later times, despite the physical time moving forward at the same rate.

4 Results

As the parameter space of our magneto-plastic evolution models is rather large we will explore some characteristic cases to illustrate possible avenues for the magnetic field evolution through electron MHD and plastic flow. In our study we focus on the following main parameters: magnetic field strength and structure, magnitude of the plastic viscosity $\nu$, failure criterion and evolution of the reference state. The simulations we have performed are detailed in table 1.

For all models with a given magnetic field strength and structure we run a ‘Hall’ simulation not permitting any plastic flow (H0-H6), a ‘plastic’ simulation with plastic flow and where the reference state is rescaled (P0-P9) with the overall magnetic field decay, and a ‘non-rescaled’ plastic simulation (N0-N9). We have varied the maximum strength of the magnetic field and its structure; see the simulations with numbers 0 through 6. We take $\nu_{max}=1.6\times 10^{38}$g cm*-1* s*-1* as our standard value for the viscosity of the crust’s plastic phase, but given the large uncertainties in the physical properties of the crust we have also experimented by decreasing its value by one or two orders of magnitude, to $\nu_{max}\sim 10^{37}$g cm*-1* s*-1* (runs P7 and N7) and $\nu_{max}\sim 10^{36}$g cm*-1* s*-1* (runs P8 and N8). (From now on, we suppress mention of the clumsy $1.6$ prefactor on the maximum plastic viscosity.) We have also explored the possibility of collective failure occurring only once the von Mises criterion is satisfied deep within the crust, i.e. in the inner half of the simulated crust slab $z<(R_{cc}+R_{nd})/2$ (runs P9, N9), as opposed to our standard assumption that failure occurs when any gridpoint satisfies the failure criterion – effectively meaning a shallow failure.

We simulate all models for a timescale of $5$ kyrs to allow for sufficient time for the system to evolve fully. We restrict the timescale to this level, as beyond this timescale the evolution practically saturates to a Hall-attractor (Gourgouliatos & Cumming, 2014) state.

4.1 Crust failure

We find that for the weakest magnetic field strength simulated (H0, P0, N0), the von Mises criterion is never satisfied, thus no failure occurred during the simulation. An interesting point to notice for this system is that even if the magnetic field were to be completely removed by the evolution, the von Mises criterion would still barely be satisfied.

Setups with stronger initial magnetic field magnetic field (H1, P1, N1) demonstrate the importance of applying the rescaling to the magnetic field reference state, to compensate for the decay of the magnetic field. In run P1 the von Mises criterion is not satisfied and the system does not initiate a plastic flow. On the contrary, in run N1, where the reference state is not rescaled, a plastic flow is initiated after $300$ years. This is mainly due to the fact that approximately $15\%$ of the magnetic field energy has dissipated and consequently the von Mises criterion is more easily satisfied. Similarly, in simulation P2, stresses hardly exceeds $\tau_{el}$ anywhere in the domain, whereas in simulation N2 where no rescaling is applied the stresses are locally $2.5$ times higher than $\tau_{el}$. To show the effect of the rescaling on the evolution, we compare models with and without this effect in Fig. 2, finding that plastic flow is more effective in reducing stresses beyond the von-Mises limit when rescaling is included – a physically more reasonable result than the stresses persisting indefinitely.

The family of simulations consisting of a double loop confined inside the star where the maximum magnetic field is $2\times 10^{15}$ G (P2-N2, P7-N7, P8-N8) initiate a plastic flow within the first 40 years. By constrast, simulations with a double loop where the initial state has open field lines develop a plastic flow much earlier at around 10 years (P3-N3). Simulations with a single loop reach the plastic regime more slowly: at around 80 years if the initial flux is contained within the crust (P4-N4) and around 50 years if there are also open field lines (P5-N5). We can understand this by noting that while the single loop initial conditions have the same maximum value of the magnetic field the electric currents are weaker, due to the simpler structure of the initial field and thus the overall stresses. Similarly, the models where the initial magnetic field has field lines emerging from the surface of the crust, yield earlier due to the stronger initial field near the surface. The simulations P9 and N9, where we demand the von Mises criterion to be satisfied in the inner half of the crust before allowing for failure (i.e. plastic flow), take 80 years to fail – twice as long as the same family of simulations where the von Mises criterion need only be satisfied at shallow depths (P2-N2).

Finally, we note that the simulations where the reference state is not rescaled will continuously fulfil the von Mises criterion once they fail, since they cannot get back to the initial conditions because of the magnetic energy that has since been lost; this is the case for all N-simulations except for the N0 run. On the contrary, some of the simulations where the reference state is rescaled (P2, P3, P7, P8 and P9) tend after some time, ranging from 500-2500 years, to stop fulfilling the von Mises criterion (Eq. 14). This is primarily due to the fact that the maximum stress they developed was only a factor of about 1.5 that of the (evolving) $\tau_{el}$, and it became possible for the plastic flow to return the magnetic-field structure to a state where its stresses could be fully absorbed by the crust again. On the contrary, in simulations P4, P5 and P6, the von Mises criterion is satisfied with the right-hand-side being a factor of 2-5 higher than $\tau_{el}$. In this case, the magneto-plastic evolution cannot push the system sufficiently close to the reference state to choke the plastic flow.

4.2 Plastic flow

In all simulations where a plastic flow develops we find that it tends to slow the pure electron MHD evolution. In particular, the plastic velocity that develops tends to oppose the electron flow velocity, whenever this is possible – as $\boldsymbol{v}_{\mathrm{pl}}$ is subject to more constraints than ${\boldsymbol{v}}_{e}$. In all simulations we have deliberately set $B_{x}(t=0)=0$, thus this component can be thought of as a proxy for the efficiency of electron MHD evolution: the initial electron velocity distribution winds up the poloidal $B_{y}$–$B_{z}$ components of the field and develops a toroidal $B_{x}$ field. If the plastic flow is in the opposite direction from the electron flow, the $B_{x}$ component that will develop will be weaker than pure electron MHD.

In Fig. 3 we compare the state of the magnetic field at $t=500$ yr in simulations H2, P2, P7, P8 and P9. Since all these have identical initial conditions we can deduce the impact of the plastic flow and the role of the value of plastic viscosity. The system that evolves entirely via electron MHD develops a maximum for $B_{x}$ of about $1.2\times 10^{15}$ G (Fig. 3 b). If a plastic flow is permitted, then the maximum value of the field drops, with the drop being dramatic for decreasing plastic viscosity: for $\nu\sim 10^{38}\textrm{g\ cm}^{-1}\textrm{s}^{-1}$ the maximum value is $10^{15}$ G (Fig. 3 c), whilst decreasing this parameter by one and two orders of magnitude results in the maximum field strength dropping to $6.5\times 10^{14}$ G (Fig. 3 d) and $2\times 10^{14}$ G (Fig. 3 e) respectively. Allowing for plastic flow only when the inner part of the crust fails (Fig. 3 f) gives only a marginal difference from the case where the failure criterion is applied globally. In all cases we notice that plastic flow leads to smoother magnetic field lines and consequently weaker electric currents.

As expected, the plastic flow relieves part of the stress on the crust, as can be seen by comparing the normalised stresses $\tau/\tau_{\rm el}$ of models H2, P2 , P7 and P8 – see Fig. 4. Note that it is not always the system with the lowest $\nu$ that has least stress – in these particular snapshots, it is the system with $\nu\sim 10^{37}\textrm{g\ cm}^{-1}\textrm{s}^{-1}$.

The velocities of the plastic flow depend mainly on the plastic viscosity chosen and the initial magnetic field strength, Fig. 5 and 6. Comparing simulations P2, P6, P7 and P8 we notice that a drastically lower plastic viscosity leads to a higher velocity, e.g. compare the maximum values of plastic velocities of $10$cm/year for P2 to $130$cm/year for P8. Similarly, stronger magnetic fields leads to faster plastic flows: doubling the magnetic field strength from simulation P2 to P6 we notice an increase by a factor of 3 on the plastic flow velocity. These maximum velocities occur at the upper layers of the crust, where the plastic viscosity is lower. Broadly speaking, from the snapshots of Fig. 5 we see a strong azimuthal plastic flow $v_{\rm pl}^{x}$ close to the surface, and a weaker internal flow in the opposite direction deep within the crust, though the detailed properties of the flow evolve in time. The $v_{\rm pl}^{y}$ flow is simpler, an artefact of our need to perform an average for consistency with the plane-parallel geometry we have imposed, and is constant for a given depth within the crust. On a sphere it would represent a flow in the $\theta$-direction; within our local domain it simply flows from left to right, or vice-versa, with respect to the viewpoint of the plots of Fig. 6.

4.3 Energy dissipation

We established in equation (40) that a magneto-plastic evolution must feature either non-zero radial plastic velocity or non-zero radial magnetic field at the surface in order to be conservative. The former quantity is always zero by assumption in our work, and the latter is also zero for the initial closed-loop configurations. However, a non-zero radial field develops by the time the plastic regime is entered, so we can expect the flow not to dissipate magnetic energy. The same is also true (without caveats) for Hall drift (Goldreich & Reisenegger, 1992). Despite these two processes being conservative, they may still affect the overall magnetic field decay due to the Ohmic term, so we check this in Fig. 6. It is well known that the Hall effect can generate smaller-scale structures that are more readily dissipated by Ohmic decay, and so are indirectly responsible for magnetic-energy loss. The addition of a highly viscous plastic flow results in a mildly faster dissipation of magnetic energy, presumably numerical in origin, since the evolution is now more complex. Interestingly though, if a lower plastic viscosity is chosen (P8), the magnetic field energy decays slower than in the pure-Hall case. We believe this is because plastic flow acts to iron out some of the smaller-scale structures generated by Hall drift, and thus makes Ohmic decay less effective.

5 Discussion

In this paper we believe we have performed the most quantitative, self-consistent simulations to date of the magneto-plastic dynamics of a neutron star’s crust beyond the elastic yield limit. Our main conclusions are the following:

(1) Magnetic fields in the inner crust that are below $10^{15}$ G everywhere barely lead to any plastic flow, as the deviations of the Maxwell stress tensor from its initial state (at which the crust was relaxed) can only ever be marginally above the yield stress.

(2) Broadly speaking, plastic flow opposes Hall evolution, slowing it down and limiting the development of sharp gradients in the magnetic field. A very low plastic viscosity practically chokes Hall evolution completely, and rather fast plastic velocities develop – up to hundreds of centimetres per year.

(3) Increasing the magnetic-field strength, for a given (moderate) value of plastic viscosity does not completely halt the Hall effect, but merely slows it down. Since plastic velocities should scale like $B^{2}$, as opposed to the linear $B$-scaling for Hall velocities, an extrapolation of our results suggests that the evolution of magnetic fields between $5\times 10^{15}-10^{16}$ G will become dominated by the plastic-flow term.

(4) Once the crust fails, we find that the plastic flow could last for $500-2000$ years for our domain size, as long as we account for the global decrease in magnetic energy in calculating the crust’s reference relaxed state. Afterwards it reaches some modified version of a ‘Hall attractor’ state.

Our work has many interesting repercussions for the physics of highly-magnetised neutron star crusts, but before discussing these we need to be aware of various caveats. These are mostly related to the local nature of our simulations, where a small slab of crust is modelled under the assumption of plane-parallel geometry. Firstly, the effectiveness of plastic flow in counteracting Hall drift may be limited by the fact that we place more assumptions on $\boldsymbol{v}_{\mathrm{pl}}$ than on ${\boldsymbol{v}}_{e}$: in particular, we assume the former velocity to be incompressible, to have no $z$-component and to have a $y$-component which is obtained by an averaging procedure. Secondly, the small size of our domain means the characteristic lengthscale of the magnetic field is far shorter than the $1$-km scale one might expect for a global simulation. As a result all of our evolutions are considerably faster than would be expected for a full-crust evolution.

By looking at the results of our evolutions and whether they allow for the expected range of magnetar activity, we can finally place some physical bounds on the viscosity of the crust’s plastic phase. Let us demand that surface motions are sufficiently fast to implant a substantial, persistent twist into the magnetosphere, as required by the corona model of Beloborodov & Thompson (2007); Beloborodov (2009); this is only possible for a viscosity $\nu\lesssim 10^{37}\ \textrm{g cm}^{-1}\textrm{s}^{-1}$. Equally, features related to Hall drift tend to be washed out for $\nu\lesssim 10^{36}\ \textrm{g cm}^{-1}\textrm{s}^{-1}$ – were the viscosity to be so low, we would be led to the undesirable need to dismiss all the literature associating features of magnetar evolution and emission with Hall drift. The noniconoclastic conclusion, avoiding conflict with previous work, is that $\nu\sim 10^{36}-10^{37}\ \textrm{g cm}^{-1}\textrm{s}^{-1}$.

We have shown on theoretical grounds, and from our simulations, that a magneto-plastic evolution generically results in substantial surface motions of the kind that would produce an azimuthal twist in the magnetosphere. This is a non-trivial result, as the flow could have been dominantly further into the crust or only weakly azimuthal. Dropping this assumption, we anticipate at least one qualitatively new feature: with $v_{\rm pl}^{y}$ less restricted, surface motions involving a circulation of plastically-deformed crustal matter. These will immediately render the magnetosphere non-axisymmetric, with a characteristic azimuthal structure dictated by the relative size of $v_{\rm pl}^{x}$ and $v_{\rm pl}^{y}$.

Our simulations, and indeed the consensus about neutron-star crustal failure, relies on the macroscopic plastic flow being a local effect, even though it is seen to be collective on microscopic scales (Chugunov & Horowitz, 2010; Hoffman & Heyl, 2012). This allows for local twisting of magnetospheric field lines and the generation of a localised corona, and larger bursts or flares can be attributed to overtwisting and the resultant violent reconnection (Lyutikov, 2003; Parfrey et al., 2013; Akgün et al., 2018). There is perhaps one dissenting voice to this. Thompson et al. (2017) consider a toy model for core-field evolution and based on its results argue that crustal failure is truly global: by imposing a large-scale stress at the crust-core boundary, a localised yielding event in a small region can propagate across the whole crust. Their simulations show some interesting features that resemble observed magnetar activity, like the presence of lower-energy ‘aftershocks’ after a giant flare. The model makes, however, some inevitable simplifications, which impede a quantitative comparison with neutron-star observables. In particular, their core field is not a hydromagnetic equilibrium solution, but rather modelled by an elastic spring, and their evolution – changing the spring constant – does not seem reflective of any known mechanism for core-field evolution.

We do, however, believe a viable picture for the giant flares of magnetars can be produced from the ideas of Beloborodov & Levin (2014) and Thompson et al. (2017), and using some insight from our simulations. If the plastic flow were not strictly conservative but instead involved some viscous heating, the dynamics would be modified: excessive heating would probably prevent significant surface motion, but mild heating would reduce the local viscosity of the plastic phase and accelerate the evolution. A sustained episode of accelerated plastic flow would implant a magnetospheric twist growing on a timescale far shorter than that of exterior dissipation, and eventually might trigger an overtwisting instability (Mikic & Linker, 1994; Parfrey et al., 2013). Since the requisite balance between the plastic heating and flow rate may be rather delicate, it would not surprising to find that such energetic events are very rare. Certainly, some mechanism is needed to ensure a localised failure – like that we simulate here – results in motions sufficiently widespread to implant a large-scale twist in the magnetosphere.

Electron MHD is, mathematically, a clean model for neutron-star crustal evolution. The range of input parameters is mostly quite constrained, results are insensitive to initial conditions, and the natural next step towards more realistic evolutions has typically been quite clear. Many of the considerable challenges have instead been numerical: for example, dealing with the formation of shocks due to Hall drift (Viganò et al., 2012), or moving from two to three spatial dimensions (Wood & Hollerbach, 2015; Gourgouliatos et al., 2016). Plastic flow introduces significant new uncertainties into modelling, and complicates the association of observables with the crust’s underlying physics. We have limited understanding of how and when the crust fails and enters the plastic phase, and of the material response of the crust to high stress. The initial conditions and the crust’s seismic history become important. On the other hand, this new physics is rich and interesting, and provides us with a possibility of using observations to constrain some of these unknown mechanical properties of the neutron-star crust.

Acknowledgements

We thank Chris Thompson and the anonymous referee for comments which helped to clarify the presentation of this work. SKL acknowledges support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 665778, via fellowship UMO-2016/21/P/ST9/03689 of the National Science Centre, Poland. Both authors thank the PHAROS COST Action (CA16214) for partial support.

This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.