Global Crustal Dynamics of Magnetars in Relation to Their Bright X-Ray Outbursts

We begin by listing several open questions relating to magnetar activity, describing our approach to them, and offering the conclusions we have drawn.

The basic picture developed here views the magnetar crust as being composed of many interacting elastic units, which experience intermittent plastic creep over a wide range of rates. Maxwell stresses imparted to the crust at its lower boundary source nonlocal solid stresses. The creep rate in a given patch is determined in part by the collective stress imparted by its neighbors, and it is regulated by conductive cooling, neutrino emission, and the shearing of an embedded magnetic field.

Can the crust experience runaway deformation to dynamical rates, comparable to (shear wave speed)/(crustal scale height)? Or are its deformations restricted to relatively slow plastic creep? We show that a slow increase in applied stress, sourced by distant Maxwell stresses, can drive rapid creep within a small part of the crust. The localized rapid creep does not feed back immediately on the imposed global stress. The size of this zone, in which inertial forces start to counterbalance the Lorentz force and elastic force, depends on the growth rate of the applied stress. Locally the dense solid crust of a neutron star can momentarily behave like a fluid.

What triggers a growth in creep rate? Molecular dynamics simulations of dense Coulomb solids (Chugunov & Horowitz 2010; Hoffman & Heyl 2012) indicate that the creep rate $\dot{\varepsilon }$ is much more sensitive to changes in stress σ than it is to changes in temperature T (Figure 1). This opens up the possibility of a chain reaction within the solid strain pattern, in which creep in one part of the magnetar crust changes the stress in another part, thereby triggering even faster creep.

Zoom In Zoom Out Reset image size

What is the time profile of the dissipation? The net plastic dissipation rate in the crust peaks at a wide range of values, depending on the details of how stresses are configured. In some cases, we find energy release profiles similar to those observed in giant magnetar flares: very fast (millisecond) rise, ∼0.1 s duration, followed by an extended power-law decay of plastic creep.

How does fast creep couple to the magnetosphere? Magnetic twist or shear can be directly ejected on a vertical shear wave timescale, consistent with some early theoretical ideas (Thompson & Duncan 1995). This happens preferentially in zones where plastic flow has previously deposited a horizontal magnetic field. The efficiency of this process, as measured by the escaping Poynting flux, is greatly enhanced by the presence of a background magnetic twist extending into the magnetosphere. The energy that powers an SGR burst is not primarily stored in an elastic wave, which only couples relatively slowly to magnetospheric modes (Blaes et al. 1989; Link 2014). The ejected magnetic displacement may or may not oscillate in the magnetosphere, depending on the length of the excited field lines.

What is the origin of the QPO activity seen in large magnetar flares? Stresses are redistributed throughout the magnetar crust by elastic forces, modified by an embedded magnetic field (Duncan 1998; Piro 2005). Changes in the stress balance in one part of the crust will have a larger influence on a neighboring part than changes in temperature, both because of the greater sensitivity of creep rate and because elastic waves propagate much more rapidly than temperature fluctuations. We find that modulation of the local creep rate by a global elastic mode can feed back positively on the global mode.

How does a temperature increase influence subsequent yielding? As has been explored in models of terrestrial faults (Turcotte & Schubert 2002), in situ plastic heating raises the rate of plastic creep; it also enhances the duration and extent of the creep, by reducing the equilibrium stress. These effects have been quantified by Chugunov & Horowitz (2010) for neutron star crusts. They are especially important in layers with a strongly temperature-dependent specific heat, e.g., where the temperature is well below the Debye temperature and the dripped neutrons are strongly paired.

Rapid temperature feedback in a magnetar crust does not automatically follow from transport of the magnetic field. Transport in the crust (by electron drift) or in the core (by electron–proton drift) is generally slow compared with the conduction of heat in the crust. As a result, heat deposited by plastic flow may be conducted away faster than creep can grow; growth of the Maxwell stress may be compensated by the solid response.

Beloborodov & Levin (2014) posit that creep has already reached high enough rates for the heat generated to be trapped in a part of the magnetar crust. In a one-dimensional calculation, they consider the conduction of heat away from the yielding zone. This triggers the creation of a propagating front and a continuing release of both elastic and Maxwell stresses. Because heat transfer is needed to keep the process going, energy is released on a timescale intermediate between the local dynamical time and the conduction time. A one-dimensional simulation of Hall drift, in which the vertical electric current is assumed to vanish, produces small intermittent episodes of thermal feedback in the upper crust (Li et al. 2016).

On the other hand, creep triggered by a growing global stress does not depend on sharp temperature gradients. It is also worth emphasizing that a thermal runaway feeds off the local elastic stress, not the local Maxwell stress. The yielding behavior uncovered here cannot be understood in terms of a local "burning" process: the horizontal (x–y) component of the embedded Maxwell stress typically increases in response to rapid plastic creep: in other words, the embedded field can act as a "fire retardant." We also find that fast dynamical creep develops more easily if the initial temperature lies below ∼5 × 10⁸ K. This means that the crust of a neutron star has a much "jumpier" response to magnetic stresses at temperatures well below the melt temperature, supporting the conjecture of Ruderman (1991).

What is the role of electron thermal conduction in regulating plastic creep? We find that thermal conduction generally reduces the creep rate, and in some circumstances it can prevent runaway creep. Equilibrium solutions for a plastic deformation "spot" (zero-dimensional) and "fault" (one-dimensional) can be constructed in which the heat generated by plastic deformation is conducted away and/or radiated to neutrinos.

What is the relative importance of core and crustal magnetic fields to "breaking" the crust? Having fixed the magnitude of the magnetic field, one obtains a larger solid stress from an imbalance in the x–z and y–z components of the Maxwell stress at the crust–core boundary, as compared with the stress locally imparted by magnetic irregularities within the crust (Thompson & Duncan 2001). The ratio of the two is R_⋆/δR_c, where R_⋆ is the stellar radius and δR_c the crustal thickness. This Maxwell stress imbalance has a nonlocal effect, increasing the solid stress at distant points in the crust (Figure 2).

Zoom In Zoom Out Reset image size

How important is the role of hydromagnetic instability in driving the largest magnetar flares? Solid stresses in the crust can resist only a limited imbalance between hydrostatic, gravitational, and elastic stresses. We find that giant flare energies larger than ∼10⁴⁵ erg depend on the onset of a hydromagnetic instability in the magnetar core. We also show that large stress imbalances may develop spontaneously in the core on timescales as short as a day. This is fast enough for the heat deposited in crustal fault zones to remain trapped. A similar effect could be produced by current-driven instabilities in the magnetosphere.

Can fault-like features form in the magnetar crust? We find that if the crust is stressed from below in a slightly anisotropic manner (corresponding, e.g., to a large-scale winding of the core magnetic field), then concentrated zones of strong shear develop naturally (Figure 3). Adding a horizontal magnetic field to the crust does not suppress this effect, but it can modify the width of the shear zone and the peak strain rate within it.

Zoom In Zoom Out Reset image size

The formation of open cracks is inhibited by the extreme ratio of hydrostatic pressure to shear modulus in the crust of a neutron star (Jones 2003). An earlier analysis by Levin & Lyutikov (2012) started by positing the spontaneous formation of a crack near a surface of maximal stress, which was assumed to run perpendicular to the magnetic field, and followed the acceleration of the magnetized solid near the crack. Yielding was restricted to a narrow layer whose thickness is set by ohmic diffusion, and the energy release was found to be minimal.

The fault-like features presented here are driven by continued stress growth, which, being nonlocally sourced, does not stop as creep begins to pick up. Much of the expansion of creep is associated with a nonlocal redistribution of stress in the neighborhood of plastically deforming zones, which is an inherently two- or three-dimensional effect. As a result, plastic flow develops over a much wider band than the ohmic diffusion length. These processes clearly have a strong influence on the spectral signature of localized heating.

In a thin-shell geometry, the shell thickness defines a characteristic scale for extended linear bands of peak stress, below which the strain field is no longer essentially two-dimensional. We evolve the global solid stress pattern by dividing the crust into $\sim 4\pi {({R}_{\star }/\delta {R}_{c})}^{2}\sim {10}^{4}$ elastic units.

Why are all magnetars not bursting SGRs, and why are SGRs not bursting all the time, given the extreme nonlinearity and temperature sensitivity of the creep rate? We find that shearing of an embedded magnetic field feeds back more strongly on peak creep rate than on the net energy dissipated. The Coulomb solid attains a state of equilibrium creep (at rates intermediate between the dynamical rate and the inverse age of the star) if the magnetic field transverse to the direction of shear reaches a substantial strength, greater than 10¹⁵ G at the base of the crust (corresponding to a Maxwell stress ≳0.1 the shear modulus). As this embedded field weakens, the crust becomes more sensitive to slow changes in applied stress. We also find that a tangled (as opposed to laminar) field is most effective at limiting runaway large-scale creep. Hall drift is a plausible source of small-scale magnetic irregularities (Goldreich & Reisenegger 1992), and in this way it may actually limit SGR burst activity.

In contrast with the approach taken here, a number of calculations of magnetic field evolution assume that the magnetic field is anchored entirely in the magnetar crust, setting aside any effects of a poloidal field threading crust/core boundary (e.g., Perna & Pons 2011; Beloborodov & Levin 2014; Gourgouliatos et al. 2016, and references therein). This is done, in part, for reasons of computational facility: one can then neglect the response of the liquid core to changing Maxwell stresses.

It is therefore worth asking whether a purely crustal magnetic field is a realistic configuration in a neutron star. Considerations of the origin of this field suggest that it is not. The progenitor star experiences a series of hydrodynamic instabilities that successively modify the magnetic field, starting with main-sequence core convection and growing in power to the most extreme core collapse phase (Thompson & Duncan 1993). Post-collapse accretion following the brief stalling of a bounce shock deposits some 3–10 times the mass of the eventual solid neutron star crust. This allows a transient activation of the magnetorotational instability in a rarefied mantle (Akiyama et al. 2003), which is followed by more persistent linear winding of the magnetic field as the mantle collapses. Rapid neutrino heating drives buoyant motions of the toroidal magnetic field if the seed field exceeds a critical value (Thompson & Murray 2001), thereby potentially distinguishing magnetars from neutron stars with a combination of slower rotation and weaker seed fields.

The dilution of the core magnetic field by entrainment with outward-drifting superfluid vortices (Ruderman et al. 1998) is impeded in comparison with ordinary radio pulsars by a relatively low ratio of rotational to magnetic energy. Given that the core superfluid transition takes place at age ∼300 yr in a more weakly magnetized neutron star such as the Cas A remnant (Ho et al. 2015), one sees that a magnetar must reach a spin period >1 s before superfluid vortices appear in its core. Heating of the cores of the most active magnetars can delay this transition further; indeed, an active lifetime of ∼10³–10⁴ yr was used to infer a core superfluid transition temperature of ≲6 × 10⁸ K (Arras et al. 2004).

We therefore adopt the following simplified magnetic field configuration: the field is strongly tilted or twisted in the core but becomes essentially vertical in the crust. Of course, a hybrid configuration is plausibly attained in magnetars, with the axis of the magnetic field winding tilted with respect to the local direction of gravity. But the chosen configuration is adequate to demonstrate the effects we are interested in.

Empirically there appears to be a strong correlation between the detection of dense episodes of short burst activity from magnetars and their ability to generate a giant flare. This suggests that the same structures responsible for the giant flares also give rise to the shorter bursts. The giant flares, in turn, appear to depend on large-scale unbalanced stresses in the neutron star, which partly motivates our approach that starts with a core magnetic field making a transition out of magnetostatic equilibrium.

The plan of this paper is as follows. Section 2 analyzes the plastic response of a small patch of Coulomb solid to an externally imposed stress.

Section 3 sets up the problem of plastic flow and elastic response in two dimensions, in response to an inhomogeneous Maxwell stress that is applied at the base of the crust.

Section 4 analyzes the results of global time-dependent calculations of such a stressed model crust. Section 5 presents outburst light curves and shows how the energy dissipation is distributed between magnetic, plastic, and elastic components and how the peak of dissipation circulates around the crust.

