The game of life on a magnetar crust: from $\gamma$-ray flares to FRBs

This paper presents a model to unify the diverse range of magnetar activity, through the building and release of elastic stress from the crust. A cellular automaton drives both local and global yielding of the crust, leading to braiding of coronal loops and energy release. The model behaves like a real magnetar in many ways: giant flares and small bursts both occur, as well as periods of quiescence whose typical duration is either $\lesssim 1$ yr or $\sim 10-30$ yr. The burst energy distribution broadly follows an earthquake-like power law over the energy range $10^{40}-10^{45}\,{\rm erg}$. The local nature of coronal loops allows for the possibility of high-energy and fast radio bursts from the same magnetar. Within this paradigm, magnetar observations can be used to constrain the poorly-understood mechanical properties of the neutron-star crust.


INTRODUCTION
Magnetars, a class of restless neutron star characterised by energetic outbursts, drive a wide range of astrophysical phenomena: from short X-ray bursts of 10 41 erg, through storms of bursts and prolonged intermediate events, up to rare giant γ-ray flares of 10 44 − 10 46 erg, among the most violent events in the Universe (Turolla, Zane, & Watts 2015; Kaspi & Beloborodov 2017). Very recently, observations (Mereghetti et al. 2020;Bochenek et al. 2020;CHIME/FRB Collaboration et al. 2020;Li et al. 2021;Ridnaia et al. 2021) have shown us that magnetars are also the central engines for at least some fast radio bursts (FRBs); that some short γ-ray bursts are the result of giant flares from extragalactic magnetars (Burns et al. 2021); and that certain ultralong-period radio emitters may be old magnetars (Hurley-Walker et al. 2022;Caleb et al. 2022;Beniamini et al. 2022). In all its activity, the magnetar's solid crust plays a key role: it stores an enormous amount of energy in the form of elastic stress τ built up as the intense internal magnetic field B evolves; seismic events then release some of this energy into the corona, ultimately leading to the activity we observe (Ruderman 1991;Thompson & Duncan 1995;Perna & Pons 2011;Lander et al. 2015;Dehman et al. 2020).
Compared with the detailed quantitative simulations of crustal magnetic-field evolution (Pons & Viganò 2019;Gourgouliatos, De Grandis, & Igoshev 2022), our modelling of how the crust releases elastic energy is rudimentary, and faces conceptual challenges (Thompson, Yang, & Ortiz 2017). Furthermore, even the qualitative picture of magnetar activity is disjointed, with e.g. short bursts and giant flares generally treated as being of different physical origin, impeding any attempts to probe the underlying crustal physics. Here, by contrast, we show how the full spectrum of magnetar activity can be interpreted by taking a new approach: a single, physically-motivated model of the crust as a cellular automaton that drives coronal activity and thus the observed bursting behaviour. The outer crust of a magnetar, with a density ρ < 4 × 10 11 g cm −3 , is relatively weak and will be partially molten for younger stars; although it may be the source of weak bursts  we neglect its effect here, and concentrate on the more universal role of the inner crust, an immensely strong crystalline structure that resists and responds elastically to any imposed force, until it reaches its elastic yield stress τ el . When the crust eventually yields, the high pressure inhibits the formation and propagation of voids through the crustal lattice, so instead of a brittle fracture the crust is expected to flow plastically (Jones 2003), releasing elastic energy 1 .
To study crustal failure quantitatively, we will first need profiles of the mass and charge density ρ, ρ e and composition throughout the crust. These are found by solving the TOV stellar structure equations together with the SLy4 equation of state, as in Lander & Gourgouliatos (2019). Fixing the mass at 1.4M ⊙ gives us a model of radius 11.7 km, whose inner crust is 550 m thick, which we adopt throughout this paper.
To estimate the free energy reservoir of the crust, we evaluate the formula for τ el from Chugunov & Horowitz (2010) and volume-integrate it over the inner crust, to find a maximum elastic energy E el ∼ τ el dV = 4π τ el (r) r 2 dr ≈ 10 47 erg. The energy of the magnetic field threading the crust is comparable with E el ; if this is also tapped during crustal failure, the total energy reservoir becomes ∼ 2 × 10 47 erg. This is a factor of ∼ 10 greater than the most powerful known giant flare (Palmer et al. 2005); crustal energy alone is therefore able, in principle, to explain all magnetar activity observed to date.
Crustal yielding is an essential part of the magnetar paradigm, as it drives the transfer of energy, via the motion of embedded coronal field footpoints, out to the corona, from where it is released in the activity we observe (Lyutikov 2006). A major challenge in neutron-star physics is how yielding occurs; microscopic molecular-dynamics simulations exhibit collective local failure (Horowitz & Kadau 2009), but if every small group of crustal ions were to yield as soon as τ = τ el , all resultant bursts would be undetectably small. Equally, the crust's stress distribution is likely to be highly anisotropic, with regions where τ ≈ τ el and others φ θ v pl Figure 1. The magnetar crust as an array of cells. The left-hand bunch of coronal field lines are untwisted, embedded in a static cell in its elastic phase. The right-hand bunch are braided by a plastic flow v pl circulating around a cell whose elastic yield stress has been exceeded.
with τ ≪ τ el , so it is energetically disfavoured for every local failure to grow into a global one. Furthermore this would lead to a scenario, contradicted by observations, where a magnetar would be unable to produce repeated small bursts. An additional piece of physics must, therefore, set the characteristic lengthscale for crustal failure.
2.2. The crust as an array of cells First-principles macroscopic simulations of crustal failure are currently out of reach, but there are many clues to guide us to the origin of magnetar activity. Observationally, numerous studies have shown how bursting activity appears to come from relatively small patches of the crust (Palmer 1999;Younes et al. 2022), often at locations across the stellar surface (Scholz & Kaspi 2011;Younes et al. 2020), meaning that crustal failure must often be a local phenomenon. Both high-energy bursts (Cheng et al. 1996) and FRBs (Wadiasingh & Timokhin 2019) exhibit a power-law distribution of number vs energy, like that of earthquakes (Bak et al. 2002). From the theory side, 3D numerical evolutions of the crustal magnetic field show the development of ∼ 1 km 2 patches with strong B φ (Gourgouliatos, Wood, & Hollerbach 2016; Igoshev et al. 2021); this significant change from the initial B induces high stress (Lander & Gourgouliatos 2019), and so elastic failure is likely to occur within such patches, but not necessarily spread beyond them. On the other hand, a crustpowered giant flare requires a larger-scale failure, to explain the amount of energy released and how much quieter the star becomes afterwards; the model must therefore also allow for this possibility.
Guided by these considerations, we split the inner crust into an array of semi-autonomous cells with fixed boundaries, each of surface area 1 km 2 . To fix the cell depth, we note that since τ el increases by a factor of ∼ 1000 from the top to the base of the inner crust, we do not expect every local failure of a cell to propagate to the full 550 m depth of the inner crust. Instead we fix the cell depth at 200 m, over which τ el varies by less than an order of magnitude, and assume that the crystalline structure in such a cell always fails collectively. As a result, we can ignore variations of the stress within a cell and assume it is given by a single spatially constant τ that evolves with time. This evolution is dependent on other parameters (see Eq. 3), so for consistency we therefore assume there is no spatial variation in any physical quantity within a cell; we use the value of each from the base of the cell. Physical quan-tities do, however, vary from cell to cell across the crust.
When a cell's stress exceeds the elastic yield value, τ > τ el , only a fraction of it is expected to be relieved in the ensuing plastic flow (Lander & Gourgouliatos 2019); if we assume for definiteness a 10% reduction from τ el , the corresponding energy release is 3 × 10 40 erg -similar to a fairly powerful short X-ray burst, and providing a sanity check of the cell model.
We regard a cell's stress as being sourced by the magnetic field B alone, and define a scalar stress τ ≡ B 2 /4π. Neglecting the effect of Ohmic decay -reasonable for young magnetars (Pons & Viganò 2019) -the evolution of B in a crust stressed beyond τ el is dictated by an interplay of Hall drift and advection due to plastic flow v pl : (1) We model the crust as a Bingham plastic: below τ el its response to stress is purely elastic, with v pl = 0 and only the Hall drift term present in (1), whilst above τ el the crust behaves as a viscoplastic with flow velocity v pl ∝ (τ − τ el ). Specifically, we use a scalar version of the slow viscous-flow model of Lander & Gourgouliatos (2019), produced by replacing spatial derivatives where ν is the viscosity of crustal matter in its plastic phase. Now, from (1) we can derive an approximate scalar equation for τ , by using the identity ∂B 2 /∂t = 2B · ∂B/∂t, and eliminating B, v pl using (2) and the relation τ = B 2 /4π: where we have swapped the signs on the two right-hand-side terms to reflect the tendency of Hall drift to increase τ and plastic flow to reduce it. Simulations (Lander & Gourgouliatos 2019;Gourgouliatos & Lander 2021) show that stresses substantially higher than τ el may form before plastic flow has a chance to relieve them; to mimic this we model a cell's response as elastic, and allow τ to grow under the Hall effect, until it reaches a critical value τ = 1.1τ el . At this point failure occurs and in principle there will be an interplay between the Hall and plastic terms. To understand this, let us use (3) to estimate characteristic timescales t Hall , t pl for the evolution of stress under Hall drift and plastic flow: using typical values: L = 200 m, τ = 1.1τ el , ρ e = 1.4 × 10 26 esu cm −3 , ν = 10 36 poise. The Hall effect is always active, but may reasonably be neglected during the plastic phase, since t pl ≪ t Hall . During the elastic phase, only the Hall effect operates. Therefore, we model a single cell's evolution as alternating phases of growth of τ until the value 1.1τ el is reached, followed by a reduction of τ down to τ el under the action of plastic flow 2 .
Time-integrating (3), we find closed-form expressions τ = τ (t, t swap ) for both the Hall and plastic phases, where t swap is the time at which the current phase began. The Hall phase takes 75 yr, and the plastic phase ∼ 1 − 5 yr -somewhat shorter than the timescale estimates of (4), since here τ is only reduced by 10%.
The crust's temperature T has a minimum 'ambient' value T amb defined at the start of each simulation. We wish to consider a range for T amb that encompasses all neutron stars that might plausibly display bursting activity, from young magnetars whose high surface T appear to require a heat source in the crust , to old sources (> 10 4 yr) that have experienced no heating since birth. These considerations lead us to adopt, based on the cooling evolutions of Ho, Glampedakis, & Andersson (2012), the range 3 T amb,9 ≡ T amb /(10 9 K) = 0.05 − 0.5. Next, v pl causes heating at rate Q pl , which we model using an approximation to equation (21) from Li, Levin, & Beloborodov (2016): integrating from T = T amb at the start of the plastic phase, and where C V is the specific heat capacity. 60% of this heat is assumed to stay in the cell and the rest to diffuse to its surroundings, to mimic the effect of thermal conductivity. Once the plastic phase ends, T → T amb linearly over a timescale ≈ C V T /Q ν ≈ 13 yr (setting T 9 = 1 for the plasticallyheated cell, and using values for C V and neutrino emissivity Q ν from Gnedin, Yakovlev, & Potekhin (2001), evaluated at the base of the cell as usual). Finally, we need an expression for ν. Since this is unknown from first principles, numerical experimentation has been required to understand the range of plausible values that produce an interplay between Hall drift and plastic flow, and therefore allow for magnetar activity (Lander 2016).
Following Lander & Gourgouliatos (2019), the density-dependence of ν is taken to be the same as for τ el . The T -dependence of viscous fluids is often approximated by the Andrade equation (Andrade 1934), where ν ∝ e 1/T . Modifying this to avoid the divergent behaviour (problematic for a solid) as T → 0, and adjusting constants to match previous work (Lander & Gourgouliatos 2019;Gourgouliatos & Lander 2021;Kojima, Kisaka, & Fujisawa 2021), we arrive at the phenomenological relation ν(ρ, T ) = 5 × 10 5 τ el (ρ)e 5/(1+T9) poise, giving a possible range 10 34 ν 10 36 poise for a cell. Heating reduces ν, affecting the crust's activity in two opposing ways: on the one hand it increases v pl and the rate of coronal twisting; on the other hand it shortens the plastic phase and reduces the chances of several plastic cells having time to join up into a cluster, which in turn makes future deep failures and giant flares less likely. Any other relation for ν of the same order of magnitude, and reducing with T , would lead to broadly similar results. those with whom it shares an edge). It is well known that complex physics can arise from simple cellular automata (von Neumann 1968;Berlekamp, Conway, & Guy 2004). At the same time, these rules need to be linked to underlying physics as rigorously as possible to have any predictive power, and results cannot be artefacts of a fine-tuned model, but should be robust and generic: self-organised criticality (SOC) (Katz 1986;Bak, Tang, & Wiesenfeld 1988). There is evidence that both X-ray bursts and FRBs are driven by the same SOC process (Wei et al. 2021), motivating the present study.
Crustal magnetic field lines thread multiple cells, and are dragged around locally by v pl in a cell. This exerts a shearing force on its neighbours, but cannot cause them all to failotherwise every localised v pl could quickly propagate across the entire crust. Instead, we encode this effect through a cell rule: a cell in its elastic phase normally switches to its plastic phase at τ = 1.1τ el , but for every plastic neighbour, the cell's yield stress is lowered by 0.025τ el ; nearby plastic flow thus hastens, rather than triggers, a cell's failure. This is the key rule that leads to SOC-like behaviour of the model. A contiguous cluster of plastic cells is regarded as a single physical entity, with v pl circulating across the entire cluster with some average velocityv pl .
Shallow failures, down to the base of the cells, do not release enough energy to explain larger magnetar events, so the deeper crust -which will also be close to its τ el -must also fail sometimes. Our criterion for this to occur is that the cell itself, and at least three neighbours, must all be in a plastic phase simultaneously. Such a 'deep' failure below a cell releases all stored elastic energy in that region, down to the crust-core boundary. Thereafter we assume τ = 0 for that cell (replenishing τ back to τ el through the entire inner crust would take substantially longer than the previous estimate (4) of t Hall ≈ 800 yr); in this way, magnetar crustal dynamics are analogous to forest-fire models (Drossel & Schwabl 1992).

