Sleeping beasts: strong toroidal magnetic field in quiescent magnetars explains their large pulsed fraction
Andrei P. Igoshev, Rainer Hollerbach, Toby Wood, Konstantinos N. Gourgouliatos
SSleeping beasts: strong toroidal magnetic field inquiescent magnetars explains their large pulsedfraction
Andrei P. Igoshev * , Rainer Hollerbach * , Toby Wood † and Konstantinos N. Gourgouliatos ‡ October 20, 2020
Magnetars are neutron stars (NSs) with extreme magnetic fields of strength × − G. They exhibit transient, highly energetic events, such as shortX-ray flashes, bursts and giant flares, all of which are powered by their enormousmagnetic energy . Quiescent magnetars have X-ray luminosities between and erg / s, and are further classified as either persistent or transient magnetars.Their X-ray emission is modulated with the rotational period of the NS, with atypical relative amplitude (so-called pulsed fraction) between 10-58 per cent, im-plying that the surface temperature is significantly non-uniform despite the highthermal conductivity of the star’s crust. Here, we present the first 3D magneto-thermal MHD simulations of magnetars with strong toroidal magnetic fields. Weshow that these models, combined with ray propagation in curved space-time, ac-curately describe the light-curves of most transient magnetars in quiescence andallow us to further constrain their rotational orientation. We find that the pres-ence of a strong toroidal magnetic field explains the observed asymmetry in thesurface temperature, and is the main cause of the strong modulation of thermalX-ray emission in quiescence. Soft X-ray emission from magnetars in quiescence originates from their surface,either at the top of their solid outer crust or their atmosphere. Magnetic fields deeper inthe crust control the surface temperature distribution and consequently the X-ray emis-sion. The magnetic field provides heating through its Ohmic decay, and also governshow this heat flows through the crust, by inhibiting di ff usion perpendicular to mag-netic field lines. Regions of open field lines cool rapidly, while heat remains trapped inregions of closed field lines . The field is generated by dynamo action during the proto-NS phase, and is expected to have both poloidal and toroidal components , althoughthe energy of the toroidal component could be ten times larger . Only the poloidal * Department of Applied Mathematics, University of Leeds, Leeds LS2 9JT, UK † School of Mathematics, Statistics and Physics, Newcastle University, Newcastle upon Tyne, NE1 7RU,UK ‡ University of Patras, Department of Physics, 26504, Patras, Greece a r X i v : . [ a s t r o - ph . H E ] O c t eld can be measured directly, via the neutron star spin-down , but there is also ob-servational evidence of a strong toroidal field. The toroidal component is responsiblefor magnetospheric twisting and, therefore, the transient behaviour of magnetars . TheX-ray spectra of many magnetars are best described if a large toroidal magnetic field isassumed .When interpreting X-ray observations, the surface thermal pattern resulting frommagneto-thermal evolution is approximated empirically as a collection of circular re-gions with di ff erent temperatures . Originally, these regions were placed at the mag-netic poles of an assumed dipolar field , but such a configuration cannot produce alarge pulsed fraction . Therefore modern interpretations allow for regions that do notcoincide with the magnetic poles, and have varying sizes and temperatures . To ex-plain the formation, location and shape of these hot surface regions requires a detailedthree-dimensional model of the temperature and magnetic field in the crust.Here we investigate for the first time the formation and evolution of hotter andcolder regions at the surface of a quiescent magnetar, using three-dimensional mag-netohydrodynamic (MHD) simulations in a spherical shell performed with a modifiedversion of the PARODY code (see also Methods Section 1). We simulate the magneto-thermal evolution for two field configurations that have strong toroidal fields containing90% of the total magnetic energy: in model A the poloidal and toroidal componentsare aligned, and in model B the toroidal magnetic field is inclined by 45 ◦ with respectto the poloidal dipole. The initial dipole magnetic field is 1 . × G; the maximumvalues of magnetic field in the crust at the beginning of the simulations are 1 . × G.Figure 1 shows the surface temperature distribution for these models after about 20 Kyrof evolution. The isothermal, purely magnetic properties of these models have previ-ously been studied in detail . The filamentary pattern of hot and cold regions visi-ble in Figure 1 reflects the magnetic field structure arising from an instability of thetoroidal field . Both models exhibit north-south asymmetry: model A has a hot zonethat wraps around the dipole axis, whereas model B has a single hot spot. The sizeof these hot zones are consistent with observations of quiescent magnetars by X-rayspectroscopy.We further compute the light-curves produced by each of these models (see Meth-ods Section 2 for details) taking into account relativistic e ff ects. We assume that theNS has radius R =
12 km and mass M = . M (cid:12) . Because magnetars rotate relativelyslowly, we use approximations for ray propagation in the Schwarzschild metric. Wefind that models A and B have soft X-ray luminosities of 0 . − × erg / s and pulsedfraction ranges from 16 to 53 per cent, which is consistent with observations of tran-sient magnetars. By contrast, models that have weak toroidal magnetic fields have atemperature distribution that is very symmetric with respect to the magnetic equator ,and typically have a maximum pulsed fraction of ≈ κ the anglebetween the dipole axis and the rotation axis, i the angle between the observer’s line ofsight and the rotation axis, and ∆Φ the phase shift. We fit our light-curves to the foldedsoft X-ray emission of seven transient magnetars in quiescence with L X (cid:46) erg / s.The details of the observational reduction and the fitting procedure can be found inMethods Sections 3 and 4 respectively. Briefly, we analyse old observations of magne-tars in quiescence and produce period-folded light-curves in the soft X-ray range 0.3-22 s T s Fig.
1: Thermal maps obtained in 3D magneto-thermal simulations. Left panel: modelA for NS with aligned poloidal and toroidal magnetic fields, age 18 Kyr. Right panel:model B for NS with poloidal and toroidal magnetic fields inclined by angle of 45 ◦ ,age 24 Kyr. The surface temperatures are in units of MK.Source name κ i ∆Φ Age Model χ / d.o.f. L . − x ( ◦ ) ( ◦ ) ( ◦ ) (Kyr) 10 erg / sSGR 0418 + ±
26 274 ±
22 217 ± /
13 0.00771E 1547.0-5408 106 ± ± ± /
13 19CXOU J164710.0-455216 206 ±
33 69 ±
22 32 ± /
13 5.5XTE J1810–197 153 ± ± ± /
13 5.8Swift J1822.3-1606 193 ±
12 284 ±
13 217 ± /
13 0.81SGR 0501 + /
13 3 . + / ∼ . + N H is unknown, so the X-ray luminosity is indicative of typical N H .KeV. Our fits are weakly sensitive to the assumed radius and mass of the NS.The parameters that produce the best fit in each case are summarised in Table 1.For the four magnetars SGR 0418 + . We briefly discuss two cases that are described less successfully by ourmodel in Methods Section 5. Our estimate for the angle κ = ◦ for 1E 1547.0-5408is somewhat di ff erent from the value inferred from radio observations κ = ◦ . Thisdi ff erence could be caused by a non-dipolar magnetic field with complicated magneticstructure. Alternatively, fitting against model B at the age 24 Kyr gives a very di ff erent3ngle, κ ≈ ◦ ( χ = C o un t s ModelObservations 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00Phase50060070080090010001100 C o un t s ModelObservations
Fig.
2: Folded soft X-ray light-curve (300-2000 eV) for magnetars. Left panel: SGR0418 + σ confidence intervals. The solid black lines are the theoreticallight-curve for the most favourable orientation.Our model only describes the magnetic field evolution in the NS crust. Furtherwork is needed to better understand the magnetic field evolution in the NS core, whereambipolar di ff usion might play an important role , particularly in very young NSs.Whether field evolution in the core is significant for quiescent magnetars is unknown.The toroidal magnetic field, which is the main source of the magnetospheric twist,could also cause crust yielding. We have shown that the same toroidal magnetic fieldcan naturally explain the X-ray emission of quiescent magnetars. Therefore, possiblythe main di ff erence between a magnetar and a strongly magnetised neutron star whichshows no magnetar-like behaviour is the strength of the toroidal magnetic field in thecrust. With the revolutionary insight obtained by the NICER telescope for recycledpulsars , it is becoming increasingly clear that the magnetic field structure of NSsis complicated, so it is extremely important to explore the process of magnetic fieldevolution and formation for NSs.In summary, our results provide strong support that an intense crustal toroidal fieldis an essential ingredient, not only for the magnetospheric behaviour, but also the ther-mal radiation originating from the crust. As well as providing heat through Ohmicdecay, it is also responsible for the formation of thermal spots. Our simulation resultsnot only produce qualitative agreement with the observational data, but also provideconstraints on the strength and geometry of magnetar magnetic fields.4 ethods1 MHD and Thermal simulations We integrate the two coupled equations describing magnetic induction and heat transferwithin the NS crust: ∂(cid:126) B ∂ t = − c ∇ × (cid:40) π en e ( ∇ × (cid:126) B ) × (cid:126) B + c πσ ∇ × (cid:126) B − e S e ∇ T (cid:41) , (1) C V ∂ T ∂ t = ∇ · ( k · ∇ T ) + |∇ × (cid:126) B | c π σ + (cid:18) c π e (cid:19) T ∇ S e · ( ∇ × (cid:126) B ) . (2)Here (cid:126) B is the magnetic field, T is the temperature, c is the speed of light, e is the ele-mentary charge, n e is the electron density, S e is the electron entropy, σ is the electricalconductivity, C V is the crust heat capacity, and k is the thermal conductivity tensor. Weuse the equation of state for a degenerate, relativistic Fermi gas, and the Wiedemann–Franz law: S e = (cid:32) π n e (cid:33) / k B Tc (cid:125) , and ( k − ) i j = e π k B T (cid:32) σ δ i j + ε i jk B k ecn e (cid:33) (3)where k B is Boltzmann’s constant, and (cid:125) is Planck’s constant.The induction equation (1) describes the evolution of the magnetic field due to theHall e ff ect, Ohmic decay, and the Biermann battery. Our previous work includedonly this equation, and without the Biermann battery term. The heat equation (2), in-cluded here for the first time, describes the evolution of temperature due to anisotropicheat di ff usion, Ohmic heating, and electron entropy advection. In both equations thefinal term is generally small, but is included for completeness. On the timescales ofinterest the heat capacity of the crust is negligible, but for numerical convenience weinclude a small heat capacity C V that is proportional to σ T . We adopt the density andconductivity profiles used in .Equations (1) and (2) are solved within a spherical shell with 9 km < r <
10 kmusing the pseudo-spectral code
PARODY . We use 128 numerical cells in the radialdirection and spherical harmonics up to degree l = ff usion, and Adams–Bashforth for the remaining terms. We use vacuum bound-ary conditions for the magnetic field at the upper boundary, and perfectly conductingboundary conditions at the lower boundary, assuming for simplicity that all magneticflux is expelled from the core. The upper boundary condition for the temperature is thestandard thermal-blanket relation − (cid:126) r · k · ∇ T | b = σ S T s (4)where σ S is the Stefan–Boltzmann constant. We employ a simple relation between thesurface temperature T s , and the temperature at the top of the crust T b : (cid:18) T b K (cid:19) = (cid:18) T s K (cid:19) (5)5he core is assumed to have a fixed temperature of 10 K.The model physics is simplified in two respects: (1) We neglect any cooling byneutrinos, both in the core and in the crust. Neutrino cooling is important for thelong-term temperature evolution, and for bursting behaviour, but is less relevant toquiescent emission in young magnetars. (2) The electrical conductivity is assumed tobe independent of temperature. In the outer part of the crust the conductivity is knownto depend on temperature, but we note that the magnetic field evolution is primarilydetermined by the Hall term rather than by the conductivity. These limitations will beaddressed in future work.
To compute the corresponding light-curve from a thermal map we use a numericalmethod with angles i and κ , where i is the angle between rotational axis and line ofsight, and κ is the angle between the original magnetic dipole and rotational axis. Co-ordinates at the NS surface are computed with respect to the magnetic pole as θ, φ . Thisis di ff erent from where the hot spots are assumed to coincide with magnetic poles.This is not the case in our simulations, where hot regions are extended and located ata significant separation from magnetic poles. In a few cases we tried to optimise theNS radius and mass as well, but due to the low photon counts (maximum 10 ) andslow rotation of magnetars the light-curve only depends weakly on the exact valuesof NS compactness. We therefore kept these parameters fixed during the optimisationprocess.We convert the temperature obtained using the upper boundary condition to in-tensity of X-ray emission from a particular element at the NS surface using a simpleblackbody model. We use a beaming factor proportional to cos α , where α is the an-gle between the direction where a photon is emitted and normal to the surface at theemission point. This curve roughly follows the numerical beaming function takinginto account vacuum polarisation e ff ects.To produce the light-curve, we integrate the flux which reaches the observer overthe whole visible hemisphere for each rotational phase. We normalise the light-curveby mean luminosity of the source seen for this particular orientation. We provide the observational IDs of dataset for magnetars in quiescence in Table 2;these are old observations . To analyse the Chandra observations we use the softwarepackage CIAO 4.12 together with the calibration database CALDB 4.9.0. The obser-vations are reprocessed with help of chandra repro package. During the analysis theMcGill magnetar catalogue was used extensively . Only events from a region centredat the source (according to the catalogue) with radius of 4 (cid:48)(cid:48) were extracted. Becausewe are interested in thermal quiescence emission, we filter out all photons outside ofthe 300-2000 eV energy interval. All times of arrival for events are transformed to http: // / pulsar / magnetar / main.html / mode Obs IDSGR 0418 + Chandra / TE 13148, 13235, 13236SGR 0501 + Chandra / TE 14811, 155641E 1547.0–5408
XMM Newton / PN 0604880101CXOU J164710.0–455216
XMM Newton / PN and MOS 0404340101XTE J1810–197
Chandra / TE 13746, 13747, 15870, 15871Swift J1822.3–1606
Chandra / TE 14819, 15988, 15989, 15992, 159933XMM J185246.6 + XMM Newton / MOS 0550671301, 0550671801, 0550671901Table 2: Data sets analysedthe solar system baricentre using axbary tool together with the
DE-405 solar systemephemeris and orbital information provided by the
Chandra data archive. We alsovisually inspected source and background light-curve to verify an absence of flares.We search for the magnetar period using the fast Fourier transform and period-folding ( pfold package) for each individual observation and compared with ephemeriscomputed based on measurements of period and period derivative collected by di ff erentauthors. If the rotational period is not seen in a particular observation, we disregard thisdataset. If an observational period is hard to determine to four significant digits fromindividual observation, we use the ephemeris value. After this a folded light-curve with16 phase-bins is produced.The first folded light-curve is phase-shifted to place minimum photon count atphase 0. If the magnetar was observed multiple times, the following folded light-curvesare produced following exactly the same procedure, but at the last step the phase-shiftbetween di ff erent observations is determined using correlation function. The result-ing light-curve is produced by summation of total number of photons in bins seen indi ff erent observations taking into account the phase-shift.Working with the XMM-Newton observations we use heasoft
SAS
RATE<0.4 for energy range 10-12 KeV. We further extract events with energies in the range 300-2000 eV centred at the source position with extraction radius of 20 arcsec. Only singleand double photon events
PATTERN<=4 for PN and PATTERN<=12 for
MOS1 and
MOS2 are selected at this stage. All arrival times are transformed to the baricentre of the solarsystem using barycen task. We prefer to analyse the PN observations, but if a smallnumber of photons is registered, we also added results from both MOS1 and
MOS2 cam-eras. As also noted in , in the case of 3XMM J1852, we had to rely only on MOS1 and
MOS2 observations. When the light-curve is extracted, we follow the same procedureas in the case of the
Chandra data and sum counts in individual phase bins, taking intoaccount possible phase shift between observations.The thermal X-ray luminosities are estimated in the spectral range 0 . − srcflux program with the mean photon energy 1 . N H values from the McGill catalogue. In the case of XMM-Newton observations we used xspec to analyse the spectra and flux and cflux command to estimate the flux. 7
Statistical analysis
After we obtain an observational folded X-ray light-curve, we perform optimisation ofthe model searching for the most probable values of three continuous parameters κ , i and ∆Φ . To do so, we use the maximum likelihood technique with likelihood in form ofC-statistics . The optimum value is found using the Nelder-Mead algorithm . Whenthe most probable values are found, we try thermal maps produced for alternative modeland for later ages and perform optimisation again. We choose the model and age whichcorrespond to the lowest value of the C-statistics. We additionally check the quality ofthe final fit using the χ test. The confidence intervals are computed for each parameter κ, i and ∆Φ by fixing the other two parameters and searching for a new value of χ statistics which di ff ers from original value by 3.84 (95% probability for χ with a singlevariable). In two cases, SGR 0501 and 3XMM J185246.6 + .In the case of 3XMM J185246.6 + N H value is unknown, therefore this object couldbe even brighter than 2 × erg / s, so it might be a persistent magnetar. Data Availability Statement
The data that support the plots within the paper and other findings are availablefrom the corresponding authors upon reasonable request.
Code Availability Statement
The codes that were used to prepare our models within the paper are available fromthe corresponding authors upon reasonable request.
References [1] Victoria M. Kaspi and Andrei M. Beloborodov. Magnetars.
Ann. Rev. Astron.Astroph. , 55(1):261–301, Aug 2017.8 .00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00Phase200300400500600 C o un t s ModelObservations 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00Phase20406080100120140160 C o un t s ModelObservations
Fig.
3: Folded soft X-ray light-curve (300-2000 eV) for magnetars. Left panel:SGR0501, right panel 3XMM J185246.6 + σ confidence intervals.[2] Christopher Thompson and Robert C. Duncan. The Soft Gamma Repeaters asVery Strongly Magnetized Neutron Stars. II. Quiescent Neutrino, X-Ray, andAlfven Wave Emission. Astrophys. J. , 473:322, Dec 1996.[3] Jos´e A. Pons and Daniele Vigan`o. Magnetic, thermal and rotational evolutionof isolated neutron stars.
Living Reviews in Computational Astrophysics , 5(1):3,December 2019.[4] Lilia Ferrario, Andrew Melatos, and Jonathan Zrake. Magnetic Field Generationin Stars.
SSRv , 191(1-4):77–109, October 2015.[5] Jonathan Braithwaite and Hendrik C. Spruit. A fossil origin for the magnetic fieldin A stars and white dwarfs.
Nature , 431(7010):819–821, October 2004.[6] J. Braithwaite and Å. Nordlund. Stable magnetic fields in stellar interiors.
Astron.Astrophys. , 450(3):1077–1095, May 2006.[7] James E. Gunn and Jeremiah P. Ostriker. Magnetic Dipole Radiation from Pul-sars.
Nature , 221(5179):454–456, February 1969.[8] C. Thompson, M. Lyutikov, and S. R. Kulkarni. Electrodynamics of Magne-tars: Implications for the Persistent X-Ray Emission and Spin-down of the SoftGamma Repeaters and Anomalous X-Ray Pulsars.
Astrophys. J. , 574(1):332–355, July 2002.[9] Maxim Lyutikov and Fotis P. Gavriil. Resonant cyclotron scattering and Comp-tonization in neutron star magnetospheres.
Mon. Not. Roy. Astron. Soc. ,368(2):690–706, May 2006.[10] N. Rea, S. Zane, R. Turolla, M. Lyutikov, and D. G¨otz. Resonant CyclotronScattering in Magnetars’ Emission.
Astrophys. J. , 686(2):1245–1260, October2008. 911] Chin-Ping Hu, C. Y. Ng, and Wynn C. G. Ho. A systematic study of soft X-ray pulse profiles of magnetars in quiescence.
Mon. Not. Roy. Astron. Soc. ,485(3):4274–4286, May 2019.[12] D. Vigan`o, N. Rea, J. A. Pons, R. Perna, D. N. Aguilera, and J. A. Miralles. Unify-ing the observational diversity of isolated neutron stars via magneto-thermal evo-lution models.
Mon. Not. Roy. Astron. Soc. , 434(1):123–141, September 2013.[13] Feryal ¨Ozel, Dimitrios Psaltis, and Victoria M. Kaspi. Constraints on ThermalEmission Models of Anomalous X-Ray Pulsars.
Astrophys. J. , 563(1):255–266,December 2001.[14] T. S. Wood and R. Hollerbach. Three Dimensional Simulation of the MagneticStress in a Neutron Star Crust.
Phys. Rev. Lett. , 114(19):191101, May 2015.[15] Konstantinos N. Gourgouliatos and Rainer Hollerbach. Magnetic Axis Drift andMagnetic Spot Formation in Neutron Stars with Toroidal Fields.
Astrophys. J. ,852(1):21, January 2018.[16] K. N. Gourgouliatos and Jos´e A. Pons. Nonaxisymmetric Hall instability: A keyto understanding magnetars.
Physical Review Research , 1(3):032049, December2019.[17] Konstantinos N. Gourgouliatos, Rainer Hollerbach, and Robert F. Archibald.Modelling neutron star magnetic fields.
Astronomy and Geophysics , 59:5.37–42,2018.[18] F. Camilo, J. Reynolds, S. Johnston, J. P. Halpern, and S. M. Ransom. The Mag-netar 1E 1547.0-5408: Radio Spectrum, Polarimetry, and Timing.
Astrophys. J. ,679(1):681–686, May 2008.[19] Andrea Passamonti, Taner Akg¨un, Jos´e A. Pons, and Juan A. Miralles. On themagnetic field evolution time-scale in superconducting neutron star cores.
Mon.Not. Roy. Astron. Soc. , 469(4):4979–4984, August 2017.[20] T. E. Riley, A. L. Watts, S. Bogdanov, P. S. Ray, R. M. Ludlam, S. Guillot,Z. Arzoumanian, C. L. Baker, A. V. Bilous, D. Chakrabarty, K. C. Gendreau,A. K. Harding, W. C. G. Ho, J. M. Lattimer, S. M. Morsink, and T. E. Strohmayer.A NICER View of PSR J0030 + Astrophys. J. Lett. , 887(1):L21, December 2019.[21] A. V. Bilous, A. L. Watts, A. K. Harding, T. E. Riley, Z. Arzoumanian, S. Bog-danov, K. C. Gendreau, P. S. Ray, S. Guillot, W. C. G. Ho, and D. Chakrabarty.A NICER View of PSR J0030 + Astrophys. J. Lett. , 887(1):L23, December 2019.[22] Konstantinos N. Gourgouliatos, Toby S. Wood, and Rainer Hollerbach. Magneticfield evolution in magnetar crusts through three-dimensional simulations.
Pro-ceedings of the National Academy of Science , 113(15):3944–3949, April 2016.1023] Emmanuel Dormy, Philippe Cardin, and Dominique Jault. MHD flow in a slightlydi ff erentially rotating spherical shell, with conducting inner core, in a dipolarmagnetic field. Earth and Planetary Science Letters , 160(1-2):15–30, July 1998.[24] Julien Aubert, Jonathan Aurnou, and Johannes Wicht. The magnetic structureof convection-driven numerical dynamos.
Geophysical Journal International ,172(3):945–956, March 2008.[25] E. H. Gudmundsson, C. J. Pethick, and R. I. Epstein. Structure of neutron starenvelopes.
Astrophys. J. , 272:286–300, September 1983.[26] Andrei M. Beloborodov. Gravitational Bending of Light Near Compact Objects.
Astrophys. J. Lett. , 566(2):L85–L88, Feb 2002.[27] M. van Adelsberg and D. Lai. Atmosphere models of magnetized neutron stars:QED e ff ects, radiation spectra and polarization signals. Mon. Not. Roy. Astron.Soc. , 373(4):1495–1522, Dec 2006.[28] S. A. Olausen and V. M. Kaspi. The McGill Magnetar Catalog.
Astrophys. J.Supp. Ser. , 212(1):6, May 2014.[29] W. Cash. Parameter estimation in astronomy through application of the likelihoodratio.
Astrophys. J. , 228:939–947, Mar 1979.[30] J. A. Nelder and R. Mead. A Simplex Method for Function Minimization.
TheComputer Journal , 7(4):308–313, 01 1965.
Correspondence
Correspondence should be addressed to Andrei Igoshev and Rainer Hollerbach.
Acknowledgements
This work was supported by STFC grant No. ST / S000275 /
1. The numerical sim-ulations were carried out on the STFC-funded DiRAC I UKMHD Science Consortiamachine, hosted as part of and enabled through the ARC3 HPC resources and supportteam at the University of Leeds.
Contributions
All authors contributed to the simulation design, interpretation and writing themanuscript. A.P.I. carried out the X-ray data reduction, the MHD simulations and themodel fitting. T.S.W. adapted the PARODY code to solve magneto-thermal equations.