In Section 6 we consider the vertical redistribution of stresses in a solid with a strong density stratification, appropriate to a neutron star crust. Fast ejection of magnetic twist is encountered when a solid crust experiences runaway creep in response to an externally applied stress. This provides, we believe, a promising foundation for understanding the short (∼0.1 s) SGR bursts.

The slow damping of free vertical Alfvénic oscillations in liquified patches of the crust is described in Section 7. This suggests that the high-frequency QPOs observed in magnetar outbursts (Strohmayer & Watts 2006) may involve the vertical excitation of the magnetic field in melted parts of the crust (e.g., within heated fault zones).

The interaction between the local Maxwell stresses and global elastic oscillations is addressed in Section 8: we show that the global mode is overstable if the localized stresses have a "patch"-like geometry.

The concluding Section 9 makes a comparison with independent theoretical approaches and summarizes the implications of our work for the physics of magnetars. The two Appendices describe resolution tests of the two-dimensional crustal model, and present a simple example of a plastic spot.

We consider the solid stress to be the sum of local and global contributions,

$$\begin{eqnarray}&&\sigma ={\sigma }_{l}+{\sigma }_{g}={\varepsilon }_{\mathrm{el},l}\,\mu +{\sigma }_{g}.\end{eqnarray} \tag{ 4 }$$

The local stress σ_l compensates a local Maxwell stress that grows in magnitude in response to plastic shear. The applied global stress σ_g is taken to grow on a timescale t_growth. It is sourced by Maxwell stresses at a distance L_x ≫ δR_c from the plastic zone and varies slowly within the plastically deforming patch.

We work in a frame where the plastic flow velocity v_∥ vanishes at a given point in the solid, so that

$$\begin{eqnarray}&&{v}_{\parallel }=({\dot{\varepsilon }}_{\mathrm{pl}}+{\dot{\varepsilon }}_{\mathrm{el},l}){x}_{\perp }.\end{eqnarray} \tag{ 5 }$$

Here $\perp \text{and}\,\parallel$ are label coordinates perpendicular and parallel to the direction of plastic flow. Here there is an additional contribution to the solid displacement from the elastic response to an imbalance between σ_l and ${B}_{\perp }{B}_{\parallel }/4\pi$, which develops as the result of plastic creep. The system moves through a series of equilibrium states in which ${\partial }_{\perp }({B}_{\perp }{B}_{\parallel }/4\pi +{\sigma }_{l}+{\sigma }_{g})\,={\partial }_{\perp }({B}_{\perp }{B}_{\parallel }/4\pi +{\sigma }_{l})=0.$ Hence,

$$\begin{eqnarray}&&\displaystyle \frac{{B}_{\perp }{\partial }_{t}{B}_{\parallel }}{4\pi }+{\partial }_{t}{\sigma }_{l}=\displaystyle \frac{{B}_{\perp }{\partial }_{t}{B}_{\parallel }}{4\pi }+{\dot{\varepsilon }}_{\mathrm{el},l}\,\mu =0,\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&{\partial }_{t}{B}_{\parallel }={\partial }_{\perp }({v}_{\parallel }{B}_{\perp })=({\dot{\varepsilon }}_{\mathrm{pl}}+{\dot{\varepsilon }}_{\mathrm{el},l}){B}_{\perp }.\end{eqnarray} \tag{ 7 }$$

Combining these equations, one sees the local stress changes in a negative manner with respect to the global stress, so as to partly compensate its growth,

$$\begin{eqnarray}{\partial }_{t}{\sigma }_{l} & = & -\,\displaystyle \frac{{B}_{\perp }^{2}/4\pi \mu }{1+{B}_{\perp }^{2}/4\pi \mu }{\dot{\varepsilon }}_{\mathrm{pl}}\,\mu \equiv -{f}_{\mathrm{local}}\cdot {\dot{\varepsilon }}_{\mathrm{pl}}\,\mu ;\\ {\partial }_{t}{\sigma }_{g} & = & \displaystyle \frac{{\sigma }_{g,0}}{{t}_{\mathrm{growth}}}-(1-{f}_{\mathrm{local}})\displaystyle \frac{{\dot{\varepsilon }}_{\mathrm{pl}}{\rm{\Delta }}}{{L}_{x}}\mu .\end{eqnarray} \tag{ 8 }$$

The equation for σ_g includes a second term representing the decrease in large-scale elastic stress caused by the net displacement ${v}_{\parallel }\sim \,({\dot{\varepsilon }}_{\mathrm{pl}}+{\dot{\varepsilon }}_{\mathrm{el},l}){\rm{\Delta }}$ $=\,(1-{f}_{\mathrm{local}}){\dot{\varepsilon }}_{\mathrm{pl}}{\rm{\Delta }}$ across the plastic patch of width Δ.

The energy equation is evolved at density 10¹⁴ g cm⁻³, where the global stress is concentrated (Section 3),

$$\begin{eqnarray}&&{C}_{V}{\partial }_{t}T={\dot{\varepsilon }}_{\mathrm{pl}}(\sigma ,T)\sigma -{{\ell }}^{-2}\,{K}_{\mathrm{cond}}(T-{T}_{\mathrm{base}})-{\dot{Q}}_{\nu }(T,B).\end{eqnarray} \tag{ 9 }$$

In Equation (9), C_V is the specific heat, which is approximated by the contribution of the ion lattice (Chabrier 1993), with the effective ion mass raised by the entrainment of a fraction 0.3 of the superfluid neutrons (Onsi et al. 2008). Then C_V ∼ T³ at temperatures below the Debye temperature, meaning that the creep rate is more sensitive to heat input from in situ dissipation at T ∼ 10⁸ K (appropriate for transient magnetars in their low states) as opposed to (5–10) × 10⁸ K (appropriate for persistently bright magnetars). The contribution of dripped neutrons to the specific heat varies rapidly with density: it is strongly suppressed by superfluidity where the pairing temperature is high compared with T. Runaway creep is therefore concentrated within the superfluid part of the lower crust.

Two terms are added to the right-hand side of Equation (9), representing thermal conduction out of the "spot" over a characteristic size (temperature gradient scale) ℓ, and neutrino cooling. The thermal conductivity varies weakly with temperature near 10¹⁴ g cm⁻³, K_cond ≃ 1 × 10²⁰ erg (cm s K)⁻¹ (Potekhin et al. 2015). Neutrino emission is dominated by thermal bremsstrahlung and synchrotron radiation by relativistic electrons and calculated using formulae given in Yakovlev et al. (2001). The conductive cooling length is taken to be ℓ = 0.3 km, approximately the vertical density scale height at the base of the crust.

Here ohmic effects can be neglected: we are studying smooth plastic flow on macroscopic scales, driven by a continued growth in applied stress. Sharper gradients in the strain field and magnetic field are encountered if a zero-stress discontinuity is introduced near the stress maximum, and the solid evolves passively thereafter (Levin & Lyutikov 2012). If instead the global stress is allowed to grow, then strengthening the transverse magnetic field does not limit the energy deposition in a shear band, but instead pushes it to larger scales (Section 2.2). The regulating effect of horizontal thermal conduction perpendicular to the fault is addressed in the one-dimensional model of Section 2.3.

A simple numerical experiment evolves the "spot" starting from an initial temperature T_base, with the initial total stress ${\sigma }_{0}={\sigma }_{g,0}+{\sigma }_{l,0}$ corresponding to a very long creep time (>1 Myr) at this temperature. The transverse magnetic field is normalized by the parameter f_local in Equation (8), and we consider a range of t_growth for the global stress ranging from about a day up to 10⁵ yr (longer than a typical magnetar lifetime).

Figure 4 shows the evolution of creep rate, total stress, and temperature in the spot. As the applied stress ramps up more quickly, one finds, not surprisingly, a faster and sharper response of the plastic flow in the spot. The initial response, over a timescale (0.01–0.1)t_growth, is driven by the nonlinear stress dependence. Thereafter, there is a brief interval where thermal feedback takes over, corresponding to the sharpest rise in T. The creep rate ${\dot{\varepsilon }}_{\mathrm{pl}}$ peaks and then subsides as the transverse magnetic field is sheared out, which allows the local Maxwell stress to cancel out further increase in the global stress. Most of the drop in total stress occurs after peak creep is reached, over a timescale (0.1–0.3)t_growth. Smaller initial temperatures correspond to higher peak creep rates, a consequence of the T³ scaling of the specific heat.

Zoom In Zoom Out Reset image size

The stress, temperature, and creep rate all eventually reach a plateau. Dynamically this can arise in two distinct ways. There are two negative terms on the right-hand sides of Equation (8) arising from (i) field line stretching and (ii) direct elastic feedback of creep on σ_g. The first dominates when

$$\begin{eqnarray}&&\displaystyle \frac{{B}_{\perp }^{2}}{4\pi \mu }\gtrsim \displaystyle \frac{{\rm{\Delta }}}{{L}_{x}},\end{eqnarray} \tag{ 10 }$$

and one has

$$\begin{eqnarray}&&{\dot{\varepsilon }}_{\mathrm{pl}}\simeq \left(1+\displaystyle \frac{4\pi \mu }{{B}_{\perp }^{2}}\right)\displaystyle \frac{{\sigma }_{}}{\mu {t}_{\mathrm{growth}}}.\end{eqnarray} \tag{ 11 }$$

The key feature here is that ${\dot{\varepsilon }}_{\mathrm{pl}}\propto {B}_{\perp }^{-2}$ unless the transverse magnetic field is extremely strong. This gives ${\dot{\varepsilon }}_{\mathrm{pl}}{t}_{\mathrm{growth}}\gg 1$, and there is significant amplification of the creep rate compared with the "seed" rate σ/μt_growth.

The plateaus also represent a balance between in situ plastic heating and neutrino cooling for all except the longest t_growth and the lower values of T_base. One observes some bimodality in final temperature and stress in the latter cases, representing a switch between a warmer, neutrino-cooled state and a cooler state in which plastic heating is balanced by conductive cooling. Large swings in creep rate are observed in a few cases where t_growth is comparable to 1 yr, the conductive cooling time across ℓ ∼ 0.3 km.

It is also instructive to compare the behavior of this zero-dimensional model with the one-dimensional thermal Hall front proposed by Beloborodov & Levin (2014); this gives a clearer sense of how the different physical effects interact with each other. As is noted by these authors, runaway damping of the solid stress can readily be found by initializing the stress to a high enough value that creep is faster than thermal conduction. But this immediately raises the question of how the minimal creep is established in the first place. Instead of introducing a nonlocal source term for the solid stress (as done here), Beloborodov & Levin (2014) consider an increase in temperature driven by thermal conduction from a neighboring yielding layer. The front that forms is much slower than an elastic wave and involves sharp gradients in temperature and Lorentz force. Another distinction is that the Maxwell stress decreases in magnitude behind such a thermal-magnetic front, instead of increasing so as to cancel growth in the global stress.

Returning to the present model, one infers from Figure 4 that there is a critical growth rate ${t}_{\mathrm{growth}}^{-1}$ for the applied stress above which ${\dot{\varepsilon }}_{\mathrm{pl}}$ approaches the dynamical value V_μ/ℓ, where inertial forces must be taken into account and the flow rate must saturate. Here V_μ = (μ/ρ)^1/2 is the shear wave speed in an unmagnetized Coulomb solid. Figure 5 shows that the crust is most susceptible to runaway creep at lower temperatures: slower growth of the applied stress (higher t_growth) is required. Alternatively, strengthening the transverse magnetic field makes runaway creep more difficult. Indeed, at fixed T_base one can consider there to be a threshold value of f_local leading to explosive creep growth. For example, at T_base = 5 × 10⁸ K the critical t_growth drops from 10¹¹ s (comparable to the magnetar lifetime) down to ∼1 s as f_local rises from 0.05 to 0.1.

Zoom In Zoom Out Reset image size

The preceding zero-dimensional creep model did not specify the width Δ of the plastic layer. Consider now the case where the applied global stress has a broad maximum at x = 0 and falls away slowly on some length scale L_x ≫ δR_c to $| x| \gt 0$ (Figure 2). As the amplitude of the global stress σ_g rises, the width of the layer that achieves a certain stress grows.