Corona
The most readily induced kind of plastic flow satisfies ∇ · v pl = 0 and has no radial component (Lander 2016).
Restricting v pl to a cell, the only permissible motion is a θ − φ circulation of matter in loops around the cell; the global crustal motion is inherently non-axisymmetric. The footpoints of external magnetic field lines are embedded in the cell, so v pl causes a braiding of these; see Fig 1.v pl within a contiguous cluster of plastic cells plays two roles: it increases the average twist ψ of the associated coronal loop, dψ/dt =v pl , and is also taken to represent the rate of transfer of elastic to coronal energy E clus for the cluster, so that: where E shallow , E deep are the sums of energy releases from all shallow and deep failures, respectively. The total coronal energy E corona is then the sum of all E clus . Magnetar bursts are linked to magnetic reconnection in the corona (Lyutikov 2006), a complex process that is not well understood; any attempt to implement a detailed prescription risks introducing several new poorly-constrained parameters, making the model harder to constrain or falsify afterwards. Instead we simply assume that when a plastic cluster ceases to exist, its remaining associated E clus is emitted as a single burst. In the special case where at least one cluster's (average) twist reaches a peak value ψ > 0.3 rad, we impose a 'hightwist' prescription where all coronal braids reconnect, ψ is reset to zero for each, and the total E corona is ejected at once in a 'giant flare'.

