3D Magneto-thermal Simulations of Tangled Crustal Magnetic Field in Central Compact Objects
Andrei P. Igoshev, Konstantinos N. Gourgouliatos, Rainer Hollerbach, Toby S. Wood
DDraft version January 22, 2021
Typeset using L A TEX twocolumn style in AASTeX63
3D Magneto-thermal Simulations of Tangled Crustal Magnetic Field in Central Compact Objects
Andrei P. Igoshev , Konstantinos N. Gourgouliatos , Rainer Hollerbach , and Toby S. Wood Department of Applied Mathematics, University of Leeds, LS2 9JT Leeds, UK Department of Physics, University of Patras, 26504 Patras, Greece School of Mathematics, Statistics and Physics, Newcastle University, Newcastle upon Tyne, NE1 7RU, UK (Received June 1, 2019; Revised January 10, 2019; Accepted January 22, 2021)
Submitted to AJABSTRACTCentral compact objects are young neutron stars emitting thermal X-rays with bolometric lumi-nosities L X in the range 10 – 10 erg/s. Gourgouliatos, Hollerbach and Igoshev recently suggestedthat peculiar emission properties of central compact objects can be explained by tangled magneticfield configurations formed in a stochastic dynamo during the proto-neutron star stage. In this casethe magnetic field consists of multiple small-scale components with negligible contribution of globaldipolar field. We study numerically three-dimensional magneto-thermal evolution of tangled crustalmagnetic fields in neutron stars. We find that all configurations produce complicated surface thermalpatterns which consist of multiple small hot regions located at significant separations from each other.The configurations with initial magnetic energy of 2 . − × erg have temperatures of hot regionsthat reach ≈ . ≈ . ≤ −
11 %.Therefore, the tangled magnetic field configuration can explain thermal emission properties of somecentral compact objects.
Keywords: magnetohydrodynamics — stars: neutron — stars: magnetic field — methods: numerical INTRODUCTIONCentral compact objects (CCOs) are neutron stars(NSs) found in close proximity to the geometric center ofsupernova remnants (SNRs) (Pavlov et al. 2004; de Luca2008; De Luca 2017). CCOs have soft X-ray luminositiesof 10 − erg/s, incompatible with cooling luminosi-ties of young NSs with ages of a few kyr. Their X-rayspectra generally contain two black-body components:the first arising from bulk NS thermal emission withtemperature ≈ Corresponding author: Andrei P. [email protected], [email protected] netars (with the exception of 1E 161348-5055 in SNRRCW 103, Rea et al. 2016, which might be a magnetarwith unusual formation path). There is no indicationthat they are binaries. CCOs are numerous: up to 30%of NSs associated with SNRs at distances up to 5 kps areCCOs (de Luca 2008), which is difficult to reconcile withthe Galactic NS birthrate (Keane & Kramer 2008). It isunclear how CCOs age. CCO-like objects are not foundamong nearby cooling NSs (Turolla 2009) with ages of >
100 kyr. Bogdanov et al. (2014) and Luo et al. (2015)searched for radio pulsars which show CCO-like X-rayemission and found none.Long-term X-ray observations have revealed the spinperiod and period derivative for three CCOs (Halpern& Gotthelf 2010; Gotthelf et al. 2013). The spin periodsare: RX J0822-4300 in SNR Puppis A has P = 112 ms(Gotthelf & Halpern 2009), CXOU J185238.6+004020in Kes 79 has P = 105 ms (Gotthelf et al. 2005) and1E 1207.4-5209 in PKS 1209-51/52 has P = 424 ms a r X i v : . [ a s t r o - ph . H E ] J a n Igoshev et al. (Zavlin et al. 2000). From these measurements andtheir period derivatives, the magnetic dipole inferredfor these CCOs can be estimated using the formula B p ≈ . × (cid:112) P ˙ P G (Lorimer & Kramer 2004). Theinferred values of B p are in the range 10 – 10 G,which is one or two orders of magnitude smaller thanthe 10 G dipole fields of young radio pulsars. The typ-ical spin-down luminosity ˙ E for CCOs is ∼ erg s − ,which means that magnetospheric currents arising frompulsar spin-down cannot be responsible for formation ofcompact hot regions at the CCOs’ surfaces.CCOs with measured periods and period derivativesare located among old radio pulsars in the P – ˙ P plane.While these pulsars are prominent radio sources, CCOsshow no non-thermal radio emission, neither persistentnor periodic. Coherent pulsar radio emission is stronglybeamed and is generated in the magnetosphere by somekind of plasma instability, or as suggested recently bynon-stationary plasma discharge (Philippov et al. 2020).Even if we miss the beam of radio emission due to thepulsar orientation, the presence of a pulsar wind neb-ula might be expected. Young, rapidly rotating neutronstars embedded in supernova remnants are often sur-rounded by a compact pulsar wind nebula (Gaensler &Slane 2006). This nebula is powered by the particle windgenerated by rotating, magnetised NSs. No evidence ofPWN is found in relation to CCOs.Therefore, observational properties of CCOs pose mul-tiple questions for researchers: (1) what is the sourceof additional heating in these NSs, (2) why are theiremitting areas so small, (3) what are the descendantsof CCOs, (4) why is no radio emission detected, and(5) how different CCOs really are from magnetars. Thisis a part of a larger problem with variety of differentNS classes. A grand unification of NSs is a scenariowhich tries to explain the difference between variousclasses of NSs based on properties and evolution of theirmagnetic fields (Kaspi 2010; Igoshev et al. 2014). Inthe framework of this scenario, CCOs might evolve intosome other category of NSs, such as radio pulsars attimescales of 10 – 10 years.Several scenarios have been proposed to answer thesequestions and to explain the origin and evolution ofCCOs. These scenarios attribute the additional sourceof energy to strong hidden magnetic fields of NSs.Among these scenarios are fall-back accretion and astochastic dynamo (Gotthelf & Halpern 2007). Mag-netic fields of NSs are multi-component; they can includeboth poloidal and toroidal parts, as well as small-scalefields. Only the large-scale poloidal dipole magnetic fieldemerging from the surface can be measured using thetiming technique. In the fall-back accretion scenario, some material ex-pelled by the supernova explosion does not have enoughenergy to escape from the system indefinitely, and is sub-sequently accreted back onto the NS (Chevalier 1989).This accretion continues after the NS crust has solidified.New infalling material buries the magnetic field, cover-ing it with a new crustal layer (Bernal et al. 2010). Theburied magnetic field gradually re-emerges through thesurface over time (Ho 2011; Vigan`o & Pons 2012; Bernalet al. 2013). In this scenario, the excessive heating isexplained by the amplification and dissipation of themagnetic field, which is compressed between the super-conducting NS core and new crust formed from fall-backmaterial. The absence of radio emission is explained bythe suppression of the surface dipolar magnetic field,preventing the NS from producing enough plasma in itsmagnetosphere. The same explanation holds for the lackof a pulsar wind nebula. Igoshev et al. (2016) studied theevolution of small-scale magnetic field after a fall-backand found that the field (both small- and large-scale)re-emerges, so these NSs should start operating as nor-mal radio pulsars after 10 – 10 years. It is still unclearif the heat produced by the magnetic field causes a for-mation of small emitting regions at the surface in thisscenario. If a strong fall-back occurs at the newly-bornmagnetar, a CCO-magnetar could be produced. Thisis a CCO which has an extremely strong magnetic fieldtrapped deep in the crust which should occasionally ex-hibit flares.In the stochastic dynamo scenario, the NS is born witha predominantly small-scale magnetic field, with only aweak dipole component. This small-scale field is formedduring the proto-NS stage, which lasts for tens of sec-onds (Pons et al. 1999) and ends with the solidification ofthe crust. During this stage there are two regions wherethe matter is unstable according to the Ledoux convec-tion criterion. One of these regions is located arounddensities of 10 g cm − of proto-NSs (Thompson &Murray 2001). Nagakura et al. (2020) analysed state-of-the-art 3D simulations of supernova explosions for pro-genitor masses ranging from 9 M (cid:12) to 20 M (cid:12) and foundthat convection always develops in proto-NSs. This con-vection zone is an ideal place for a dynamo which couldgenerate magnetic fields in the range 10 – 10 G.The generation of a large-scale magnetic field isthought to require rapid rotation of the proto-NS (seee.g. recent simulations by Raynaud et al. 2020). If theproto-NS rotates slowly, a dynamo could still operate,but producing only a strong small-scale field (Thompson& Murray 2001). Such a field with surface B s ≈ Gis hidden for a distant observer, because it hardly con-tributes to the NS spin-down. A stochastic dynamo
CO II ∼
10 kyr of evolu-tion, this object might start operating as a radio pul-sar, which solves the problem regarding the absence ofCCO descendants. Because in this case CCOs turn intonormal radio pulsars with no unusual properties after aperiod of time.This article is structured as follows. In Section 2 wedescribe the method which we use in our simulations; inSection 3 we present results of three-dimensional simu-lations and corresponding light curves. We discuss allresults and conclude in Sections 4 and 5 respectively. METHODSWe perform simulations in two steps. First we com-pute the surface temperature distribution using an up-dated version of the MHD code
PARODY (Dormy et al.1998; Wood & Hollerbach 2015; Gourgouliatos et al.2016). Next we compute the light curves and pulsedfraction for different orientations of rotating NSs usingthe ray-tracing in general relativity.2.1.
Magneto-thermal evolution
The details of the magneto-thermal code are sum-marised in De Grandis et al. (2020); Igoshev et al. (2020). We solve magnetic induction and thermal diffu-sion equations in the solid crust of NS using the electron-MHD approximation in which only electrons carry elec-tric charge and thermal flux. Basically, we solve twocoupled partial differential equations. The first is theinduction equation: ∂ (cid:126)B∂t = − c(cid:126) ∇ × (cid:26) πen e ( (cid:126) ∇ × (cid:126)B ) × (cid:126)B + c πσ (cid:126) ∇ × (cid:126)B − e S e ∇ T (cid:27) (1)where (cid:126)B is the magnetic field, n e is the electron numberdensity, e is the elementary charge, c is the speed oflight, σ is the electric conductivity, S e is the electronentropy, and T is the temperature. The second is theheat equation: C V ∂T∂t = (cid:126) ∇ · (ˆ k · (cid:126) ∇ T ) + | (cid:126) ∇ × B | c π σ + (cid:16) c πe (cid:17) T (cid:126) ∇ S e · ( (cid:126) ∇× (cid:126)B ) (2)where ˆ k is the thermal conductivity tensor. The first twoterms in our induction eq. (1) are the ones introducedby Goldreich & Reisenegger (1992), and represent theHall evolution and Ohmic dissipation, respectively. Thethird term describes the so-called Biermann battery ef-fect, i.e. magnetic field generation due to temperaturegradients (Blandford et al. 1983). This effect is closelyrelated to electron baroclinicity, which acts as a sourceof electron circulation, and hence magnetic induction.In the heat eq. (2), the first term representsanisotropic thermal diffusion. While heat flows freelyalong magnetic field lines, heat transfer is inhibited inthe direction orthogonal to field lines. The second termrepresents heating by Ohmic field decay, and the thirdterm represents the entropy carried by electric currents.Unless strong temperature gradients are maintained ex-ternally, e.g. by magnetospheric heating (De Grandiset al. 2020), the last terms in eqs. (1,2) are generallyvery small in comparison to other terms, and play anegligible role in the field evolution.We write the thermal conductivity tensor in the fol-lowing form:(ˆ k − ) ij = 3 e π k B T (cid:18) σ δ ij + (cid:15) ijk B k ecn e (cid:19) (3)where δ ij is the Kronecker delta, (cid:15) ijk is the Levi-Civitasymbol, k B is the Boltzmann constant, and B k is thefield component along the k -axis. Igoshev et al.
Electron entropy is related to temperature and chem-ical potential µ ( r ) of a degenerate, relativistic Fermi gasas follows: S e = π k B Tµ (4)Similar to Gourgouliatos & Cumming (2014) we usethe following radial profiles for the chemical potential,electron density and conductivity: µ ( r ) = µ (cid:18) − r . (cid:19) / (5) n e ( r ) = n (cid:18) − r . (cid:19) (6) σ ( r ) = σ (cid:18) − r . (cid:19) / (7)where µ = 2 . × − erg, n = 2 . × cm − and σ = 1 . × s − . For computational convenience weassume that the heat capacity is proportional to the tem-perature in the form (see e.g. Page et al. 2004) C V ∝ σT (8)so that eq. (2) becomes a linear equation in T .2.2. Boundary and initial conditions
For the magnetic induction at the inner boundary, weassume the Meissner condition, i.e. zero radial compo-nent of magnetic field B r = 0 and zero tangential com-ponent of the current J t = 0 at the crust-core boundary,see De Grandis et al. (2020) and Hollerbach & R¨udiger(2004) for more details. Physically this corresponds toa scenario in which the field is expelled from the corebefore the crust solidifies. In reality, the core might stillcontain some magnetic field, although the timescale onwhich the core field evolves is uncertain, see e.g. Gusakovet al. (2020). For the magnetic field at the outer bound-ary we use vacuum boundary conditions, i.e. we assumethat (cid:126) ∇ × (cid:126)B = 0 in the region outside the star. Thiscondition is implemented numerically as: ∂B m p ,l ∂r + l + 1 r B m p ,l = 0 (9)and B m t ,l = 0 (10)where B m p ,l and B m t ,l represent the coefficients of thepoloidal and toroidal magnetic potentials, which are ex-panded in series of spherical harmonics of degree l andorder m . Physically these boundary conditions corre-spond to neglecting any electric currents in the star’stenuous plasma atmosphere. We note that CCOs haveno radio emission, and no pulsar wind nebula was ever detected around a CCO, so their magnetospheres mighthave less plasma than is typical for normal radio pulsars.For the temperature, at the crust core interface weassume no cooling, yielding the boundary condition ∂T /∂t = 0. In reality the core cools at a rate deter-mined by neutrino emission, which is practically insen-sitive to conditions within the crust. Since our focusin this study is on small-scale temperature anomaliesproduced within the crust we neglect the core coolingfor simplicity. At the outer boundary we assume thatthe heat flux out of the computational domain is subse-quently emitted as black-body radiation from the star’ssurface, i.e. − (cid:126)r · ˆ κ · ∇ T | b = σ S T s (11)where σ S is the Stefan-Boltzmann constant, T b is thetemperature at the top of the computational domain,and T s is the effective surface temperature. We assumethat these two temperatures are related by the thermalblanket equation (cid:18) T b K (cid:19) = (cid:18) T s K (cid:19) , (12)similar to a relation suggested by Gudmundsson et al.(1983), but in a more simplified form. This assumptionis necessary in our calculations because the thermal re-laxation timescale of the envelope is much shorter thanany other evolution timescale involved in our simula-tions. As was shown by Gudmundsson et al. (1982),the effective surface temperature depends essentiallyonly on the temperature in deep layers (below densities10 g cm − ) and surface gravity.We use the same magnetic field configurations as inGourgouliatos et al. (2020), so the initial magnetic fieldconsists of several high-order multipoles with 10 ≤ l ≤
20. We choose phases randomly for these fields. Theamount of initial magnetic energy E tot , in this field isset at the beginning of simulations. A dipolar poloidalmagnetic field is added to the stochastic magnetic fieldconfiguration with fixed B dip , value before the simula-tions are started.We compute the same first five models as in Gour-gouliatos et al. (2020). In models 6 and 7 we decreasethe total magnetic energy in the crust by a factor of 4.We make this change because for larger energies we en-counter certain numerical problems. We summarise thedetails of our simulations in Table 1.2.3. Light curves and pulsed fraction
In order to compute the light curves we follow thesame technique as described in Igoshev et al. (2020).Namely, NS orientation is described using three angles: i is the angle between line of sight and rotational axis of CO II Name E tot , B dip , (cid:104) B (cid:105) p max , PF m (erg) (G) (G) %Model 1 2 . × × . × × . × × . × × . × × × × × × Table 1.
Numerical models computed here. p max , is thevalue of total magnetic energy decay exponent at its firstmaximum. PF m is the maximum pulsed fraction computedfor all models at age 3.5 kyr. NS, κ is the angle between magnetic dipole moment androtational axis, and ∆Φ is the initial phase. These threeangles uniquely describe NS orientation at the start ofrotation. To simulate rotation of NS, we use simplifiedequations presented by Beloborodov (2002) combinedwith: µ = sin i [cos Φ cos θ sin κ − sin θ (sin φ sin Φ − cos κ cos φ cos Φ) ]+ cos i [cos κ cos θ − cos φ sin κ sin θ ] (13)where θ and φ are coordinates at the NS surface withrespect to original dipolar magnetic poles, µ = cos ψ isthe direction toward the position with coordinates θ, φ at infinity, and Φ ∈ [0 , π ] is the angular phase. Oureq. (13) reduces to eq. (5) of Beloborodov (2002) if ahot spot is located at the magnetic pole, i.e. θ = 0 = φ .We choose the photon beaming to be proportional tocos θ (cid:48) (DeDeo et al. 2001) because of the following ar-gument. In an atmosphere which is not affected by mag-netic field, the darkening is described by the Hopf func-tion, see e.g. Chandrasekhar (1950). The intensity emit-ted at angle 90 ◦ is ≈ / ◦ the emission is 0.1 of its value at 10 ◦ , andit goes to 0 at 90 ◦ (see figure 15 of van Adelsberg & Lai2006), which is roughly similar to a cos θ (cid:48) behaviour.For CCOs the pulsed fractions are typically small, so itis better to overestimate the pulsed fraction in our mod-els to check if the model has enough predictive power.The flux emitted by each surface element of NS isintegrated over the visible hemisphere (which is slightlymore than half of the star due to the light bending ina strong gravitational field) to produce the visible fluxfor each rotational phase. We use 16 phases to simulatea light curve. For a few orientations we compute the pulsed fraction (PF) as: P F = F max − F min F max + F min , (14)where F max and F min are maximum and minimum fluxesfor a fixed NS rotational orientation. We select the max-imum pulsed fraction among different orientations andshow them in Table 1.In Igoshev et al. (2020), we fitted observed light curvesof magnetars in quiescence to estimate the NS rotationalorientation. This was possible because magnetars seemto have regular large-scale magnetic fields. This is notthe case for CCOs. If their magnetic fields are indeedformed as a result of a stochastic dynamo, the dipo-lar component is weak and changes location relativelyquickly (see Gourgouliatos et al. 2020). Additional loopsof magnetic field are located randomly; therefore, it isnecessary to describe the location of each individual loopof magnetic field to reproduce a thermal map. Thus,there is no use in exact fitting of the light curve, sincemultiple slightly different magnetic field configurationswill result in very similar light curves but slightly differ-ent NS rotational orientations. RESULTSWe computed models 1-3 up to 50 Kyr, and model4,5,6 and 7 up to 100 Kyr. We show the surface temper-ature distributions for models 2, 4 and 6 in Figures 1, 2and 3. The surface temperature distributions in the caseof model 1 and 3 are very similar to model 2. The sur-face temperatures of model 5 and 7 are very similar toones produced in model 4 and 6, respectively, so we donot plot these maps. Basically the surface temperaturedistributions are defined only by the value of the initialmagnetic energy, since the initial field topology is thesame in all models. The initial strength of the dipolarcomponent weakly affects the values of surface temper-atures and the location of hot and cold regions.In all cases the small-scale structure of the magneticfield leads to formation of complicated thermal patternwhich includes alternating small hot and cold regions.Individual hot spots have simple shapes but tempera-tures of groups of spots differs in different parts of NS,making the pattern even more complicated. Overall, at3.5 Kyr the surface temperature map is composed ofa few relatively small hot regions and hot spot associ-ations, and a much more extended colder continuum.We measure the linear sizes of several hot regions in thecase of model 6 at 1.2 kyr. Two brightest hot regionshave linear sizes of 2.2 and 2.6 km respectively, assumingthat the NS radius is 10 km. At later stages of evolution( >
20 kyr) the temperature and size of hot regions seen
Igoshev et al. B r ( G) (a) B r ( G) (b) B r ( G) (c) T s ( K) (d) T s ( K) (e) T s ( K) (f) Figure 1.
Surface magnetic fields (upper row) and surface temperatures (lower row) for model 2. Individual panels correspondto different ages: (a, d) 3.5 Kyr, (b, e) 9.5 Kyr, (c,f) 37.7 Kyr. B r ( G) (a) B r ( G) (b) B r ( G) (c) T s ( K) (d) T s ( K) (e) T s ( K) (f) Figure 2.
Surface magnetic fields (upper row) and surface temperatures (lower row) for model 4. Individual panels correspondto different ages: (a, d) 3.5 Kyr, (b, e) 9.5 Kyr, (c,f) 37.7 Kyr.
CO II B r ( G) (a) B r ( G) (b) B r ( G) (c) T s ( K) (d) T s ( K) (e) T s ( K) (f) Figure 3.
Surface magnetic fields (upper row) and surface temperatures (lower row) for model 6. Individual panels correspondto different ages: (a, d) 3.5 Kyr, (b, e) 6.8 Kyr (c,f) 38 Kyr.
Igoshev et al. T m a x ( K ) E t o t ( e r g ) Figure 4.
Evolution of maximum surface temperature (bluedots) and total magnetic energy (dashed red line) for model4. in all models decays. The hot regions become completelyisolated and cool down.The hot regions are much hotter in the case of model 4and 6 with larger initial magnetic energy than in model2. In the case of model 6, which has the largest magneticfield energy, the hot spots are reaching temperaturesof 2 . × K, which is ≈ .
19 keV. Gotthelf et al.(2013) mentioned that it is necessary for a successfulCCO model to show variations of the temperature by afactor of two over the surface, because this temperaturedifference is seen in observations of RX J0822-4300. Ourmodels 4, 5, 6 and 7 show temperature variations witha factor of 1 . . ≈
20 kyr, the initial total magnetic energy E tot decays from 2 . × erg to ≈ × erg with typicalpower of ≈ × erg/s. In real NSs this power is emit-ted through the surface photon emission and neutrinos.In our simulations this power is partly emitted throughthe surface boundary condition and partly absorbed bythe inner boundary condition.We estimate the instantaneous exponent which de-scribes the decay of total magnetic energy: E ( t ) = (cid:90) E ( k, t ) dk (15)in a way similar to work by Brandenburg (2020): p ( t ) = − d log E ( t ) d log t (16)We show the results of this evolution in Figure 5. Thebehaviour of p ( t ) is not monotonic. It reaches the first maximum with values ranging 0.39-0.58 (values for in-dividual models can be found in Table 1). This levelis reached at different physical times because the Halltimescale differs in these models. After this initial satu-ration, the p ( t ) value briefly declines and starts growingagain reaching values of 1-3 at the end of our simula-tions, see Figure 5. It means that by the end of oursimulations the energy of magnetic field decays expo-nentially due to the Ohmic decay.The values of 0.39-0.58 corresponds to the value 2 / E = E t − / which is typical for helical config-urations.Another way to describe a complicated field using asingle value is by using the root mean square of mag-netic field B rms , similar to what Brandenburg (2020)obtained using the PENCIL code (Brandenburg et al.2020). We compute B rms over the whole volume of theNS crust. We show the evolution of the maximum tem-perature as a function of B rms in Figure 6. We seethat similar values of B rms might correspond to differ-ent maximum temperatures depending on the evolutionstage. On the other hand, the value of emerging dipolar,poloidal magnetic field seems to correlate with B rms , asseen in Figure 7. The initial value of the dipolar compo-nent gets erased quickly by the Hall evolution, and thefield which emerges afterwards is related to B rms . In thecase of model 2, the emergent field is 2 × G. In thecase of models 6 and 7, the dipolar field reaches valuesof ≈ G, quite independently of the initial dipolarmagnetic field.It is interesting to note that the Hall evolution cre-ates vortices similar to turbulence seen in the surfacetemperature maps, especially in Figure 1, age 9.5 kyr.Typically the Hall evolution creates power-spectra some-what similar to turbulence, but with a different slope of l − , see e.g. Wareing & Hollerbach (2009); Goldreich &Reisenegger (1992). When the Hall evolution starts withinitial condition consisting of a strong global dipole, itmostly preserves the dipole component and redistributesa part of its energy into small-scale magnetic fields. Theintensity of these fields is small in comparison to thedipole, and they are not distinguishable as spatial struc-tures, see e.g. Igoshev et al. (2020). In simulation 2 wesuppress the initial dipole and choose the small-scaleharmonics to be strong at the beginning of the simula-tion. As the result we see that the magnetic energy isredistributed among small-scale harmonics in such a way CO II Age (kyr)0.00.51.01.52.02.53.0 p ( t ) Model 2Model 4Model 6Model 7a (a) Age (kyr)10 ( t ) ( e r g ) Model 2Model 4Model 6Model 7a (b)
Figure 5.
Evolution of the instantaneous energy exponent eq. (16) (a) and total magnetic energy (b) for different models.Model 7a corresponds to model 7 in Gourgouliatos et al. (2020). B rms (G)1.01.21.41.61.82.02.2 m a x T s ( K ) timeModel 2Model 4Model 6Model 7 Figure 6.
Evolution of maximum surface temperature as afunction of B rms . that spatial vortices emerge in the thermal map, whichshows that under certain conditions the Hall turbulenceis very similar to regular turbulence. The system evo-lution of model 2 in scales larger than the crust scale-height H is affected by the geometry of the system. Thisroughly corresponds to modes with (cid:96) ≈ πR NS /H , with (cid:96) being the spherical harmonic decomposition mode and R NS the radius of the star. For a value of H ∼ (cid:96) >
30 where the behaviourof the turbulence will be more profound and the fieldwill not be affected by the geometry of the system.Since the hot regions occur at multiple isolated loca-tions at the surface, it means that a few of these regionsare seen simultaneously, independently of the NS orien- B rms (G)10 B p ( G ) timeModel 2Model 4Model 6Model 7 Figure 7.
Evolution of dipolar poloidal component of themagnetic field as a function of B rms . tation. This effect is amplified further when the curvedpath of photons in the NS gravitational field is takeninto account. We show the light curves for models 4and 6 in Figure 8. Despite the factor two difference intemperature between hot and cold regions, the pulsedfraction is limited to 9-11%, similar to real CCOs.Overall, the light curves have a few components whichcorrespond to different hot regions. The shape of lightcurves is simple because the hot regions are compact.There is no use in discussing how much one light curvepeak is higher/lower than another peak, because it is de-fined by the small-scale structure of the magnetic field.Slightly different configurations of the field could pro-duce a very different light curve shape.0 Igoshev et al. F / F m e a n = F / F m e a n = F / F m e a n = (a) F / F m e a n = F / F m e a n = F / F m e a n = (b) Figure 8.
Soft X-ray light curves for model 4 (left panel) and model 6 (right panel) at age 3.5 Kyr. Each panel corresponds toa different value of κ (written in the lower right corner); black dotted line shows i = 30 ◦ , solid blue line corresponds to i = 60 ◦ and red dashed line is for i = 90 ◦ .4. DISCUSSIONThe main caveats of our simulations are as follows:(1) we do not take NS neutrino cooling into account,(2) we use a simplified treatment of the NS atmosphere,(3) we assume a strong beaming of the thermal photonemission and (4) no magnetic field evolution in the coreis assumed. Below we briefly explain how addition ofthese effects could change the results of our simulations.NS cools down and the core temperature decays be-low 10 K assumed in our simulations within first 10 kyrof the NS life. Therefore, the surface temperatures willbe smaller in extended simulations. Nevertheless, themagnetic field isolates certain regions from the core andheats them up, so the temperature difference betweenhot and cold regions might increase in realistic simu-lations with NS neutrino core cooling. Another factorwhich is not taken into account is the neutrino coolingof the crust. It might remove the excess heat from hotspots, limiting their temperatures by some threshold.In our simulations we assume a temperature depen-dence between the deep crust and surface in form ofeq. (12). In reality the temperature depends on mag-netic fields as well, see e.g. Potekhin & Yakovlev (2001).This might change the surface temperature thermal mapand could affect the light curves. We assume a strongbeaming of the thermal emission proportional to cos θ (cid:48) ,which approximately follows the exact numerical calcu- lations of van Adelsberg & Lai (2006). For weak mag-netic fields the beaming might be weaker which will fur-ther decrease the maximum pulsed fraction which we seein our models.Recently, Gusakov et al. (2020) modelled the ambipo-lar diffusion in the NS core. They noticed that magneticfield in the core could evolve significantly on timescalesof NS cooling i.e. 10 – 10 years (Igoshev & Popov 2014,2015, 2020). Such a magnetic field evolution proceedingat the lower boundary condition could noticeably changeboth the final magnetic field configuration arising as aresult of magnetic field evolution in the crust and thesurface thermal pattern. CONCLUSIONSWe perform three-dimensional magneto-thermalelectron-MHD simulations of the NS crust to studythe tangled magnetic field configurations. These config-urations were suggested to be a possible mechanism toexplain the peculiar emission properties of CCOs. Wefound that these magnetic field configurations lead to: • Formation of small hot regions with typical size of ≈ • The heating correlates with initial total magneticenergy. Most of the crustal heating due to the
CO II × erg s − .Only a fraction of this power is emitted as thermalX-ray radiation from the NS surface. • In our simulations with initial total magnetic en-ergy of E tot , = 2 . × erg, the hot regions havetemperatures 1.7 times larger than the bulk sur-face temperature of NS. In our simulations with E tot , = 10 erg, the hot regions have temper-atures which are 2.2 times larger than the bulksurface temperature. These factors are compat-ible with ones seen in the X-ray observations ofCCOs. • The maximum surface temperature decays expo-nentially on a timescale of ≈
15 kyr. The max-imum temperature approaches the bulk surfacetemperature already after ≈
20 kyr. • The resulting light curves show modulations withmaximum PFs of 9 −
11 % in the case of models 4,5, 6 and 7. These small PFs can be explained bythe fact that hot regions are located at multiplepositions on the NS, and it is impossible to choosesuch an orientation where none are present. This is a reason why larger E tot , does not necessarilylead to an increase in PF. Regions become hotter,but a few of them are still seen simultaneously. • The final poloidal magnetic field correlates withroot mean square magnetic field.Overall, the hidden magnetic field could provideenough energy to explain enhanced thermal emission ofCCOs. A part of the tangled magnetic field energy isreleased through thermal emission from NS surface viaemission of small hot spots. The small size of these spotsis related to a typical size of magnetic field loops as-sumed as the initial condition in our simulations. Thesesizes are comparable to the NS crust depth and could beproduced in a stochastic dynamo. Therefore, the maindifference between CCOs and magnetars in this scenariois the typical size of the magnetic field. In magnetars themagnetic field is large-scale, while in CCOs it is mostlysmall-scale. ACKNOWLEDGEMENTSThis work was supported by STFC grantST/S000275/1. This work was undertaken on ARC3and ARC4, part of the High Performance Computingfacilities at the University of Leeds, UK.REFERENCES