Here we must distinguish between the width Δ_eq of the layer that achieves the equilibrium creep rate (11), where growth in the applied stress is canceled off, and the (instantaneous) width Δ_T of the layer experiencing thermal runaway. The second is much narrower than the first, except when the applied stress grows on subsecond timescales. This means that Δ_eq provides a first estimate of the width of the dissipative zones that persist for months to decades (or longer) in transient magnetars and the AXPs. On the other hand, we observe rapid and extensive growth of fault-like features during giant outbursts in the global two-dimensional yielding calculations described in Section 4. In this situation t_growth is shorter than ∼0.1 s, implying that Δ_eq and Δ_T are comparable, as we now show.

If the transverse magnetic field B_⊥ is very small, then Equation (10) is only satisfied as long as ${\rm{\Delta }}\lesssim {{\rm{\Delta }}}_{\mathrm{eq}}\,={L}_{x}({B}_{\perp }^{2}/4\pi \bar{\mu })\ll \delta {R}_{c}$. Beyond that stage of growth, the long-term, equilibrium creep rate becomes large enough to compensate any further growth of σ_g, stress growth stalls outside the fault, and Δ_eq saturates.

On the other hand, if ${B}_{\perp }^{2}/4\pi \gtrsim (\delta {R}_{c}/{L}_{x})\mu$, then magnetic shearing sufficiently limits ${\dot{\varepsilon }}_{\mathrm{pl}}$ that Δ_eq can reach (or exceed) δR_c as the stress builds up. In a neutron star crust with δR_c ∼ 0.03R_⋆ this corresponds to B_⊥ ≳ 0.6 × 10¹⁵ G at a density of ∼10¹⁴ g cm⁻³, well within the magnetar range. Even weaker fields embedded in the lower-density parts of the crust will facilitate faults of width Δ_eq ∼ δR_c.

Now consider Δ_T, which we define as the width of the layer experiencing peak strain rate ${\dot{\varepsilon }}_{\mathrm{pl},\max }$. The latter quantity depends in a highly nonlinear manner on initial temperature, t_growth, and B_⊥, so we simply normalize it to the dynamical rate V_μ/L_x ∼ 10² s⁻¹. As is seen in the top left panel of Figure 4, peak creep involves a rapid increase $\delta ({B}_{\perp }{B}_{\parallel }/4\pi )\sim {10}^{-2}\mu$ in the local Maxwell stress (in the case of low t_growth and low initial temperature T_base). Taking f_local ≪ 1, the duration of peak creep is seen from Equation (8) to be

$$\begin{eqnarray}&&\delta {t}_{\max }\sim \displaystyle \frac{\delta ({B}_{\perp }{B}_{\parallel }/4\pi )}{{\dot{\varepsilon }}_{\mathrm{pl},\max }{f}_{\mathrm{local}}\mu }\sim 2\times {10}^{-3}\,{B}_{\perp ,15}^{-2}{\left(\displaystyle \frac{{\dot{\varepsilon }}_{\mathrm{pl},\max }}{{V}_{\mu }/{L}_{x}}\right)}^{-1}\ {\rm{s}},\end{eqnarray} \tag{ 12 }$$

where B_⊥ = B_⊥,15 × 10¹⁵ G. The width of the layer experiencing peak creep is

$$\begin{eqnarray}&&{{\rm{\Delta }}}_{T}\sim \displaystyle \frac{\delta {t}_{\max }}{{t}_{\mathrm{growth}}}{L}_{x}=0.2\,{B}_{\perp ,15}^{-2}{\left(\displaystyle \frac{{\dot{\varepsilon }}_{\mathrm{pl},\max }}{{V}_{\mu }/{L}_{x}}\right)}^{-1}{\left(\displaystyle \frac{{t}_{\mathrm{growth}}}{0.1{\rm{s}}}\right)}^{-1}\ \mathrm{km},\end{eqnarray} \tag{ 13 }$$

which is similar to the crustal thickness δR_c.

Widening of a fault is also caused by the redistribution of stress into the surrounding solid. This is a two-dimensional effect, which is associated with curvature of the fault and also with the expansion of stress around the tip of a growing fault. It cannot be captured by the (effectively one-dimensional) Cartesian geometry adopted here. A nice example of nonlocal dissipation around a curved fault is shown in Figure 22.

We now turn to a steady, one-dimensional model of plastic creep in a fault-like geometry. The maximum of temperature and stress sits in a vertical plane situated (say) at horizontal coordinate x = 0. The injection of heat by plastic creep is balanced by volumetric neutrino cooling and thermal conduction away from the stress maximum. Similar solutions were constructed some time ago by Schubert & Yuen (1978), with a goal of explaining the formation of faults in Earth's lithosphere.

The vertical scale height δR_c provides a characteristic scale below which the temperature gradient is mainly horizontal, and vertical conductive losses can be neglected, so that T = T(x). We first consider a weakly magnetized fault (i.e., neglect the calming effect of an embedded magnetic field within the creep layer).

We construct steady profiles of temperature, strain rate, and heat flux, with a range of background temperature T(δR_c) = T_base, from the coupled equations

$$\begin{eqnarray}{F}_{\mathrm{cond}} & = & -{K}_{\mathrm{cond}}\displaystyle \frac{{{dT}}_{\mathrm{cond}}}{{dx}};\\ \sigma {\dot{\varepsilon }}_{\mathrm{pl}}(T,\sigma ) & = & -\displaystyle \frac{{{dF}}_{\mathrm{cond}}}{{dx}}-{\dot{Q}}_{\nu }(T,B).\end{eqnarray} \tag{ 14 }$$

Each value of T_base defines a sequence of solutions with different peak temperatures T_mid ≡ T(0) at the middle of the fault and different total stress σ = σ_g = const. The conductive energy flux F_cond vanishes at x = 0, corresponding to a solution that is symmetric about the temperature maximum: T(−x) = T(x) and F_cond(−x) = −F_cond(x).

The global stress relaxes on a timescale

$$\begin{eqnarray}&&{t}_{\mathrm{creep}}\sim \displaystyle \frac{{R}_{\star }}{{\rm{\Delta }}{v}_{\parallel }},\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}&&{\rm{\Delta }}{v}_{\parallel }={\int }_{0}^{\delta {R}_{c}}{\dot{\varepsilon }}_{\mathrm{pl}}[T(x),\sigma ]{dx}\end{eqnarray} \tag{ 16 }$$

is the net differential creep speed across the fault. Figure 6 shows the sequence of curves relating t_creep to the magnitude of the stress. Each curve is labeled by T_base, with lower stresses corresponding to higher midplane temperatures. One also sees a tight relation between T_mid and σ, reflecting the extreme sensitivity of ${\dot{\varepsilon }}_{\mathrm{pl}}$ to σ and T.

Zoom In Zoom Out Reset image size

The magnetic field within the magnetar core is slowly reconfigured by ambipolar drift and beta reactions (Haensel et al. 1990; Goldreich & Reisenegger 1992; Pethick 1992) over its active lifetime of t_NS ∼ 10¹¹ s. The feedback of heating on the rate of modified URCA reactions allows the diffusion time to scale with the lifetime of the star while the core neutrons remain in a normal state (Thompson & Duncan 1996). Slow transport of the magnetic field through a fluid star allows the field to evolve into a nonaxisymmetric configuration (Braithwaite & Spruit 2006) that becomes susceptible to a global hydromagnetic instability.

Many details of how this happens cannot yet be captured by global hydromagnetic simulations, in which numerical diffusion is typically rapid in the outer parts of the model star. Nonetheless, simple considerations of continuity suggest that modes of a very low growth rate must develop during a gradual transition away from stability. We make the simplest approximation of treating the core field as a one-dimensional spring, with a spring constant $\propto {\,\omega }_{0}^{2}$. Stable equilibrium corresponds to ω₀ ∼ V_A/R, where V_A is the Alfvén speed within the core.

Now we let the spring constant evolve linearly through zero over the diffusion time ∼t_inst (which initially may be comparable to t_NS),

$$\begin{eqnarray}&&{\omega }_{0}^{2}\sim {\left(\displaystyle \frac{{V}_{{\rm{A}}}}{R}\right)}^{2}\left(1-\displaystyle \frac{t}{{t}_{\mathrm{inst}}}\right).\end{eqnarray} \tag{ 17 }$$

Near t = t_inst, one finds that $| {\omega }_{0}|$ is much smaller than the spin frequency Ω (typically ∼1 s⁻¹ for an X-ray-bright magnetar), meaning that the growth of a hydromagnetic instability is limited by the Coriolis force. Balancing this force against internal magnetic stresses gives ${\rm{\Omega }}| v| \sim | {\omega }_{0}^{2}| {R}_{\star }$. The hydromagnetic displacement at speed v can then be converted to a growth rate

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{eff}}\sim -\displaystyle \frac{{\omega }_{0}^{2}}{{\rm{\Omega }}}.\end{eqnarray} \tag{ 18 }$$

Instability occurs where ${\omega }_{0}^{2}\lt 0$, but in a fluid star it can develop only when Γ_effΔt > 1, where ${\rm{\Delta }}t=t-{t}_{\mathrm{inst}}$. This gives a minimum growth time

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\mathrm{eff}}^{-1}\gtrsim {(\frac{R}{{V}_{{\rm{A}}}})}^{}\left({\rm{\Omega }}\,{t}_{\mathrm{inst}}^{1/2}\right),\end{eqnarray} \tag{ 19 }$$

which is much shorter than the neutron star lifetime, but still longer than the Alfvén crossing time R/V_A ∼ 0.1 s. Taking t_inst ∼ t_NS ∼ 10¹¹ s gives ${{\rm{\Gamma }}}_{\mathrm{eff}}^{-1}\sim 3\times {10}^{4}\,{\rm{s}}$.

These considerations apply to a fluid star. The next question that arises is whether the crustal elastic stress should force the transition from a decaying to a growing mode away from the vanishing core spring constant (ω₀ = 0). That would not happen if, at each instant in the evolution of the star, the core could be placed in magnetohydrostatic equilibrium, with the vanishing shear component of the Maxwell stress at the base of the crust. But the magnetar will start at some early time t with the shear stress above the yield stress. When the core is otherwise close to magnetohydrostatic equilibrium, the shear Maxwell stress can only relax to the point where the plastic creep rate has dropped to $\dot{\varepsilon }\sim {\varepsilon }_{b}/t$. Here ${\varepsilon }_{b}\equiv \sigma ({\dot{\varepsilon }}_{\mathrm{pl}})/\mu \sim 0.02\mbox{--}0.1$ is the yield strain at creep rate ${\dot{\varepsilon }}_{\mathrm{pl}}$ as described by Equations (1) and (3).

This suggests a numerical experiment in which the creep rate starts at a very low value in the plastic parts of the crust. In the initial state, the lattice strain exceeds the yield strain only in small patches. We allow the Maxwell stress to grow very slowly and then see whether the crust is capable of a much faster response. Only then do we vary the sign of the core spring constant.

We idealize the crust as a planar solid layer that is threaded by a uniform vertical magnetic field B_z, which is tilted by some angle at the lower boundary (Figure 2). The horizontal magnetic field at the lower boundary can be expressed in terms of a solenoidal displacement field ${{\boldsymbol{\xi }}}_{B}$ at the base of the core (depth h_core ∼ R_⋆)

$$\begin{eqnarray}&&{{\boldsymbol{B}}}_{\perp }^{-}=-\displaystyle \frac{{{\boldsymbol{\xi }}}_{B}}{{h}_{\mathrm{core}}}{B}_{z},\end{eqnarray} \tag{ 20 }$$

with ${\partial }_{x}{B}_{x}^{-}+{\partial }_{y}{B}_{y}^{-}={\partial }_{x}{\xi }_{B,x}+{\partial }_{y}{\xi }_{B,y}=0$. We define ${{\boldsymbol{\xi }}}_{B}$ as the net displacement before any compensating horizontal motion of the crust. The applied Maxwell stress drives both elastic and plastic displacements in the crust,

