Three-dimensional simulations of the magnetic stress in a neutron star crust
aa r X i v : . [ a s t r o - ph . S R ] J un Three-dimensional simulations of the magnetic stress in a neutron star crust
T. S. Wood ∗ School of Mathematics and Statistics, Newcastle University,Newcastle upon Tyne, NE1 7RU, United Kingdom
R. Hollerbach
Department of Applied Mathematics, University of Leeds, Leeds, LS2 9JT, United Kingdom (Dated: July 18, 2018)We present the first fully self-consistent three-dimensional model of a neutron star’s magnetic field,generated by electric currents in the star’s crust via the Hall effect. We find that the global-scalefield converges to a dipolar Hall-attractor state, as seen in recent axisymmetric models, but thatsmall-scale features in the magnetic field survive even on much longer timescales. These small-scalefeatures propagate toward the dipole equator, where the crustal electric currents organize themselvesinto a strong equatorial jet. By calculating the distribution of magnetic stresses in the crust, wepredict that neutron stars with fields stronger than 10 G can still be subject to starquakes morethan 10 yr after their formation. I. INTRODUCTION
Neutron stars are of interest not only for the exoticstates of matter they contain, but also for their magneticfields, which are the strongest in the universe. In thecase of pulsars, the rotation of the magnetic field pro-duces beams of non-thermal radiation that can be de-tected thousands of light-years away. The field can alsoaffect the dynamics of the star itself, through the mag-netic stress it exerts on the star’s solid outer layer, or“crust”. In neutron stars with especially strong magneticfields (known as magnetars) the magnetic stress can bestrong enough to fracture the crust, producing a star-quake [1]. This is believed to be the mechanism behindthe gamma-ray flares and X-ray outbursts detected inthese objects [2].The magnetic fields of pulsars, inferred from spindownmeasurements, can be anywhere in the range 10 –10 G.However, the spindown rate depends only on the large-scale component of the field at the magnetic poles, sothese measurements may underestimate the actual fieldstrength in the star. This could explain why some neu-tron stars produce magnetar-like emissions yet have slowspin-down rates [3–5]. There is also evidence that at leastsome neutron stars have stronger magnetic features onsmaller scales [6–9]. To interpret these observations it isnecessary to develop a self-consistent, three-dimensional(3D) model of neutron star magnetic fields.The external field of a neutron star is generated byelectric currents flowing within its crust and core. Theions in the crust form a rigid lattice, and the currentsthere arise purely through the flow of electrons, whosedynamics depend primarily on the Hall effect [10]; thissituation is commonly referred to as electron magneto-hydrodynamics (EMHD). Recently there have been nu-merous studies of EMHD in neutron star crusts [11–19]. ∗ [email protected]
However, because of difficulties solving the EMHD equa-tions computationally in 3D, in all of these studies themagnetic field was assumed to be axisymmetric. Al-though there have been 3D studies of EMHD in simpli-fied periodic-box geometry [20–22], the results cannot beused to predict the global-scale field morphology in a realneutron star. Currently, the only predictions regardingthe global-scale field come from axisymmetric models.Using one such axisymmetric model, Gourgouliatosand Cumming [18] found that an initially dipolar mag-netic field evolves towards a quasi-steady configurationthat they called a “Hall attractor”. After ∼ yr theelectric currents in the crust are concentrated in a narrowjet around the dipole equator, producing a strong belt ofpoloidal field. However, it is unknown whether this ax-isymmetric attractor would be stable to 3D perturbationsin a real neutron star.Here, we present for the first time global 3D numericalsimulations of the magnetic field in a neutron star crust.We use a pseudo-spectral code that allows us to performsimulations that are not only fully 3D, but also higherresolution than any of those previously presented even in2D. We describe how the Hall attractor is modified in 3D,and discuss the magnitude and distribution of magneticstresses within the crust. II. THE MODEL
We work in a reference frame that corotates with theneutron star’s crust. Because the ions are fixed withinthe crust, in this frame the electric current, J , dependsonly on the electron fluid velocity, v . The current is alsodirectly related to the magnetic field, B , via Amp`ere’sLaw, and we have (in Gaussian cgs units) J = − e n v = c π ∇ × B , (1)where n is the electron number density, e is the elemen-tary charge, and c is the speed of light. If the crust iseither sufficiently cool, or close to isothermal, then themagnetic field is frozen to the electron fluid [10], andevolves according to the equation ∂ B ∂t = ∇ × (cid:18) c B π e n × ( ∇ × B ) (cid:19) − ∇ × ( η ∇ × B ) , (2)where η is the magnetic diffusivity. The two terms on theright-hand side of Equation (2) represent the Hall effectand Ohmic dissipation, respectively.In this study we approximate the structure of the crustas fixed and spherically symmetric, and so the electrondensity n and magnetic diffusivity η are fixed functionsof radius r . We do not therefore take account of anydeformations in the crust arising from magnetic stresses.However, we can use our results to calculate the magneticstress and determine whether it would be large enoughto induce fractures in the crust. We take the top andbottom of the crust to be the spherical surfaces r = R and r = 0 . R respectively, where R = 10 km is the radiusof a typical neutron star. For n ( r ) and η ( r ) we adoptthe same analytical profiles used by Gourgouliatos andCumming [23]. Specifically, n = n (cid:18) − r/R . (cid:19) and η = η (cid:18) − r/R . (cid:19) − / , where n = 2 . × cm − and η = 4 . × − cm s − are the values at r = R . These profiles are only roughapproximations to the (highly uncertain) profiles in realneutron stars, but fortunately our results are not sensi-tive to the specific profiles used.Finally, we must impose boundary conditions on themagnetic field at the top and bottom boundaries of thecrust. At the top we impose vacuum boundary condi-tions, i.e., we match to a current-free field outside thestar. At the bottom we impose either the idealizedboundary conditions used by Gourgouliatos and Cum-ming [23], which are B r = 0 and J r = 0, or the morerealistic boundary conditions used by Hollerbach andR¨udiger [12], which model the star’s core as a type-I su-perconductor. From here on we refer to these two setsof conditions as GC and HR respectively. In order toimplement the HR boundary conditions, it is convenientto make a small modification to the electron density pro-file n ( r ), as described by Hollerbach and R¨udiger [12], tomake 1 /n vanish at the bottom of the crust. We there-fore use a slightly modified density profile ˜ n ( r ), definedas 1 / ˜ n ( r ) = 1 /n ( r ) − /n (0 . R ).The relative importance of the Hall effect and Ohmicdissipation terms in Equation (2) depends on thestrength of the magnetic field, B say, and is quantifiedby the Hall parameter, H ≡ cB / (4 π e n η ). In a typicalmagnetar, with B = 10 G, we have H ≃
50, implyingthat the Hall effect dominates the dynamics of the mag-netic field. In that case we expect the field to evolve onthe Hall timescale, t Hall ≡ π e n R cB ≃ . , (3) where we have assumed that the characteristic length-scale for the magnetic field is R = 10 km. If the field hasstrong, small-scale features, then these will evolve on ashorter timescale.To solve the EMHD equation (2) we have adapted the3D MHD code PARODY, developed by Dormy et al. [24]and Aubert et al. [25]. The code is pseudo-spectral, anduses spherical harmonic expansions in latitude and lon-gitude, and a discrete grid in radius, making it perfectlysuited to solving problems in spherical-shell geometry.The results that we present here have a resolution of 128grid-points in radius, and spherical harmonics up to de-gree l = 100. We have benchmarked the code againstpreviously published axisymmetric results [12, 23, 26] andfind excellent agreement in all cases. III. RESULTSA. Robustness of the Hall attractor
To determine whether the “Hall attractor” seen in ear-lier axisymmetric simulations is robust against 3D per-turbations, we have repeated one simulation of Gourgou-liatos and Cumming [23], using the same boundary con-ditions and dipolar “Ohmic eigenmode” initial conditionfor the axisymmetric component of the magnetic field. Tothis initial axisymmetric poloidal field we add a low am-plitude, small-scale 3D perturbation with both poloidaland toroidal components, comprising spherical harmon-ics of degree 20 l
40. As shown in Figure 1, the mag-netic field evolves towards a state broadly similar to theaxisymmetric Hall attractor. However, significant non-axisymmetric features persist in the simulation even onlong timescales (comparable to the global Hall timescale t Hall ). These features would rapidly decay by the actionof Ohmic dissipation alone, so their longevity can be at-tributed to the Hall effect. In fact it has previously beenshown that, in the presence of a strong density gradient,the Hall effect can sustain or even amplify small-scale fea-tures in the magnetic field [27]. The length-scale of thelongest-lived non-axisymmetric features is comparable tothe largest scales present in the initial 3D perturbations(degree l = 20), suggesting that they form by an up-scaletransfer of magnetic energy.The convergence of the global-scale magnetic field tothe attractor state occurs on the Hall timescale (3), afterwhich the evolution proceeds on the slower timescale ofOhmic diffusion. This is in agreement with results fromaxisymmetric simulations [12, 18]. (In the simulations ofMarchant et al. [19], by contrast, the field was initializedalready in a Hall equilibrium state, and therefore evolvedon the Ohmic timescale throughout.) The poloidal mag-netic flux becomes concentrated into a “belt” around theequator of the magnetic dipole, implying a strong equa-torial jet of electrons within the crust. The angular ve-locity in this jet is approximately constant along eachpoloidal field line (see Figure 1). We do not find that the PSfrag replacements t = 0 . t = 0 . t = 1 . PSfrag replacements t = 0 . t = 0 . t = 1 . PSfrag replacements t = 0 . t = 0 . t = 1 . -0.6-0.4-0.2-0.00.20.40.6 PSfrag replacements t = 0 . t = 1 . B r / G -0.15-0.10-0.05-0.000.050.100.15 PSfrag replacements t = 0 . t = 1 . B r / G(a) (b)
FIG. 1. (a) Lines of the azimuthally averaged poloidal field and electron angular velocity at successive times, illustratingconvergence to the Hall attractor. The dashed lines indicate the boundaries of the crust. (b) The radial component of themagnetic field at the surface at early and late times in the same simulation.
Hall effect enhances the diffusion of the magnetic field, asoriginally suggested by Goldreich and Reisenegger [10].However, by concentrating the poloidal flux around theequator, the Hall effect does lead to a more rapid declinein the strength of the radial field at the poles.The magnetic field is primarily poloidal, with a weakertoroidal component in the inner crust. However, the GCboundary conditions used in this simulation impose zerotoroidal field both at the top and bottom of the crust.In a real neutron star, currents can flow between thecrust and core, generating significant toroidal fields. Wehave therefore repeated this simulation using the morerealistic HR boundary conditions, which allow for finitetoroidal field at the lower boundary. Figure 2 comparesthe azimuthally averaged magnetic field at the same timein both simulations, revealing that although the mor-phology of the poloidal field is modified, a similar Hallattractor state still exists. The toroidal field in the sec-ond simulation is much stronger, as expected, but re-mains weaker than the poloidal field. This toroidal fieldis associated with a meridional flow of electrons that isequatorward at the star’s surface, and drags the poloidalmagnetic field lines. As a result, the equatorial belt ofpoloidal flux, and the corresponding electron jet, is evenstronger in the simulation with HR boundary conditions,and is pushed slightly deeper into the crust. Figure 3illustrates the development of the equatorial belt, whichmanifests as bands of positive and negative B r on eitherside of the equator at the star’s surface. B. The effect of field strength
Spin-down measurements of neutron stars tell us onlythe strength of the large-scale magnetic field near thepoles. However, it is clear from Figure 3 that the fieldstrength elsewhere in the crust can greatly exceed thisobserved value. At the time plotted in Figure 2b, for ex-ample, the strength of the surface magnetic field variesbetween ≃ × G at the dipole axis and ≃ × Gnear the equator. Near the bottom of the crust, the fieldis stronger by a further order of magnitude than at the
PSfrag replacements B φ / G -1.0-0.8-0.6-0.4-0.20.0 PSfrag replacements B φ / G PSfrag replacements B φ / G FIG. 2. Poloidal fieldlines and contours of toroidal field at t = 1 . surface. The Hall effect may therefore be even more sig-nificant in neutron stars than previously thought, and sowe have repeated the same simulation shown in Figure 3with a stronger initial magnetic field (i.e., a larger Hallparameter H ). Because the evolution of the global-scalemagnetic field occurs on the Hall timescale (3), we expectthat a stronger field will converge more rapidly to theHall attractor, before the small-scale magnetic featureshave been dissipated by resistivity. This is confirmed inFigure 4, which shows the surface radial field in threesimulations with increasing magnetic field strengths. Ineach case the field is plotted at time t = 0 . t Hall , whichcorresponds to t ≃
1, 0.5, and 0.25 Myr respectively, andan equatorial jet of electrons has already formed. In thecase with the strongest magnetic field, the jet is broaderand more spatially disordered; at later times in the samesimulation, the jet becomes increasingly laminar, but re-mains broader than in the other simulations.
C. The crustal magnetic stress
If the magnetic shear stress within the crust exceedsthe breaking stress of the ionic lattice then it will inducea crustal fracture. Molecular dynamics models indicatethat the breaking stress is approximately 5 × − ofthe electron pressure, P e [28]. At each point within thecrust, the strongest magnetic shear stress is exerted on t = 0 .
21 Myr t = 0 .
67 Myr t = 1 .
40 Myr
FIG. 3. Evolution of the magnetic field in the simulation with HR boundary conditions. Top row: B r at r = R . Bottom row: B φ at r = 0 . R . Colorbars are adjusted to the maximum values in each plot, and use a linear scale. -0.1 -0.3-0.2 PSfrag replacements B r /B FIG. 4. The surface radial field in three simulations with B = 10 G, 2 × G, 4 × G (i.e., H = 50 , , a surface that makes an angle of 45 ◦ to the local mag-netic field direction, and this stress is exactly equal tothe magnetic pressure, P m . The condition for fracturingis therefore P m & × − P e , where P m = | B | / (8 π )and P e ≃ (3 π n ) / ℏ c/ (12 π ). Fractures are most likelynear the surface of the crust, where the density and pres-sure are lowest, and the energy released in a near-surfacefracture can directly power flares and outbursts. In Fig-ure 5 we plot the ratio of magnetic pressure to breakingstress at the surface of the crust in the same simulationshown in the last panel of Figure 4, revealing that patchesaround the dipole equator are susceptible to crustal frac-turing. These patches are localized in both latitude andlongitude, as a consequence of the three-dimensionalityof the magnetic field. They persist on long timescales, oforder t Hall , and propagate azimuthally in the direction ofthe electron jet.
IV. DISCUSSION
Our results indicate that neutron starquakes are mostcommon in the vicinity of the electron jet within thecrust, which typically forms a ring around the dipole axis.This jet forms on the Hall timescale, which is dependenton the overall strength of the magnetic field. At earliertimes the structure of the magnetic field is dependent onthe initial conditions for the proto-neutron star, whichunfortunately are not well known. The surface field isstrongest in localized patches around the magnetic equa-tor, where the local field strength typically exceeds thatat the poles by an order of magnitude.In our simulations, surface magnetic features migrate
FIG. 5. Lines of the magnetic field generated by currentsin the crust. The coloring indicates the ratio of magneticpressure to breaking stress at the surface of the crust. Aportion of the surface has been cut away to show fieldlinesinside the crust, whose lower boundary is the gray sphere. equatorward as a result of the toroidal magnetic field inthe crust and the associated meridional flow of electrons.This migration could potentially be reversed if a strongtoroidal field of the opposite sign were present, corre-sponding to a poleward flow of electrons at the star’ssurface [29]. There is no obvious process that can gen-erate such a strong toroidal field within the crust, butit could be a remnant from the star’s formation [30] orthe result of toroidal flux expulsion from the core [31].To fully describe the processes that can impart a strongtoroidal field to the crust, it will be necessary to coupleour crust model to a model of the superconducting core.The flow of heat from a cooling neutron star is inhib-ited across magnetic field lines, so the strong concentra-tion of poloidal magnetic flux in the equatorial belt couldtrap heat within the electron jet [32]. The breaking stressis very sensitive to temperature [28], and so the equato-rial region of the crust may be even more susceptible tofracturing than we have found here. We are currentlyextending our model to describe the thermal–magneticevolution of the crust, and its coupling to the core.
ACKNOWLEDGMENTS
T.S.W. and R.H. were supported by STFC grantST/K000853/1. Computations were performed on ARC1and ARC2, part of the High Performance Computing fa-cilities at the University of Leeds. We thank Konstanti-nos Gourgouliatos, C´eline Guervilly, David Ian Jones,Maxim Lyutikov, and the anonymous referees for helpfulcomments and suggestions. [1] M. Ruderman, Astrophys. J. , 587 (1991).[2] C. Thompson and R. C. Duncan,Astrophys. J. , 980 (2001).[3] F. P. Gavriil, M. E. Gonzalez, E. V. Gotthelf,V. M. Kaspi, M. A. Livingstone, and P. M. Woods,Science , 1802 (2008).[4] N. Rea, P. Esposito, R. Turolla, G. L. Israel, S. Zane,L. Stella, S. Mereghetti, A. Tiengo, D. G¨otz, E. G¨o˘g¨u¸s,and C. Kouveliotou, Science , 944 (2010).[5] P. Scholz, V. M. Kaspi, and A. Cumming,Astrophys. J. , 62 (2014).[6] T. G¨uver, E. G¨oˇg¨u¸s, and F. ¨Ozel,Mon. Not. R. Astron. Soc. , 2773 (2011).[7] E. V. Gotthelf, J. P. Halpern, and J. Alford,Astrophys. J. , 58 (2013).[8] A. Tiengo, P. Esposito, S. Mereghetti, R. Turolla,L. Nobili, F. Gastaldello, D. G¨otz, G. L. Israel,N. Rea, L. Stella, S. Zane, and G. F. Bignami,Nature , 312 (2013).[9] U. Geppert, J. Gil, and G. Melikidze,Mon. Not. R. Astron. Soc. , 3262 (2013).[10] P. Goldreich and A. Reisenegger,Astrophys. J. , 250 (1992).[11] D. A. Shalybkov and V. A. Urpin, A&A , 685 (1997).[12] R. Hollerbach and G. R¨udiger,Mon. Not. R. Astron. Soc. , 1273 (2004).[13] J. A. Pons and U. Geppert, A&A , 303 (2007).[14] R. Perna and J. A. Pons, Astrophys. J. , L51 (2011).[15] Y. Kojima and S. Kisaka,Mon. Not. R. Astron. Soc. , 2722 (2012).[16] D. Vigan`o, J. A. Pons, and J. A. Miralles,Computer Physics Communications , 2042 (2012). [17] J. A. Pons, D. Vigan`o, and N. Rea,Nature Physics , 431 (2013).[18] K. N. Gourgouliatos and A. Cumming,Phys. Rev. Lett. , 171101 (2014).[19] P. Marchant, A. Reisenegger, J. Alejandro Valdivia, andJ. H. Hoyos, Astrophys. J. , 94 (2014).[20] D. Biskamp, E. Schwarz, A. Zeiler, A. Celani, and J. F.Drake, Phys. Plasmas , 751 (1999).[21] J. Cho and A. Lazarian, Astrophys. J. , L41 (2004).[22] C. J. Wareing and R. Hollerbach,J. Plasma Phys. , 117 (2010).[23] K. N. Gourgouliatos and A. Cumming,Mon. Not. R. Astron. Soc. , 1618 (2014).[24] E. Dormy, P. Cardin, and D. Jault,Earth Planet. Sci. Lett. , 15 (1998).[25] J. Aubert, J. Aurnou, and J. Wicht,Geophys. J. Int. , 945 (2008).[26] R. Hollerbach and G. R¨udiger,Mon. Not. R. Astron. Soc. , 216 (2002).[27] T. S. Wood, R. Hollerbach, and M. Lyutikov,Phys. Plasmas , 052110 (2014).[28] A. I. Chugunov and C. J. Horowitz,Mon. Not. R. Astron. Soc. , L54 (2010).[29] U. Geppert and D. Vigan`o,Mon. Not. R. Astron. Soc. , 3198 (2014).[30] C. Thompson and R. C. Duncan,Astrophys. J. , 194 (1993).[31] S. K. Lander, Mon. Not. R. Astron. Soc. , 424 (2014).[32] D. Vigan`o, N. Rea, J. A. Pons, R. Perna,D. N. Aguilera, and J. A. Miralles,Mon. Not. R. Astron. Soc.434