NUMERICAL CODE
To avoid conceptual issues where the two footpoints of a coronal field line might both move in such a way that v pl does not cause any increase in ψ, we assume one footpoint of every coronal field line is in the 'active' Northern hemisphere, and the partner footpoint in a 'passive' Southern hemisphere, releasing elastic energy but not driving the motion. Covering the Northern hemisphere at a radius 11.3 km (the outer boundary of the inner crust) requires ≈ 800 cells of 1 km 2 surface area; we also want a grid which is four times as long in φ as in θ (since 0 ≤ θ < π/2, 0 ≤ φ < 2π). We therefore choose a fiducial resolution of 14×56, i.e. 784 cells. For the top row of cells, identified with the pole, we set τ = 0; those around the equator are assumed to be mirrored (for the purposes of counting numbers of plastic neighbours) with an unmodelled set of partner cells in the Southern hemisphere. Periodic boundary conditions are imposed to identify the φ = 0 and φ = 2π edges of the grid. The network of cells is evolved for 1000 yr with a C++ code which tracks the formation, evolution and extinction of plastic clusters, with a default timestep of 0.01 yr. At the start of each simulation all cells are in the Hall phase, with stresses in the range 0.9 < τ /τ el < 1.1 randomly assigned to each cell. Thus, in this paper the time t = 0 yr represents a mature magnetar's crust at age ∼ 1000 yr; the newborn crust is unstressed, so no seismic activity will occur at this stage. Differences in magnetic-field strength and topology mean that some magnetars will reach this highly-stressed state earlier than others, but are otherwise not likely to result in radically different behaviour of the model.