$$\begin{eqnarray}&&{\boldsymbol{\xi }}(x,y)={{\boldsymbol{\xi }}}_{\mathrm{el}}(x,y)+{{\boldsymbol{\xi }}}_{\mathrm{pl}}(x,y),\end{eqnarray} \tag{ 21 }$$

in response to which

$$\begin{eqnarray}&&{{\boldsymbol{B}}}_{\perp }^{-}\to \displaystyle \frac{{\boldsymbol{\xi }}-{{\boldsymbol{\xi }}}_{B}}{{h}_{\mathrm{core}}}{B}_{z}.\end{eqnarray} \tag{ 22 }$$

Since the crustal shear modulus is only ∼10⁻³ of the hydrostatic pressure, the displacement field (21) is nearly two-dimensional and incompressible. Therefore, both ${{\boldsymbol{\xi }}}_{B}$ and ${\boldsymbol{\xi }}$ can be expressed using stream functions Ψ_B, ${\rm{\Psi }}={{\rm{\Psi }}}_{\mathrm{el}}+{{\rm{\Psi }}}_{\mathrm{pl}}$, with, e.g., ${\xi }_{i}={\epsilon }_{{ij}}{\partial }_{j}{\rm{\Psi }}$. Here ${\epsilon }_{{ij}}=\pm 1$ is the antisymmetric symbol and i = x, y.

In this approximation, the wavelength of horizontal variations in Ψ_B and Ψ is constrained to be larger than or comparable to δR_c. We neglect any vertical structure in ${\boldsymbol{\xi }}$, that is, magnetic tilt interior to the crust, whose evolution will be considered in Section 6.

In some calculations, we also include a magnetic field localized to the crust. The seed horizontal magnetic field ${{\boldsymbol{B}}}_{\perp ,0}^{c}$ is uniform and configured independently of the core field. Elastic and plastic displacements of the crust perturb this field and also change the tilt of the core field,

$$\begin{eqnarray}{B}_{x}^{-} & = & \displaystyle \frac{{\partial }_{y}({\rm{\Psi }}-{{\rm{\Psi }}}_{B})}{{h}_{\mathrm{core}}};\quad {B}_{y}^{-}=-\displaystyle \frac{{\partial }_{x}({\rm{\Psi }}-{{\rm{\Psi }}}_{B})}{{h}_{\mathrm{core}}};\\ {B}_{x}^{c}-{B}_{x,0}^{c} & = & {B}_{x,0}^{c}{\partial }_{x}{\xi }_{x}+{B}_{y,0}^{c}{\partial }_{y}{\xi }_{x}\\ & = & {B}_{x,0}^{c}{\partial }_{x}{\partial }_{y}{\rm{\Psi }}+{B}_{y,0}^{c}{\partial }_{y}^{2}{\rm{\Psi }};\\ {B}_{y}^{c}-{B}_{y,0}^{c} & = & {B}_{x,0}^{c}{\partial }_{x}{\xi }_{y}+{B}_{y,0}^{c}{\partial }_{y}{\xi }_{y}\\ & = & -\,{B}_{x,0}^{c}{\partial }_{x}^{2}{\rm{\Psi }}-{B}_{y,0}^{c}{\partial }_{x}{\partial }_{y}{\rm{\Psi }}.\end{eqnarray} \tag{ 23 }$$

A tangled component of the magnetic field is easy to implement in Fourier space, as discussed below.

The equation of magnetoelastic equilibrium in the horizontal direction reads

$$\begin{eqnarray}&&{\partial }_{z}\left(\displaystyle \frac{{B}_{i}{B}_{z}}{4\pi }+\mu {\sigma }_{{iz}}\right)+\mu {\partial }_{j}({\sigma }_{{ij}}+{\sigma }_{{ij}}^{M})=0,\end{eqnarray} \tag{ 24 }$$

where the solid stress is

$$\begin{eqnarray}{\sigma }_{{ij}} & = & \mu ({\partial }_{i}{\xi }_{\mathrm{el},j}+{\partial }_{j}{\xi }_{\mathrm{el},i})+P{\delta }_{{ij}};\\ {\partial }_{j}{\sigma }_{{ij}} & = & \mu {\epsilon }_{{ik}}{\partial }_{k}({\partial }^{2}{{\rm{\Psi }}}_{\mathrm{el}})+{\partial }_{i}P,\end{eqnarray} \tag{ 25 }$$

the horizontal Maxwell stress is

$$\begin{eqnarray}&&{\sigma }_{{xx}}^{M}=\displaystyle \frac{{({B}_{x}^{c})}^{2}-{({B}_{y}^{c})}^{2}}{8\pi }=-{\sigma }_{{yy}}^{M};\quad {\sigma }_{{xy}}^{M}=\displaystyle \frac{{B}_{x}^{c}{B}_{y}^{c}}{4\pi },\end{eqnarray} \tag{ 26 }$$

and

$$\begin{eqnarray}{\partial }_{i}{\sigma }_{{xi}}^{M} & = & \displaystyle \frac{{B}_{y}^{c}}{4\pi }({\partial }_{y}{B}_{x}^{c}-{\partial }_{x}{B}_{y}^{c})\\ & \simeq & \,\displaystyle \frac{{B}_{y,0}^{c}}{4\pi }({B}_{x,0}^{c}{\partial }_{x}{\partial }^{2}{\rm{\Psi }}+{B}_{y,0}^{c}{\partial }_{y}{\partial }^{2}{\rm{\Psi }});\\ {\partial }_{i}{\sigma }_{{yi}}^{M} & = & -\,\displaystyle \frac{{B}_{x}^{c}}{4\pi }({\partial }_{y}{B}_{x}^{c}-{\partial }_{x}{B}_{y}^{c})\\ & \simeq & \,-\,\displaystyle \frac{{B}_{x,0}^{c}}{4\pi }({B}_{x,0}^{c}{\partial }_{x}{\partial }^{2}{\rm{\Psi }}+{B}_{y,0}^{c}{\partial }_{y}{\partial }^{2}{\rm{\Psi }}).\end{eqnarray} \tag{ 27 }$$

Only the elastic displacement field enters into the solid stress (25). We neglect forces arising from background gradients in ${{\boldsymbol{B}}}_{\perp ,0}^{c}$, in which case equality holds in the second line after Equation (27).

We integrate Equation (24) in vertical coordinate z from the bottom to the top of the crust, neglecting the exterior⁴ Maxwell stress ${B}_{i}^{+}{B}_{z}/4\pi$ and defining vertically averaged shear modulus and pressure via $\int \mu {dz}=\bar{\mu }\delta {R}_{c}$ and $\int {Pdz}=\bar{P}\delta {R}_{c}$. This gives

$$\begin{eqnarray}&&\displaystyle \frac{{B}_{z}^{2}}{4\pi }\displaystyle \frac{{\partial }_{y}({{\rm{\Psi }}}_{B}-{{\rm{\Psi }}}_{\mathrm{el}})}{{h}_{\mathrm{core}}}+\delta {R}_{c}\phantom{\rule{}{3.2ex}}[\bar{\mu }{\partial }_{y}{\partial }^{2}{{\rm{\Psi }}}_{\mathrm{el}}\\ &&\quad \left.+\displaystyle \frac{{B}_{y,0}^{c}}{4\pi }({B}_{x,0}^{c}{\partial }_{x}{\partial }^{2}{\rm{\Psi }}+{B}_{y,0}^{c}{\partial }_{y}{\partial }^{2}{\rm{\Psi }})+{\partial }_{x}\bar{P}\right]=0;\\ &&\displaystyle \frac{{B}_{z}^{2}}{4\pi }\displaystyle \frac{{\partial }_{x}({{\rm{\Psi }}}_{B}-{{\rm{\Psi }}}_{\mathrm{el}})}{{h}_{\mathrm{core}}}+\delta {R}_{c}\phantom{\rule{}{3.2ex}}[\bar{\mu }{\partial }_{x}{\partial }^{2}{{\rm{\Psi }}}_{\mathrm{el}}\\ &&\quad \left.+\displaystyle \frac{{B}_{x,0}^{c}}{4\pi }({B}_{x,0}^{c}{\partial }_{x}{\partial }^{2}{\rm{\Psi }}+{B}_{y,0}^{c}{\partial }_{y}{\partial }^{2}{\rm{\Psi }})-{\partial }_{y}\bar{P}\right]=0.\end{eqnarray} \tag{ 28 }$$

One observes that the pressure perturbation must be included only when ${{\boldsymbol{B}}}_{\perp ,0}^{c}$ is nonvanishing.

One solves for $\bar{P}$ by taking x- and y- derivatives of the first and second lines in Equation (28) and subtracting:

$$\begin{eqnarray}&&\bar{P}=\displaystyle \frac{{B}_{x,0}^{c}{B}_{y,0}^{c}}{4\pi }({\partial }_{y}^{2}-{\partial }_{x}^{2}){\rm{\Psi }}+\displaystyle \frac{{({B}_{x,0}^{c})}^{2}-{({B}_{y,0}^{c})}^{2}}{4\pi }{\partial }_{x}{\partial }_{y}{\rm{\Psi }}.\end{eqnarray} \tag{ 29 }$$

Substituting this back into Equation (28) gives

$$\begin{eqnarray}&&{\epsilon }_{{ij}}{\partial }_{j}{\rm{\Phi }}=0,\end{eqnarray} \tag{ 30 }$$

where

$$\begin{eqnarray}{\rm{\Phi }} & = & \displaystyle \frac{{B}_{z}^{2}}{4\pi }\displaystyle \frac{{{\rm{\Psi }}}_{B}-{\rm{\Psi }}}{{h}_{\mathrm{core}}}+\delta {R}_{c}\phantom{\rule{}{3ex}}[\bar{\mu }{\partial }^{2}{{\rm{\Psi }}}_{\mathrm{el}}\\ & & +\left.\displaystyle \frac{1}{4\pi }{({B}_{x,0}^{c}{\partial }_{x}+{B}_{y,0}^{c}{\partial }_{y})}^{2}({{\rm{\Psi }}}_{\mathrm{el}}+{{\rm{\Psi }}}_{\mathrm{pl}})\right]=0.\end{eqnarray} \tag{ 31 }$$

Uniqueness of the solution implies that Φ vanishes.

The elastic displacement stream function Ψ_el is easily solved for in Fourier space (k_x, k_y) in terms of Ψ_B and the cumulative plastic stream function Ψ_pl:

$$\begin{eqnarray}{{\rm{\Psi }}}_{\mathrm{el},{\boldsymbol{k}}} & = & -\,{\left[1+\delta {R}_{c}{h}_{\mathrm{core}}\left(\displaystyle \frac{4\pi \bar{\mu }}{{B}_{z}^{2}}{k}^{2}+\displaystyle \frac{{({\boldsymbol{k}}\cdot {{\boldsymbol{B}}}_{0}^{c})}^{2}}{{B}_{z}^{2}}\right)\right]}^{-1}\\ & & \times \,\left\{-{{\rm{\Psi }}}_{B,{\boldsymbol{k}}}+{{\rm{\Psi }}}_{\mathrm{pl},{\boldsymbol{k}}}\left[1+\delta {R}_{c}{h}_{\mathrm{core}}\displaystyle \frac{{({\boldsymbol{k}}\cdot {{\boldsymbol{B}}}_{0}^{c})}^{2}}{{B}_{z}^{2}}\right]\right\}.\end{eqnarray} \tag{ 32 }$$

To introduce a tangled magnetic field B_t with no preferred direction, we replace ${({\boldsymbol{k}}\cdot {{\boldsymbol{B}}}_{0}^{c})}^{2}\to {({\boldsymbol{k}}\cdot {{\boldsymbol{B}}}_{0}^{c})}^{2}+{k}^{2}{B}_{t}^{2}$ in Equation (32).

3.2.1. Numerical Procedure

Our procedure therefore is as follows:

1. The computational domain is a square in the x–y plane with periodic boundary conditions, and $4\pi {R}_{\star }^{2}/\delta {R}_{c}^{2}\sim {({2}^{6})}^{2}$ pixels. Runs with (2⁵)² or (2⁷)² pixels have been performed and show qualitatively similar results. The unit of length L is normalized to give a total area L² = 4π × (10 km)². The two-dimensional balance between boundary Maxwell stress and horizontal elastic stress is calculated using δR_c = 0.3 km, with density 1 × 10¹⁴ g cm⁻³ and shear modulus corresponding to the Negele & Vautherin (1973) equation of state and gravity g = 2 × 10¹⁴ cm s⁻². The specific heat is approximated in the same way as in the zero-dimensional model of Section 2, with the contribution from non-superfluid neutrons neglected. The strength and energy of the core magnetic field are calculated using an effective core depth⁵ h_core = L/(2π).

2. A Maxwell stress is imposed at the lower crust boundary with a power-law spectrum

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{B}_{x,y}{B}_{z}}{4\pi }\right)}_{k}\propto {({k}_{\parallel }^{2}+\eta {k}_{\perp }^{2})}^{-(1+\alpha )}\end{eqnarray} \tag{ 33 }$$

and random phases.⁶ This corresponds to a stress ${B}_{x}{B}_{z}/4\pi \propto {k}^{-2\alpha }$ per logarithm of wavenumber. The parameter η can be chosen to be less than unity to represent an anisotropic Maxwell stress, e.g., a nearly axially symmetric core field with predominantly latitudinal gradients. The basis (k_∥, k_⊥) can also be rotated with respect to (k_x, k_y): we typically choose a 45° rotation when η < 1, so as to eliminate the possibility of grid-alignment effects.

3. The magnitude A_z0 of the applied Maxwell stress is gradually raised until the peak creep rate is a small value (∼10⁻⁴ yr⁻¹). From this point on, the amplitude is raised linearly in time, dA_z/dt = A_z0/t_growth. Here t_growth is a free parameter, describing the speed of the transition to hydromagnetic instability in the core. We adaptively test the time step, raising or reducing it so as to maintain a maximum fractional change of ∼10⁻² in the creep rate in any pixel. As a result, the time step can vary from ∼10³ yr or longer down to a fraction of a millisecond. The creep rate, as determined by Equation (1), is capped at 0.1V_μ/δR_c, representing the feedback of inertial forces on unbalanced stresses.

4. In each time step δt, the plastic stream function is updated according to the method described in Section 3.3: ${{\rm{\Psi }}}_{\mathrm{pl},{\boldsymbol{k}}}\to {{\rm{\Psi }}}_{\mathrm{pl},{\boldsymbol{k}}}+{\dot{{\rm{\Psi }}}}_{\mathrm{pl},{\boldsymbol{k}}}\delta t$, where ${\dot{{\rm{\Psi }}}}_{\mathrm{pl},{\boldsymbol{k}}}$ is given by Equation (48). In the linear approximation, the update in the deformation field due to plastic flow can be separated from the back-reaction on the applied Maxwell stress and the self-consistent elastic response of the crust. The core magnetic stream function is updated according to ${{\rm{\Psi }}}_{B}\to {{\rm{\Psi }}}_{B}\mp {\dot{{\rm{\Psi }}}}_{\mathrm{pl},{\boldsymbol{k}}}\delta t$. Then the elastic response is recalculated according to Equation (32). The temperature is updated according to Equation (9).

5. The upper or lower sign in the update for Ψ_B corresponds to a positive or negative core spring constant. The change in core field due to creep supplements the imposed linear rise in the Maxwell stress and mediates a runaway global instability when the core spring constant is negative.

6. We separately test the presence or absence of a horizontal crustal magnetic field. In particular, a tangled crustal field has a strong effect on suppressing runaway yielding (Section 4).

Here we make use of the ab initio relation ${\sigma }_{b}({\dot{\varepsilon }}_{\mathrm{pl}})$ (Equation (1)) between breaking strain and strain rate to describe inhomogeneous plastic creep in one and two dimensions. The generalization to fully three-dimensional deformations is not addressed.

The molecular dynamics simulations of three-dimensional solids on which this is based take the imposed shear to have planar symmetry and grow at a uniform rate, e.g., $\varepsilon ={\dot{\varepsilon }}_{\mathrm{pl}}t$ with ${\dot{\varepsilon }}_{\mathrm{pl}}={\partial }_{x}{v}_{y}=\mathrm{const}$ (Chugunov & Horowitz 2010; Hoffman & Heyl 2012).

The plastic strain rate ${\dot{\varepsilon }}_{\mathrm{pl}}$ is defined as a positive scalar function of the scalar stress σ. In the planar model described here, it is given by

$$\begin{eqnarray}&&\sigma =\mu \varepsilon \equiv 2\mu {({\varepsilon }_{{xx}}^{2}+{\varepsilon }_{{xy}}^{2})}^{1/2},\end{eqnarray} \tag{ 34 }$$

where

$$\begin{eqnarray}&&{\varepsilon }_{{xy}}=\displaystyle \frac{1}{2}({\partial }_{y}^{2}-{\partial }_{x}^{2}){{\rm{\Psi }}}_{\mathrm{el}};\quad {\varepsilon }_{{xx}}=-{\varepsilon }_{{yy}}={\partial }_{x}{\partial }_{y}{{\rm{\Psi }}}_{\mathrm{el}}.\end{eqnarray} \tag{ 35 }$$

The first step is therefore to write the creep rate as a function of scalar strain, ${\dot{\varepsilon }}_{\mathrm{pl}}={\dot{\varepsilon }}_{\mathrm{pl}}(\varepsilon ),$ at each position in the solid.

The full evolution of the strain tensor by plastic creep is described by a single scalar function ${\dot{{\rm{\Psi }}}}_{\mathrm{pl}}$. Here we treat all parts of the crust as effectively plastic, but with exponentially varying creep rate. As long as the deformation rate is low enough that inertial forces can be neglected, the crust moves through a series of quasi-equilibrium states. Then the new elastic equilibrium can be obtained from Equation (32) after the plastic stream function is updated to ${{\rm{\Psi }}}_{\mathrm{pl}}\to {{\rm{\Psi }}}_{\mathrm{pl}}+{\dot{{\rm{\Psi }}}}_{\mathrm{pl}}\delta t.$

The procedure in outline is therefore

$$\begin{eqnarray}&&\varepsilon \ \mathop{\to }\limits^{{\boldsymbol{x}}}\ {\dot{\varepsilon }}_{\mathrm{pl}}\ \mathop{\to }\limits^{{\boldsymbol{k}}}\ {\dot{{\rm{\Psi }}}}_{\mathrm{pl}}\ \mathop{\to }\limits^{{\boldsymbol{k}}}\ {\dot{\varepsilon }}_{{ij},\mathrm{pl}}.\end{eqnarray} \tag{ 36 }$$

Here the label ${\boldsymbol{x}}$ or ${\boldsymbol{k}}$ refers to whether the step takes place in real space or Fourier space. As a final step, we can derive a scalar strain rate $\dot{\varepsilon }^{\prime}$ directly from ${\dot{\varepsilon }}_{{ij}}$ using Equation (34). Our approach to two-dimensional creep starts by equating $\dot{\varepsilon }^{\prime} =\dot{\varepsilon }$ and working backward.

3.3.1. Inhomogeneous One-dimensional Creep and Crack Formation

Some subtleties are involved in obtaining a self-consistent time evolution equation for Ψ_pl in a compact domain. Consider first the case of linear creep in cylindrical symmetry with one periodic coordinate,

$$\begin{eqnarray}&&{\boldsymbol{\xi }}({\boldsymbol{x}})={\xi }_{y}(x)\hat{y};\quad {\xi }_{y}(x+L)={\xi }_{y}(x).\end{eqnarray} \tag{ 37 }$$

At first it might appear straightforward to apply the continuous creep law (1) to the stress profile $\sigma =\mu | {\xi }_{y,\mathrm{el}}^{\prime }|$:

$$\begin{eqnarray}&&{\dot{\xi }}_{y}^{\prime }=\mathrm{sgn}({\xi }_{y,\mathrm{el}}^{\prime })\dot{\varepsilon }(\varepsilon ),\end{eqnarray} \tag{ 38 }$$

where $\dot{\varepsilon }={\dot{\varepsilon }}_{\mathrm{pl}}+{\dot{\varepsilon }}_{\mathrm{el}}$ and $\varepsilon \equiv | {\xi }_{y,\mathrm{el}}^{\prime }|$. But an obstruction arises from the fact that the net displacement

$$\begin{eqnarray}&&{\rm{\Delta }}{\xi }_{y}={\int }_{0}^{L}{\xi }_{y}^{\prime }{dx}=\mathrm{const}\end{eqnarray} \tag{ 39 }$$

must be invariant under continuous changes in ξ_y, whereas the integral

$$\begin{eqnarray}&&{\rm{\Delta }}{\dot{\xi }}_{y}={\int }_{0}^{L}\mathrm{sgn}({\xi }_{y,\mathrm{el}}^{\prime })\dot{\varepsilon }{dx}\ne 0\end{eqnarray} \tag{ 40 }$$

does not vanish in general.

Equivalently, Equation (38) does not self-consistently prescribe the time evolution of the zero-wavenumber mode of Ψ. To see this, we can work in Fourier space,

$$\begin{eqnarray}&&{\rm{\Psi }}=\displaystyle \sum _{r=-N/2}^{N/2-1}{{\rm{\Psi }}}_{r}{e}^{i2\pi {rx}/L}.\end{eqnarray} \tag{ 41 }$$

We first use the identity ${\varepsilon }_{{xy}}={\xi }_{y,\mathrm{el}}^{\prime }={\partial }_{x}^{2}{{\rm{\Psi }}}_{\mathrm{el}}$ and write the plastic term on the right-hand side of Equation (38) as ${\partial }_{x}^{2}{{\rm{\Psi }}}_{\mathrm{el}}\times {\dot{\varepsilon }}_{\mathrm{pl}}(\varepsilon )/{\varepsilon }_{}$ before Fourier transforming, giving⁷

$$\begin{eqnarray}&&-{r}^{2}{\dot{{\rm{\Psi }}}}_{\mathrm{pl},r}=-\displaystyle \frac{1}{N}\displaystyle \sum _{l=-N/2}^{N/2-1}{l}^{2}{{\rm{\Psi }}}_{\mathrm{el},l}{\left(\displaystyle \frac{{\dot{\varepsilon }}_{\mathrm{pl}}}{{\varepsilon }_{}}\right)}_{r-l}.\end{eqnarray} \tag{ 42 }$$

This equation is not satisfied for r = 0. This conundrum exists only if the deformation is continuous. Introducing a dislocation near the location of peak strain would remove Equation (39) as a conserved quantity and restore consistency with the local creep law (38).

3.3.2. Self-consistent Two-dimensional Creep

A different constraint on the creep law is present in a two-dimensional solid with inhomogeneous deformation. Here it is natural to consider aligning the plastic deformation tensor with the stress tensor (Hill 1998)

$$\begin{eqnarray}&&{\dot{\varepsilon }}_{{ij},\mathrm{pl}}={\varepsilon }_{{ij}}\displaystyle \frac{{\dot{\varepsilon }}_{\mathrm{pl}}}{\varepsilon }.\end{eqnarray} \tag{ 43 }$$

This solves the expression

$$\begin{eqnarray}&&{\dot{\varepsilon }}_{\mathrm{pl}}=2{\varepsilon }^{-1}({\varepsilon }_{{xx}}{\dot{\varepsilon }}_{{xx},\mathrm{pl}}+{\varepsilon }_{{xy}}{\dot{\varepsilon }}_{{xy},\mathrm{pl}}),\end{eqnarray} \tag{ 44 }$$

which is obtained by taking the time derivative of Equation (34). However, the strain tensor is derived from an incompressible displacement field and must therefore satisfy consistency rules. Equation (35) implies that

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}({\partial }_{y}^{2}-{\partial }_{x}^{2}){\dot{\varepsilon }}_{{xx},\mathrm{pl}}={\partial }_{x}{\partial }_{y}{\dot{\varepsilon }}_{{xy},\mathrm{pl}}.\end{eqnarray} \tag{ 45 }$$