RESULTS
Fig 2 shows a representative example light curve from 1000 yr of evolution. The first high-twist 'giant flare' event is seen at t = 85 yr. Fig 3 shows the crust's stress pattern before and after this event and demonstrates how every giant flare leaves behind extended patches of unstressed crust, thus reducing the chances of any future large-scale event occurring. The fractal pattern remaining after 1000 yr is characteristic of these simulations, and is seen at higher resolutions too.
E corona can be high for long periods without any individual loop developing high twist (although note that our simulations have no term for twist decay); in such a state the magnetar could be relatively quiet, but perhaps with a significantly non-dipolar spindown rate (making estimates of the external field unreliable) (Harding, Contopoulos, & Kazanas 1999;Thompson, Lyutikov, & Kulkarni 2002). Long periods of 'quiescence' (defined here as E corona < 10 42 erg) punctuated by occasional intermediate events are seen, especially at later times. Fig 4 shows durations of these quiescent periods, for two representative runs with different T amb . The distribution is roughly bimodal: our model magnetars go quiet for either a few months, or ∼ 10 − 30 yr. There is no correlation between the size of smaller events and the waiting time until the next one; though in the aftermath of a giant flare the model stars are -like real magnetars -often quieter.
The focus of this paper is on understanding the crustal dynamics and transfer of energy to the corona, a widely accepted idea for high-energy magnetar bursts that is now also a leading model to explain (at least some) FRBs. Radio emission, however, is likely to require substantially lower levels of coronal twist than high-energy emission (Wadiasingh & Timokhin 2019), making it hard to reconcile the standard globally-twisted magnetar corona model (e.g. Thompson, Lyutikov, & Kulkarni (2002)) with recent observations of contemporaneous FRBs and X-ray bursts from the magnetar SGR 1935+2154 (Mereghetti et al. 2020;Bochenek et al. 2020;CHIME/FRB Collaboration et al. 2020). By contrast, this paper's model, where crustal failure leads to a network of more localised coronal braids with varying levels of twist, naturally allows for this.
Whilst we cannot determine whether a given event from our evolutions will ultimately be seen as an X-ray burst or an FRB, since we do not study emission physics, we can still infer whether conditions in the corona are propitious for generation of a particular kind of radiation. From Fig. 2 we see that for the first 300 yr the corona always contains a lot of twist energy, making FRB emission highly disfavoured. Confirming observational results (Lin et al. 2020), we thus expect a classically active young magnetar to produce far more X-ray bursts than FRBs. Later on in our evolutions, energetic bursts of shorter duration occur, in between increasingly long periods of quiescence, and we thus anticipate -following Wadiasingh & Timokhin (2019) -that FRB emission will become more likely. ber of events does however increase with resolution, because there are more cells available to undergo elastic failure. Note that the evolution with 44 × 176 cells, ten times the fiducial value, is shown as an extreme case; the other three resolutions are likely to be more realistic (recall Section 2.2). The infrequency of E burst < 10 40 erg events is an artefact of our model, which considers the inner crust only; lower-energy bursts are likely to involve the outer crust. Our model predicts that no magnetar will produce more than ∼ 10 giant flares over its first 1000 yr of maturity and -given the results of Figs 2 and 3 -none at all thereafter. Cellular automaton models, and many different kinds of astrophysical source, generically exhibit power-law energy distributions with Γ ≈ 1.5 − 2 (Aschwanden et al. 2016); in the fractal-diffusive model for cellular automata the key variable is the spatial dimension of the cell dynamics, with Γ = 1.5 predicted for 3D models (Aschwanden 2012). This is very close to the typical value Γ = 1.6 for magnetar X-ray bursts (Cheng et al. 1996;Gögüş et al. 1999Gögüş et al. , 2000Gavriil, Kaspi, & Woods 2004), which is plotted in Fig 5 for comparison.
In corona-type cellular automata models, the cell rules encode an immediate diffusive redistribution of energy from one cell to its neighbours (Isliker, Anastasiadis, & Vlahos 2000) and lead to a tight N − E correlation (Lu & Hamilton 1991;Dȃnilȃ, Harko, & Mocanu 2015). Here the cell rules have the less direct effect of making energy release from a plastic cell's neighbours more likely rather than guaranteed, and as a result the burst statistics show more deviation from a simple power law. There is also an overrepresentation of 10 40 and 10 43 erg events (corresponding to a single cell's shallow or deep failure, respectively). Within the paradigm presented here, therefore, the degree of scatter of bursts from a power-law relation may encode valuable information about the nature of magnetar crustal failure.