The update $\delta {\varepsilon }_{{ij},\mathrm{pl}}={\dot{\varepsilon }}_{{ij},\mathrm{pl}}\delta t$ obtained from Equation (43) over a small time interval δt will satisfy Equation (45) when ${\dot{\varepsilon }}_{\mathrm{pl}}$ is constant, and also when the strain field has a planar or rotational symmetry, but not more generally.

Expressing the ansatz (43) in terms of the time derivative ${\dot{{\rm{\Psi }}}}_{\mathrm{pl}}$ gives two separate equations. Although these cannot be individually correct, because they do not satisfy the consistency rule, we will combine them by calculating the quantity ${\partial }_{x}{\partial }_{y}{\dot{\varepsilon }}_{{xx},\mathrm{pl}}+\tfrac{1}{2}({\partial }_{y}^{2}-{\partial }_{x}^{2}){\dot{\varepsilon }}_{{xy},\mathrm{pl}}=\tfrac{1}{4}{\partial }^{4}{\dot{{\rm{\Psi }}}}_{\mathrm{pl}}$. This gives

$$\begin{eqnarray}{\partial }^{4}{\dot{{\rm{\Psi }}}}_{\mathrm{pl}} & = & 2{\partial }_{x}{\partial }_{y}\left[2{\partial }_{x}{\partial }_{y}{{\rm{\Psi }}}_{\mathrm{el}}\displaystyle \frac{{\dot{\varepsilon }}_{\mathrm{pl}}}{\varepsilon }\right]\\ & & +\,({\partial }_{y}^{2}-{\partial }_{x}^{2})\left[({\partial }_{y}^{2}-{\partial }_{x}^{2}){{\rm{\Psi }}}_{\mathrm{el}}\displaystyle \frac{{\dot{\varepsilon }}_{\mathrm{pl}}}{\varepsilon }\right].\end{eqnarray} \tag{ 46 }$$

One can check that this ansatz reduces to Equation (42) when the strain field has a planar symmetry, e.g., when only ε_xy or ε_xx = −ε_yy, or some constant linear combination of these components, is nonvanishing throughout the solid.⁸

We then convert to Fourier space,

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\mathrm{pl}}(x,y)=\displaystyle \sum _{r,s=-N/2}^{N/2-1}{{\rm{\Psi }}}_{\mathrm{pl}r,s}{e}^{2\pi {ir}(x/L)}{e}^{2\pi {is}(x/L)},\end{eqnarray} \tag{ 47 }$$

giving

$$\begin{eqnarray}{\dot{{\rm{\Psi }}}}_{\mathrm{pl}r,s} & = & \displaystyle \frac{2{rs}}{{({r}^{2}+{s}^{2})}^{2}}\displaystyle \sum _{l,m}\,2{lm}{{\rm{\Psi }}}_{{lm},\mathrm{el}}{\left(\displaystyle \frac{{\dot{\varepsilon }}_{\mathrm{pl}}}{{\varepsilon }_{}}\right)}_{r-l,s-m}\\ & & +\,\displaystyle \frac{{r}^{2}-{s}^{2}}{{({r}^{2}+{s}^{2})}^{2}}\displaystyle \sum _{l,m}\,({l}^{2}-{m}^{2}){{\rm{\Psi }}}_{{lm},\mathrm{el}}{\left(\displaystyle \frac{{\dot{\varepsilon }}_{\mathrm{pl}}}{{\varepsilon }_{}}\right)}_{r-l,s-m}.\end{eqnarray} \tag{ 48 }$$

Similarly to the one-dimensional example discussed above, this equation does not evolve the (r = 0, s = 0) Fourier mode.

The two-dimensional description of the strain pattern must break down on length scales shorter than the crust scale height δR_c, and so we impose a cutoff N ∼ R_⋆/δR_c in Fourier space.

Dissipation in the crust and the associated change in internal magnetic energy are calculated as follows. The change in the plastic deformation stream function (48) is converted back to coordinate space and then used to calculate the time derivative of the strain tensor, ${\dot{\varepsilon }}_{\mathrm{pl}{xy}}=\tfrac{1}{2}({\partial }_{y}^{2}-{\partial }_{x}^{2}){\dot{{\rm{\Psi }}}}_{\mathrm{pl}}$, ${\dot{\varepsilon }}_{\mathrm{pl}{xx}}={\partial }_{x}{\partial }_{y}{\dot{{\rm{\Psi }}}}_{\mathrm{pl}}$. The dissipated energy over a time interval δt is $\int {dz}\,2({\sigma }_{{xy}}{\dot{\varepsilon }}_{\mathrm{pl}{xy}}+{\sigma }_{{xx}}{\dot{\varepsilon }}_{\mathrm{pl}{xx}})\delta t\,=$ $4\bar{\mu }({\varepsilon }_{{xx}}\cdot {\dot{\varepsilon }}_{\mathrm{pl}{xx}}+{\varepsilon }_{{xy}}\cdot {\dot{\varepsilon }}_{\mathrm{pl}{xy}})\delta t$ per unit area of crust. The change in core magnetic energy per unit area of crust is ${h}_{\mathrm{core}}\cdot {B}_{i}^{-}\delta {B}_{i}^{-}/4\pi$, where $\delta {{\boldsymbol{B}}}^{-}$ is obtained from Equation (23) by substituting $\delta {\rm{\Psi }}-\delta {{\rm{\Psi }}}_{B}={\dot{{\rm{\Psi }}}}_{\mathrm{pl}}\delta t\,+\delta {{\rm{\Psi }}}_{\mathrm{el}}-{\dot{{\rm{\Psi }}}}_{B}\delta t$. Here ${\dot{{\rm{\Psi }}}}_{B}$ describes the imposed secular growth in the Maxwell stress at the crust–core boundary, and δΨ_el is the self-consistent elastic response obtained from the sum of $\delta {{\rm{\Psi }}}_{\mathrm{pl}}={\dot{{\rm{\Psi }}}}_{\mathrm{pl}}\delta t$ and $\delta {{\rm{\Psi }}}_{B}={\dot{{\rm{\Psi }}}}_{B}\delta t$ via Equation (32).

Figure 13 shows the dissipation profiles of two outbursts (the first and 12th) from the same calculation as in Figure 3. The duration of the main pulse (∼0.1 s) is comparable to the time for unbalanced stresses to redistribute around the star at the limiting creep speed of ∼0.1 V_μ. There is a decaying power-law tail of emission, scaling as t⁻¹ and with total energy ∼10⁻² of the main pulse energy, which is powered by continuing plastic creep. This extended dissipation is associated with the growth of a fault (by an increase in length and in width). Figure 14 shows an example of fault growth.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The outburst light curve as shown does not include emission from a trapped fireball (Thompson & Duncan 1995) or from current dissipation in the magnetosphere (Thompson et al. 2002; Beloborodov & Thompson 2007). Excess decaying persistent X-ray emission with a similar power-law index (−0.7) and relative energy was observed following the 1998 giant flare of SGR 1900+14 (Woods et al. 2001). Similar power-law behavior is seen following shorter SGR bursts (e.g., Ibrahim et al. 2001; Lenters et al. 2003), but our two-dimensional simulations do not have the resolution to adequately address the time profiles of these lower-energy events.

The outburst rise time is ∼10⁻³ s, which is comparable to that observed in giant flares. There are also repeated injections of energy near the peak of the dissipation, similar to the actual burst light curves recorded by Terasawa et al. (2005) and Tanaka et al. (2007). Our calculation represents the accumulation and release of internal stresses. Although magnetospheric currents will slightly lag the internal motions, one can expect current-driven instabilities eventually to be triggered, with slightly faster rise times over a scale of ∼10 km (Thompson & Duncan 1995; Lyutikov 2003; Gill & Heyl 2010). The internal crustal dissipation would experience a faster rise if the limiting creep speed were taken to be a larger multiple of V_μ (than 0.1).

Figure 15 provides an explicit demonstration of how the energy release (both internal [black points] and external [red points]) is distributed in a nonlocal manner. (The black lines in Figure 3 only show the shifting position of peak creep.) In fact, the peak of energy release does not necessarily coincide with the initial trigger point.

Zoom In Zoom Out Reset image size

The total flare energy repeatably reaches ∼10⁴⁵ erg (including changes in external magnetic energy) but does not approach the 6 × 10⁴⁶ erg that was observed from the 2004 giant flare of SGR 1900+14 (assuming isotropic emission; Hurley et al. 2005; Palmer et al. 2005). Larger energies are possible if the magnetar core passes through a transition to hydromagnetic instability. Based on our model, we infer that this was the case for the 2004 flare.

Repeat bursts of intermediate energy appear in these simulations and are also observed following giant flares (e.g., the ∼4 s burst of energy ∼10⁴³ erg detected the day after the 1998 giant flare; Ibrahim et al. 2001). An analysis of the energy distribution of bursts is beyond the scope of this paper; a proper analysis would require introducing physical scales smaller than δR_c ∼ 0.3 km. Sizable repeat bursts appear to be too common, in comparison with the limited sample of two well-resolved giant flares. In this regard, it should be kept in mind that the moving locus of peak activity allows the burst to probe crustal zones with differing internal magnetizations and differing susceptibility to rapid creep.

We consider a vertical one-dimensional test problem that illustrates the dynamical ejection of a horizontal magnetic field. The crust is threaded by a uniform vertical magnetic field B_z = 10¹⁵ G that is bent by a variable amount in the horizontal x-direction, ${B}_{x}(z)={\partial }_{z}{\xi }_{B,x}{B}_{z}$. Our method assumes that the crust is in a state of slow plastic creep at the onset of the calculation: the strain rate ${\dot{\varepsilon }}_{\mathrm{pl}}$ is taken as a dynamical variable, and the solid stress is computed as a function of ${\dot{\varepsilon }}_{\mathrm{pl}}$. As a result, the calculation applies to parts of the crust, such as the shear bands described in Sections 2 and 4, in which the Maxwell stress B_xB_z/4π has been raised close to the solid yield stress.

We show that the timescale for the ejection of magnetic shear is at least 30 times shorter than the growth time of the global stress that triggers it. The zero-dimensional creep model of Section 2 suggests that the result will be qualitatively similar if the local Maxwell stress lies initially below the yield stress and runaway creep is preceded by a brief interval of growth in the total stress.

The creep rate is initialized to a uniform value at all depths z below the magnetar surface. (The results turn out to be similar when creep is relatively suppressed in the upper parts [lower densities] of the computational domain.) Therefore, the Maxwell stress B_xB_z/4π saturates the yield stress σ_b (Equation (3)). The corresponding transverse field, shown in Figure 16, equals B_z only at a particular depth. Everywhere else ${B}^{2}/4\pi =[({B}_{x}/{B}_{z})+({B}_{z}/{B}_{x})]{\sigma }_{b}\gt 2{\sigma }_{b}$, meaning that the field needed to induce plastic creep is generally stronger than (8πσ_b)^1/2 in the absence of an addition global stress.

Zoom In Zoom Out Reset image size

The vertical dependence of crustal density, melting temperature, and shear modulus are obtained by fitting the data from nuclear equation-of-state calculations, both below (Negele & Vautherin 1973) and above (Haensel et al. 1989) the neutron drip point, for a surface gravity g = 2 × 10¹⁴ cm s⁻². Smooth fitting functions are applied to the density and shear modulus profiles ρ(z), μ(z), which allows us to avoid complications associated with discrete jumps in these variables.

Figure 16 also shows the ratio between ${B}_{x}^{2}/(8\pi )$ and the initial elastic energy density $\tfrac{1}{2}{\sigma }_{b}^{2}/\mu$, as a function of depth in the crust. The magnetic field is the dominant energy reservoir near the base of the crust, where most of the potential energy is concentrated.

The evolution away from the initial slow creep state can be encapsulated in two variables

$$\begin{eqnarray}&&v\equiv {\partial }_{t}{\xi }_{B,x},\quad \chi \equiv {\partial }_{z}{\xi }_{B,x}.\end{eqnarray} \tag{ 54 }$$

Freezing of magnetic flux in the plastically deforming crust implies that

$$\begin{eqnarray}&&{\partial }_{t}\chi ={\partial }_{z}v,\end{eqnarray} \tag{ 55 }$$

so that v is the combined transverse velocity of magnetic field and crustal material. The transverse component of the relativistic Euler equation reads