OUTLOOK
The model presented here was constructed based on theoretical considerations, and aims to be a faithful minimal representation of the salient physics of magnetar activity. Where possible, the model is quantitative: the crustal structure is calculated with a realistic equation of state, the elastic stress from a fit to molecular dynamics simulations, and the characteristic lengthscale for the cells taken from 3D magnetoelastic simulations. The details of how the crust fails are, however, unknown, which makes the use of somewhat ad-hoc prescriptions inevitable. Nonetheless, we have performed extensive checks to confirm that the model's key results are robust to variation in cell size (see Fig. 5), the expression for ν, and the critical value at which the transition between elastic and plastic regimes occurs. In particular, we have not attempted to 'tune' any input quantities in order to better mimic the behaviour of any particular magnetar. With this concrete model for crustal failure, however, we can now directly use information from magnetar bursts to constrain the model and the star's physics. For example, because even the largest events in our current model are relatively localised and weaker than the brightest known giant flare, this indicates the presence of an additional mechanism driving the propagation of crustal failure, e.g. a thermoplastic instability (Beloborodov & Levin 2014); this may be manifested observationally as deviations from the power law for weaker bursts. The detection of a giant flare with energy 10 47 erg would point to the involvement of core magnetic field evolution. With a temperature-dependent plastic viscosity, hotter crusts produce frequent bursts, but slightly cooler ones are needed to allow time for a large plastic cluster to form and potentially power a giant flare. If the crust is too cold, however, the sluggish plastic flow is more likely to lead to long-lived multipolar coronal fields or non-dipole spindown than outbursts. It will also be easier to observe FRBs from cooler crusts, since the more sparsely distributed coronal loops will not significantly inhibit radio emission (Wadiasingh & Timokhin 2019;Suvorov & Kokkotas 2019).
Very young and very old magnetars are not expected to produce giant flares: the former because large regions of high stress have not yet developed; the latter because numerous previous events have produced a fractal low-stress region that inhibits further large-scale failure.
Magnetar activity is often linked directly to the evolution of the star's toroidal B, and although results of such simulations also inform our choice of cell geometry, our focus is instead on the distribution of evolving elastic stress in the crust. This evolution is driven by the changing B, but is insensible to its quantitative features. It is not obvious whether activity driven by an intense poloidal B would be noticeably different within our paradigm; perhaps it would lead to a different natural cell geometry. Any other source of stress could also, in principle, drive seismic activity, whether or not it then leads to characteristic magnetar behaviour. The most obvious example would be stresses developing through spindown; in this case both the cell geometry and the evolution equations would need to be revised.
The crust is theéminence grise of magnetars: it powers their activity, but in a way that is difficult to discern from ob-servations, which essentially 'see' only the corona. The goal of this work is to provide a framework to compare observations directly with theory, and so to probe the crust's poorlyunderstood mechanical properties.
It is a pleasure to thank Ersin Gögüş, George Younes and Zorawar Wadiasingh for many stimulating discussions related to this work.