$$\begin{eqnarray}&&{\partial }_{t}v={V}_{{\rm{A}}}^{2}{\partial }_{z}\chi +\displaystyle \frac{1}{{\rho }_{B}}{\partial }_{z}{\sigma }_{{xz}}.\end{eqnarray} \tag{ 56 }$$

The magnetic field dominates the effective inertia in the upper crust, ${\rho }_{B}=\rho +{B}_{z}^{2}/4\pi {c}^{2}$, so that the Alfvén speed ${V}_{{\rm{A}}}\,={B}_{z}/{(4\pi {\rho }_{B})}^{1/2}$ is limited to the speed of light c. Figure 17 shows the vertical dependence of V_A and transverse shear wave speed ${V}_{\mu }=\sqrt{\mu /{\rho }_{B}}$. The density field is extended to low values ($\rho \ll {B}_{z}^{2}/4\pi {c}^{2}$) where electromagnetic inertia dominates. The initial equilibrium state is

$$\begin{eqnarray}&&{\left(\chi +\displaystyle \frac{4\pi }{{B}_{z}^{2}}{\sigma }_{{xz}}\right)}_{t=0}={\chi }_{0},\end{eqnarray} \tag{ 57 }$$

where ${\chi }_{0}\sim { \mathcal O }(0.1)$ is the constant horizontal magnetic field (representing a shear or twist) that is imposed by the magnetosphere.

Zoom In Zoom Out Reset image size

Because we are studying the transition between magneto-elastic equilibrium and a dynamical state of the crust, we cannot assume that at each time step the crust relaxes to magneto-elastic balance. Therefore, we approximate ${\partial }_{z}v={\dot{\varepsilon }}_{\mathrm{pl}}$ and use Equation (3) to write

$$\begin{eqnarray}&&{\sigma }_{{xz}}={\sigma }_{b}(T,| {\partial }_{z}v| )\mathrm{sgn}({\partial }_{z}v).\end{eqnarray} \tag{ 58 }$$

The yield stress is insensitive to small changes in ${\dot{\varepsilon }}_{\mathrm{pl}}$ following Equation (3).

Near the top of the crust, μ is relatively small and V_A ∼ c ≫ V_μ. This means that translations of the upper crust follow passively a magnetic perturbation evolving below. This is especially clear in the case of horizontal shear deformations, e.g., a translation of the upper crust in the y-direction that is differential in x. These involve an approximate balance between the inertial force and the Lorentz force; during the first stages of the deformation, the plastic flow simply follows the magnetic field. In some circumstances the upper parts of the magnetar crust may actually be in a fluid state, but the melt depth of a quiescent magnetar is typically shallow enough that V_A ≃ c, and the solid-to-fluid transition has a small influence on the dynamics.

We impose an outgoing wave boundary condition at the upper boundary,

$$\begin{eqnarray}&&{({\partial }_{t}-{V}_{{\rm{A}}}{\partial }_{z})v| }_{z=0}=0.\end{eqnarray} \tag{ 59 }$$

This choice represents an external magnetic cavity large enough to sustain a full oscillation and is done in part to avoid dealing with complications associated with magnetospheric dissipation. In the deeper part of the crust V_A < V_μ. We keep the evolution short enough that the initial magnetoelastic equilibrium is only slightly disturbed near the lower boundary, and therefore we simply assume that the magnetic field is pinned there with constant Maxwell stress,

$$\begin{eqnarray}&&v(z=\delta {R}_{c},t)=0;\quad \chi (z=\delta {R}_{c},t)=\mathrm{const}.\end{eqnarray} \tag{ 60 }$$

To avoid numerical instabilities, the uniform creep rate ${\dot{\varepsilon }}_{\mathrm{pl}}$ is smoothly set to zero in a thin layer near the lower boundary. It is also essential to smoothly interpolate the sgn function in Equation (58), which otherwise is discontinuous when ∂_zv changes sign:

$$\begin{eqnarray}&&\mathrm{sgn}({\partial }_{z}v)\to \displaystyle \frac{{\partial }_{z}v}{\sqrt{{({\partial }_{z}v)}^{2}+{\dot{\varepsilon }}_{0}^{2}}}.\end{eqnarray} \tag{ 61 }$$

Here ${\dot{\varepsilon }}_{0}$ is a small parameter characterizing the thickness of the transition. With this modification, Equation (56) takes a form similar to a heat diffusion equation with a spatially variable thermal conduction coefficient.

Next, an external stress σ_xy is smoothly added to the system, representing, e.g., a horizontally propagating elastic wave. This means that, at a fixed vertical creep rate, the vertical stress is reduced to

$$\begin{eqnarray}&&{\sigma }_{{xz}}=\mathrm{sgn}({\partial }_{z}v)\sqrt{{[{\sigma }_{b}(T,{\partial }_{z}v)]}^{2}-{\sigma }_{{xy}}^{2}}.\end{eqnarray} \tag{ 62 }$$

Strictly speaking, ${\dot{\varepsilon }}_{\mathrm{pl}}^{2}={({\partial }_{z}v)}^{2}+{\dot{\varepsilon }}_{{xy}}^{2}$, where the second term is the horizontal strain rate. As the timescale for horizontal stress growth is longer by a factor of ∼R_⋆/δR_c, we approximate ${\dot{\varepsilon }}_{\mathrm{pl}}\approx | {\partial }_{z}v|$ during the one-dimensional evolution.

Collecting Equations (55)–(62), we obtain a complete system to model the growing stress imbalance and vertical hydromagnetic wave ejection.

Figure 18 shows the result of introducing an external "driver" wave stress

$$\begin{eqnarray}&&{\sigma }_{{xy}}={\varepsilon }_{{xy}}f(z)\sin (\omega t)\mu .\end{eqnarray} \tag{ 63 }$$

The radial eigenfunction f(z) reproduces the ${}_{2}{t}_{0}$ global elastic mode (McDermott et al. 1988), which has a frequency $\omega \simeq {2}^{1/2}{\bar{V}}_{\mu }/{R}_{\star }$, where ${\bar{V}}_{\mu }$ is the vertically averaged shear wave speed. To obtain numerically clean results, we take a driver amplitude ε_xy = 0.06. The background creep rate is $\dot{{\varepsilon }_{\mathrm{pl}}}$ = 1 s⁻¹, which is significantly slower than the ensuing evolution. The smoothing parameter ${\dot{\varepsilon }}_{0}$ in the stress formulae (61) and (62) is necessarily lower than this, and it is taken to be ∼0.1 s⁻¹, as limited by numerical precision.

Zoom In Zoom Out Reset image size

The net result is a burst of upward Poynting flux S_P that builds up significantly over a few milliseconds, tapping the horizontal magnetic field χ. The Poynting flux is proportional to χ,

$$\begin{eqnarray}&&{S}_{{\rm{P}}}=\displaystyle \frac{{B}_{z}^{2}}{4\pi }\chi v,\end{eqnarray} \tag{ 64 }$$

and is significantly enhanced in the upper crust if the shear has a baseline value χ₀ that is imposed by external magnetic twist or shear. The net electromagnetic energy ejected to the magnetosphere is

$$\begin{eqnarray}&&\int {dt}\,{S}_{{\rm{P}}}\approx \displaystyle \frac{{B}_{z}^{2}}{4\pi }\,{\chi }_{0}\,{\rm{\Delta }}x{| }_{z=0},\end{eqnarray} \tag{ 65 }$$

where Δx is the net horizontal displacement of the field lines. Taking χ₀ = 0.1, this works out to about 7% of the total energy loss over the first 10 ms of the evolution. The majority of the change in magnetic gradient energy is dissipated within the crust, at a rate ${\sigma }_{b}| \dot{{\varepsilon }_{\mathrm{pl}}}|$ per unit volume. Numerical tests show that this efficiency is insensitive to the amplitude of the driver elastic wave.

The magnetic tilt fluctuation is small near the stellar surface, χ ≃ v/V_A ∼ v/c, and consequently the plastic creep rate is also suppressed, $\dot{{\varepsilon }_{\mathrm{pl}}}\sim \chi /(10\,\mathrm{ms})\sim v/(c\,\cdot \,10\,\mathrm{ms})$. Combining this with the relatively small shear modulus, we find negligible plastic heating at shallow depths. The crust may behave like a fluid even while it is still in the solid phase.

What happens to the magnetic pulse after it is ejected from the crust depends on details such as the length of the magnetospheric flux bundle and the horizontal twist structure, which is not captured by this one-dimensional model. Effectively we have taken an infinitely long magnetospheric cavity, so that reflection of the outgoing "wave" can be neglected. This corresponds in a dipole geometry to a field line extending to a maximum radius ${R}_{\max }\sim \tfrac{1}{3}c\cdot 3\,\mathrm{ms}\sim 30\,{R}_{\star }$. Here we make the implicit assumption that twist ejection is accompanied by enough pair creation or particle ejection to sustain the MHD approximation in the magnetosphere. Then the twist propagates with a speed V_A ∼ c over a cavity of length ∼3R_max. In more compact parts of the magnetosphere, a quasi-static but strongly localized twist is established, which probably is subject to secondary current-driven instabilities. Details of the subsequent damping are not investigated here.

Stress buildup and release in the magnetar crust is a nonlocal phenomenon, in good part due to the small aspect ratio of crustal thickness to stellar radius. A yielding criterion involving a local balance between Maxwell stress and solid stress (Thompson & Duncan 1996; Perna & Pons 2011; Beloborodov & Levin 2014; Lander et al. 2015) does not offer an adequate description of runaway creep at dynamical rates. On the observational side, spectral analysis suggests that some slower transient behavior seen in the AXPs involves global response of the crust to local dissipation (Woods et al. 2004). The nonlocal yielding process described here also offers a framework for investigating small precursor events that have been detected shortly before giant magnetar flares (Hurley et al. 2005).

A consequence of the preceding item is that a standing elastic wave will modulate the creep rate at discrete plastic spots within the crust. The elastic wave becomes overstable if the covering fraction of the plastic spots is large enough, and if the crust is hydromagnetically unstable within these zones. The oscillation is not excited at a discrete time, thereafter passively decaying, but is continually reinforced.

Magnetar QPOs are therefore identified with periodically forced yielding in concentrated shear zones. This allows the emission zone of the modulated X-rays to be fairly compact, smaller than or comparable to R_⋆. Indeed, the strong rotational modulation detected in the QPO power (Strohmayer & Watts 2006) implies an occultation of the emission zone by the star and requires a compact size. A persistent oscillation of the magnetar crust also couples to extended dipolar magnetic field lines outside the star (Timokhin et al. 2008), but only relatively weakly (Thompson & Duncan 2001), and in a way that should be continuously visible over a full rotation.

This mechanism provides a promising explanation for the observed persistence of discrete QPOs in magnetar flare light curves (Israel et al. 2005; Strohmayer & Watts 2005; Watts & Strohmayer 2006). The objection of Levin (2006) to identifying magnetar QPOs with modes confined to the magnetar crust is weakened: we find that growth can be faster than the damping caused by coupling to a continuous spectrum of modes in the magnetized core. In other words, the relative concentration of the excited mode between the crust and the magnetized core (e.g., Levin 2007) may shift back toward the crust in the presence of strong forcing at a natural frequency of the crust. Active forcing also permits QPO activity to extend well beyond the core coupling time. We do find continuing plastic dissipation well after the impulsive energy release (Figure 13), but the hydromagnetic instability of the core is only treated schematically.

We have uncovered a plausible reason why super-Eddington outbursts are relatively rare events in magnetars (in comparison with slower but still dramatic changes in persistent X-ray output and spindown torque, both of which require strong magnetospheric currents). The extreme sensitivity of creep rate to applied stress and temperature can be compensated by the buildup of a reverse Maxwell stress during plastic flow. In a two-dimensional simulation, the seed magnetic field must exceed a substantial minimum strength. This need not suppress the formation of narrow zones of plastic creep, but it can help to cap the rate of creep within these zones. As the strength of this embedded field is increased, one sees a transition from a giant-flare-like phenomenon to a slower and less energetic outburst more typical of the transient magnetars. Super-Eddington outbursts may depend on a favorable magnetic geometry.

A transition from hydromagnetic stability to instability in the core is needed to explain the most energetic magnetar X-ray flares (>10⁴⁵ erg). A significant structural change in the core magnetic field must come with longer-term side effects; indeed, our simulations produce repeat bursts of lower energy. It should be taken seriously as the underlying driver of lower-energy SGR bursts as well. We show that such a transition can, even in the absence of a crustal response, produce transient effects over timescales of a day, significantly shorter than the thermal conduction time across the crust.

The ∼0.1 s duration seen in a wide range of magnetar X-ray bursts (e.g., Göǧüş et al. 2001) is a signature of stress redistribution within the crust. We demonstrate this for large (∼10⁴⁵ erg) energy releases using a global two-dimensional crust model that includes magnetic, elastic, plastic, and thermal effects. In major outbursts, the peak of plastic dissipation jumps by kilometers, often repeatedly. Models that postulate that most of the available energy is stored in a twisted magnetosphere at the onset of a burst (Lyutikov 2003) have not yet offered a quantitative explanation for the ∼0.1 s timescale.

Millisecond growth times in the dissipation rate arise from the nonlinear and temperature-dependent plastic response of the crust: they correspond to stress rebalancing over a distance ∼ δR_c ∼ 0.3 km. Rapid horizontal shear motions in the crust are accompanied by strong vertical magnetospheric currents and plausibly drive secondary current-driven instabilities (Thompson & Duncan 2001; Elenbaas et al. 2016). The shortest (submillisecond) growth seen in giant flares may point to such a magnetospheric phenomenon (Thompson & Duncan 1995; Lyutikov 2003).

Horizontal shuffling of poloidal magnetic field lines in the core is not fully mirrored by shuffling in the magnetosphere and therefore leaves behind vertical magnetic field gradients in the crust. We show that such embedded shear is rapidly (over several milliseconds) ejected from the crust when the nonlocal component of the stress builds up to a critical value. This is in line with early theoretical ideas suggesting the ejection of magnetic shear as the source of low-energy magnetar bursts (Thompson & Duncan 1995). The efficiency of this process is significantly enhanced if the background magnetic field is already sheared outside the neutron star. It will be further enhanced by secondary current-driven instabilities. The coupling of an elastic wave to the magnetosphere is slower by comparison (Blaes et al. 1989; Thompson & Duncan 2001; Link 2014).

This process starts with a decorrelation between solid and magnetic stresses throughout the crust, which is not present in a one-dimensional formulation. For example, restricting gradients to the vertical (z) direction and setting ${\partial }_{z}{T}_{{iz}}={\partial }_{z}({B}_{i}{B}_{z}/4\pi +{\sigma }_{{iz}})\,=0$, one finds a local relation between the maxima/minima of the solid stress and the Maxwell stress, ${\sigma }_{{iz}}=-{B}_{i}{B}_{z}/4\pi +\mathrm{const}.$ The constant boundary term does not change this result: it only biases up or down the local magnetic stress. By restricting the dynamics to one dimension, it is possible to consider a local extraction of magnetic energy through the deformation of the solid in which it is embedded. An example is the thermal-magnetic front studied by Beloborodov & Levin (2014) and Li et al. (2016), which operates on an intermediate timescale. Adding a growing global stress produces faster growth and larger energy release without the need for sharp temperature gradients.

Narrow fault-like structures easily form in the magnetar crust in response to unbalanced Maxwell stresses at the lower crust boundary. These structures have dissipative envelopes extending to scales much larger than δR_c and covering more than 10 percent of the surface area of the magnetar. As we show in Appendix A, they are well resolved, in the sense that the distribution of dissipation over the magnetar surface does not depend on resolution, and only a fraction of the total energy is dissipated in the highest-temperature core of the fault. We find that rapid development of extensive faults is associated with subsecond giant outbursts and is driven by the nonlocal response of the crust to rapid plastic flow in localized zones.

Fault-like structures provide a promising explanation for the spectral inference of hot spots on the surfaces of active magnetars (Ibrahim et al. 2001; Lenters et al. 2003; Woods et al. 2004; Esposito et al. 2007). They are especially prominent when the applied Maxwell stress has a rough rotational symmetry (one example being a kinked toroidal field in the outer core). The core magnetic field can otherwise be fairly smooth and need not contain current sheets.

In contrast with the analysis of Levin & Lyutikov (2012), we find that magnetic shearing within a fault need not suppress the release of energy. Instead, the resistance of the embedded magnetic field to the growing stress forces the shear layer to spread out as the global stress grows. In the presence of magnetar-strength (∼10¹⁵ G) fields, the fault width can approach the crustal thickness. This is energetically feasible because plastic creep with a narrow fault is driven by a stress accumulated over a much wider area of crust. Curvature of the fault in two dimensions is associated with a redistribution of stress outside the fault, creating wider dissipative structures.

Our two-dimensional model shows that the crust is nearly fully melted within active faults during major outbursts. We show that a standing Alfvén wave in a vertically magnetized liquid crust damps relatively slowly and oscillates with a roughly kilohertz frequency. Melted faults are promising locations for the high-frequency QPO activity seen in giant magnetar flares (Strohmayer & Watts 2006).

The same two-dimensional model shows extended plastic dissipation in the aftermath of a giant energy release, declining as a power law ∼t⁻¹ and releasing a net energy ∼10⁻² of the outburst energy. These properties are remarkably similar to the X-ray afterglow profile detected following the 1998 giant flare of SGR 1900+14 (Woods et al. 2001).

One extracts from this example a broader lesson that may apply to other instances of transient X-ray emission (e.g., to the slower outbursts of transient magnetars; Ibrahim et al. 2004; Mori et al. 2013; Alford & Halpern 2016, and references therein). The magnetospheric twist associated with transient emission need not be injected all at once, as in the model of Beloborodov (2009), and is probably localized closer to the star than is implied by the dipolar "j-bundle" construction.

Implications for the magnetar eigenvalue problem. The potential of magnetar QPOs as a diagnostic of the internal structure has inspired detailed approaches to the eigenvalue problem by several groups (Piro 2005; Glampedakis et al. 2006; Levin 2007; Samuelsson & Andersson 2007; Sotani et al. 2008; Cerdá-Durán et al. 2009; Colaiuda & Kokkotas 2011; Gabler et al. 2011; van Hoven & Levin 2011, 2012; Passamonti & Pons 2016). The eigensolutions are not only sensitive to magnetic and compositional effects in the elastic parts of the crust but also subtly modified by plastic flow. In particular, creep in localized faults may contribute to the observed QPO frequency drifts by introducing phase shifts into shear wave propagation.

Hard and persistent X-ray emission. Magnetars are copious sources of hard X-rays in quiescence, with a power greatly exceeding the spindown power (Murakami et al. 1994; Kouveliotou et al. 1999; Kuiper et al. 2006). This emission is probably a signature of very strong electric currents flowing through the closed magnetosphere. The source location has variously been ascribed to compact and mildly relativistic plasma near the magnetar surface (Thompson & Beloborodov 2005), or to outward relativistic flows of pairs along a thin bundle of closed field lines surrounding the magnetic dipole axis (Beloborodov 2013). Our two-dimensional simulations naturally produce compact currents associated with strong gradients in creep rate and distributed broadly across the crust. They therefore favor the first, more localized type of X-ray emission process. The pulsed radio emission detected from some transient magnetars (Camilo et al. 2006, 2007) also gives evidence for a dynamic magnetosphere, but it transmits much less energy than the X-rays and is much more rapidly variable, and it is consistent with emission from open (or nearly open) magnetic field lines (Thompson 2008).

Thermalization of magnetic shear energy ejected from magnetar crusts. We do not address how a time-dependent magnetic displacement will damp after injection into the magnetosphere. Emission of quasi-thermal X-rays during super-Eddington magnetar bursts depends on a transfer of energy from relatively large scale (≳δR_c) magnetic gradients to sub-relativistic particles with gyrational radii some 10⁻¹⁷ times smaller. Some type of cascade process must be involved (Thompson & Blaes 1998).

Plastic dissipation in the upper crust is found to be negligible, in agreement with Li & Beloborodov (2015). We nonetheless note that any cascade process that operates outside the star should also operate within the upper crust (where B²/8π ≫ ρc² and V_A ∼ c). The profile of heat deposition may therefore be flat compared with the pressure profile (Lyubarsky et al. 2002). Prompt burst afterglow emission is a probe of the thermalization process.

Surface Maxwell stresses. Persistent electric currents flowing through the surfaces of magnetars would cut the vertical imbalance in B_⊥B_z/4π. The outer core and inner crust generally can sustain a stronger magnetic twist than the magnetosphere, but this effect could be important in zones with modest internal twist. A reduction in magnetospheric twist by current-driven instabilities would increase the net Maxwell stress applied to the crust, thereby facilitating the growth of rapid creep. The magnetosphere could also have a nonlocal effect on crustal yielding by communicating twist along closed magnetic field lines: injected twist would be focused at the opposing magnetic footpoint, although at the cost of diluting the twist over the length of the flux bundle.

Do lower-energy SGR bursts represent a "bottom-up" or "top-down" phenomenon? Earthquakes are ultimately the consequence of large-scale convective motions and compositional divisions within Earth (Turcotte & Schubert 2002). We have shown that fault-like structures emerge naturally in magnetar crusts, but our global model does not resolve features smaller than ∼0.3 km. Therefore, we cannot offer a quantitative model for the origin of low-energy SGR bursts.

The ∼0.1 s characteristic duration of low-energy SGR bursts does highlight an important difference with earthquakes: this duration is comparable to the shear wave propagation time across the magnetar crust, pointing to strongly nonlocal energy release, whereas the terrestrial events are relatively short by the analogous measure. In spite of this distinction, the energy distributions of SGR X-ray bursts and earthquakes follow similar power laws, with smaller bursts supplying a smaller cumulative energy release than the largest flares (Cheng et al. 1996; Göǧüş et al. 2001).

Finally, we note that a role for a large-scale hydromagnetic instability in facilitating low-energy burst activity is suggested by the relative activity of the same sources that (so far) have produced giant flares.

What role does Hall drift play in driving transient magnetar phenomena? Our approach here is partly motivated by the fact that the rate of crustal Hall drift may be exceeded in young and hot magnetars by the rate by magnetic drift in the core (Haensel et al. 1990; Goldreich & Reisenegger 1992), given the strong feedback of magnetic dissipation on the temperature and therefore on the ability of the magnetic field to overcome the back-pressure induced by the radial electron-fraction gradient (Thompson & Duncan 1996).

The sign of the effect of Hall drift on magnetar activity is not clear. The time-dependent yielding calculations presented here raise the possibility that Hall drift has a calming effect, opposite to the one usually considered. Sourcing small-scale magnetic irregularities (Goldreich & Reisenegger 1992; Hollerbach & Rüdiger 2002; Wareing & Hollerbach 2009; Pons & Geppert 2010; Takahashi et al. 2011; Li et al. 2016) helps to create components of the magnetic field transverse to the local direction of shear, which we find are key to suppressing runaway creep. Alternatively, small zones of enhanced magnetic flux density driven by Hall drift could become preferred dissipation sites in the presence of time-dependent global stresses. The Lorentz force density increases with wavenumber in an electron-MHD cascade (Cho & Lazarian 2004), and localized zones of magnetic flux enhancement are seen in three-dimensional simulations of Hall transport (e.g., Gourgouliatos et al. 2016).

The relevance of Hall drift for large-scale reorganizations of the magnetic field depends on the location of the currents that support the stellar magnetic field. Large-scale simulations (Perna & Pons 2011; Gourgouliatos et al. 2016) achieve a large redistribution of the dipole magnetic flux within the active lifetimes of a magnetar (∼10³–10⁴ yr) only if the external field is anchored in the crust, so that the crustal toroidal field approaches 10¹⁶ G. Otherwise, the latitudinal electron drift is too slow to transport the poloidal magnetic field a significant distance. In the alternative approach advanced here, global twists, which can strongly modulate the spindown of the magnetar, are sourced by global imbalances in the core.