3D Magneto-thermal Simulations of Tangled Crustal Magnetic Field in Central Compact Objects

Central compact objects are young neutron stars emitting thermal X-rays with bolometric luminosities $L_X$ in the range $10^{32}$-$10^{34}$ erg/s. Gourgouliatos, Hollerbach and Igoshev recently suggested that peculiar emission properties of central compact objects can be explained by tangled magnetic field configurations formed in a stochastic dynamo during the proto-neutron star stage. In this case the magnetic field consists of multiple small-scale components with negligible contribution of global dipolar field. We study numerically three-dimensional magneto-thermal evolution of tangled crustal magnetic fields in neutron stars. We find that all configurations produce complicated surface thermal patterns which consist of multiple small hot regions located at significant separations from each other. The configurations with initial magnetic energy of $2.5-10\times 10^{47}$ erg have temperatures of hot regions that reach $\approx 0.2$ keV, to be compared with the bulk temperature of $\approx 0.1$ keV in our simulations with no cooling. A factor of two in temperature is also seen in observations of central compact objects. The hot spots produce periodic modulations in light curve with typical amplitudes of $\leq 9-11$ %. Therefore, the tangled magnetic field configuration can explain thermal emission properties of some central compact objects.


INTRODUCTION
Central compact objects (CCOs) are neutron stars (NSs) found in close proximity to the geometric center of supernova remnants (SNRs) (Pavlov et al. 2004;de Luca 2008;De Luca 2017). CCOs have soft X-ray luminosities of 10 32 −10 34 erg/s, incompatible with cooling luminosities of young NSs with ages of a few kyr. Their X-ray spectra generally contain two black-body components: the first arising from bulk NS thermal emission with temperature ≈0.1-0.2 keV, and the second produced by a small, thermally emitting area (a few percent of the NS size) with temperatures of 0.2-0.4 keV. CCO luminosities are comparable with those of some quiescent magnetars (for a review see e.g. Kaspi & Beloborodov 2017). However, CCOs show no activity similar to mag-netars (with the exception of 1E 161348-5055 in SNR RCW 103, Rea et al. 2016, which might be a magnetar with unusual formation path). There is no indication that they are binaries. CCOs are numerous: up to 30% of NSs associated with SNRs at distances up to 5 kps are CCOs (de Luca 2008), which is difficult to reconcile with the Galactic NS birthrate (Keane & Kramer 2008). It is unclear how CCOs age. CCO-like objects are not found among 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-ray emission and found none.
Long-term X-ray observations have revealed the spin period and period derivative for three CCOs (Halpern & Gotthelf 2010;Gotthelf et al. 2013). The spin periods are: RX J0822-4300 in SNR Puppis A has P = 112 ms (Gotthelf & Halpern 2009), CXOU J185238.6+004020 in Kes 79 has P = 105 ms (Gotthelf et al. 2005) and 1E 1207.4-5209 in PKS 1209-51/52 has P = 424 ms (Zavlin et al. 2000). From these measurements and their period derivatives, the magnetic dipole inferred for these CCOs can be estimated using the formula B p ≈ 3.2 × 10 19 PṖ G (Lorimer & Kramer 2004). The inferred values of B p are in the range 10 10 -10 11 G, which is one or two orders of magnitude smaller than the 10 12 G dipole fields of young radio pulsars. The typical spin-down luminosityĖ for CCOs is ∼ 10 32 erg s −1 , which means that magnetospheric currents arising from pulsar spin-down cannot be responsible for formation of compact hot regions at the CCOs' surfaces.
CCOs with measured periods and period derivatives are located among old radio pulsars in the P -Ṗ plane. While these pulsars are prominent radio sources, CCOs show no non-thermal radio emission, neither persistent nor periodic. Coherent pulsar radio emission is strongly beamed and is generated in the magnetosphere by some kind of plasma instability, or as suggested recently by non-stationary plasma discharge (Philippov et al. 2020). Even if we miss the beam of radio emission due to the pulsar orientation, the presence of a pulsar wind nebula might be expected. Young, rapidly rotating neutron stars embedded in supernova remnants are often surrounded by a compact pulsar wind nebula (Gaensler & Slane 2006). This nebula is powered by the particle wind generated by rotating, magnetised NSs. No evidence of PWN is found in relation to CCOs.
Therefore, observational properties of CCOs pose multiple questions for researchers: (1) what is the source of additional heating in these NSs, (2) why are their emitting areas so small, (3) what are the descendants of CCOs, (4) why is no radio emission detected, and (5) how different CCOs really are from magnetars. This is a part of a larger problem with variety of different NS classes. A grand unification of NSs is a scenario which tries to explain the difference between various classes of NSs based on properties and evolution of their magnetic fields (Kaspi 2010;). In the framework of this scenario, CCOs might evolve into some other category of NSs, such as radio pulsars at timescales of 10 4 -10 6 years.
Several scenarios have been proposed to answer these questions and to explain the origin and evolution of CCOs. These scenarios attribute the additional source of energy to strong hidden magnetic fields of NSs. Among these scenarios are fall-back accretion and a stochastic dynamo (Gotthelf & Halpern 2007). Magnetic fields of NSs are multi-component; they can include both poloidal and toroidal parts, as well as small-scale fields. Only the large-scale poloidal dipole magnetic field emerging from the surface can be measured using the timing technique.
In the fall-back accretion scenario, some material expelled by the supernova explosion does not have enough energy to escape from the system indefinitely, and is subsequently accreted back onto the NS (Chevalier 1989). This accretion continues after the NS crust has solidified. New infalling material buries the magnetic field, covering it with a new crustal layer (Bernal et al. 2010). The buried magnetic field gradually re-emerges through the surface over time (Ho 2011;Viganò & Pons 2012;Bernal et al. 2013). In this scenario, the excessive heating is explained by the amplification and dissipation of the magnetic field, which is compressed between the superconducting NS core and new crust formed from fall-back material. The absence of radio emission is explained by the suppression of the surface dipolar magnetic field, preventing the NS from producing enough plasma in its magnetosphere. The same explanation holds for the lack of a pulsar wind nebula. Igoshev et al. (2016) studied the evolution of small-scale magnetic field after a fall-back and found that the field (both small-and large-scale) re-emerges, so these NSs should start operating as normal radio pulsars after 10 4 -10 5 years. It is still unclear if the heat produced by the magnetic field causes a formation of small emitting regions at the surface in this scenario. If a strong fall-back occurs at the newly-born magnetar, a CCO-magnetar could be produced. This is a CCO which has an extremely strong magnetic field trapped deep in the crust which should occasionally exhibit flares.
In the stochastic dynamo scenario, the NS is born with a predominantly small-scale magnetic field, with only a weak dipole component. This small-scale field is formed during the proto-NS stage, which lasts for tens of seconds (Pons et al. 1999) and ends with the solidification of the crust. During this stage there are two regions where the matter is unstable according to the Ledoux convection criterion. One of these regions is located around densities of 10 13 g cm −3 of proto-NSs (Thompson & Murray 2001). Nagakura et al. (2020) analysed state-ofthe-art 3D simulations of supernova explosions for progenitor masses ranging from 9 M to 20 M and found that convection always develops in proto-NSs. This convection zone is an ideal place for a dynamo which could generate magnetic fields in the range 10 12 -10 15 G.
The generation of a large-scale magnetic field is thought to require rapid rotation of the proto-NS (see e.g. recent simulations by Raynaud et al. 2020). If the proto-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 ≈ 10 14 G is hidden for a distant observer, because it hardly contributes to the NS spin-down. A stochastic dynamo could also operate in convective fall-back flow (Thompson & Murray 2001), and also lead to formation of smallscale magnetic fields.
In the stochastic dynamo scenario, the excessive heating of CCOs could be explained by Ohmic decay combined with fast Hall evolution of this much stronger nondipolar magnetic field. We test it in our research using the tangled configuration of magnetic field suggested by Gourgouliatos et al. (2020). In these configurations, the magnetic fields are comparable in strength to magnetar ones, but should not produce similar stresses and therefore, should not lead to magnetar-like activity. Gourgouliatos et al. (2020); Brandenburg (2020) found that these configurations where the poloidal dipole is absent or rather small develop such a component in the course of their evolution due to an inverse Hall cascade.
The Hall cascade was first suggested for electron-MHD by Goldreich & Reisenegger (1992) with application to NS crusts. Occurrence of this cascade was confirmed in detailed numerical simulations (see e.g. Wareing & Hollerbach 2009Brandenburg 2020). Turbulence leads to redistribution of magnetic energy to both smallscale (forward cascade) and large-scale (inverse cascade) fields. In the case of tangled magnetic field configurations, an inverse cascade leads to an increase of the dipolar magnetic field with time. The small dipolar component at the beginning of NS evolution is not sufficient for the activation of the radio-pulsar mode. When the global dipolar component rises after ∼ 10 kyr of evolution, this object might start operating as a radio pulsar, which solves the problem regarding the absence of CCO descendants. Because in this case CCOs turn into normal radio pulsars with no unusual properties after a period of time.
This article is structured as follows. In Section 2 we describe the method which we use in our simulations; in Section 3 we present results of three-dimensional simulations and corresponding light curves. We discuss all results and conclude in Sections 4 and 5 respectively.

METHODS
We perform simulations in two steps. First we compute the surface temperature distribution using an updated version of the MHD code PARODY (Dormy et al. 1998;Wood & Hollerbach 2015;Gourgouliatos et al. 2016). Next we compute the light curves and pulsed fraction for different orientations of rotating NSs using the ray-tracing in general relativity.

Magneto-thermal evolution
The details of the magneto-thermal code are summarised in De Grandis et al. (2020); Igoshev et al. (2020). We solve magnetic induction and thermal diffusion equations in the solid crust of NS using the electron-MHD approximation in which only electrons carry electric charge and thermal flux. Basically, we solve two coupled partial differential equations. The first is the induction equation: where B is the magnetic field, n e is the electron number density, e is the elementary charge, c is the speed of light, σ is the electric conductivity, S e is the electron entropy, and T is the temperature. The second is the heat equation: wherek is the thermal conductivity tensor. The first two terms in our induction eq. (1) are the ones introduced by Goldreich & Reisenegger (1992), and represent the Hall evolution and Ohmic dissipation, respectively. The third term describes the so-called Biermann battery effect, i.e. magnetic field generation due to temperature gradients (Blandford et al. 1983). This effect is closely related to electron baroclinicity, which acts as a source of electron circulation, and hence magnetic induction.
In the heat eq. (2), the first term represents anisotropic thermal diffusion. While heat flows freely along magnetic field lines, heat transfer is inhibited in the direction orthogonal to field lines. The second term represents heating by Ohmic field decay, and the third term represents the entropy carried by electric currents. Unless strong temperature gradients are maintained externally, e.g. by magnetospheric heating (De Grandis et al. 2020), the last terms in eqs. (1,2) are generally very small in comparison to other terms, and play a negligible role in the field evolution.
We write the thermal conductivity tensor in the following form: where δ ij is the Kronecker delta, ijk is the Levi-Civita symbol, k B is the Boltzmann constant, and B k is the field component along the k-axis.
Electron entropy is related to temperature and chemical potential µ(r) of a degenerate, relativistic Fermi gas as follows: Similar to Gourgouliatos & Cumming (2014) we use the following radial profiles for the chemical potential, electron density and conductivity: where µ 0 = 2.9 × 10 −5 erg, n 0 = 2.5 × 10 34 cm −3 and σ 0 = 1.8 × 10 23 s −1 . For computational convenience we assume that the heat capacity is proportional to the temperature in the form (see e.g. Page et al. 2004) so that eq. (2) becomes a linear equation in T 2 .

Boundary and initial conditions
For the magnetic induction at the inner boundary, we assume the Meissner condition, i.e. zero radial component of magnetic field B r = 0 and zero tangential component of the current J t = 0 at the crust-core boundary, see De Grandis et al. (2020) and Hollerbach & Rüdiger (2004) for more details. Physically this corresponds to a scenario in which the field is expelled from the core before the crust solidifies. In reality, the core might still contain some magnetic field, although the timescale on which the core field evolves is uncertain, see e.g. Gusakov et al. (2020). For the magnetic field at the outer boundary we use vacuum boundary conditions, i.e. we assume that ∇ × B = 0 in the region outside the star. This condition is implemented numerically as: and B m t,l = 0 where B m p,l and B m t,l represent the coefficients of the poloidal and toroidal magnetic potentials, which are expanded in series of spherical harmonics of degree l and order m. Physically these boundary conditions correspond to neglecting any electric currents in the star's tenuous plasma atmosphere. We note that CCOs have no radio emission, and no pulsar wind nebula was ever detected around a CCO, so their magnetospheres might have less plasma than is typical for normal radio pulsars.
For the temperature, at the crust core interface we assume no cooling, yielding the boundary condition ∂T /∂t = 0. In reality the core cools at a rate determined by neutrino emission, which is practically insensitive to conditions within the crust. Since our focus in this study is on small-scale temperature anomalies produced within the crust we neglect the core cooling for simplicity. At the outer boundary we assume that the heat flux out of the computational domain is subsequently emitted as black-body radiation from the star's surface, i.e.
where σ S is the Stefan-Boltzmann constant, T b is the temperature at the top of the computational domain, and T s is the effective surface temperature. We assume that these two temperatures are related by the thermal blanket equation similar to a relation suggested by Gudmundsson et al. (1983), but in a more simplified form. This assumption is necessary in our calculations because the thermal relaxation timescale of the envelope is much shorter than any other evolution timescale involved in our simulations. As was shown by Gudmundsson et al. (1982), the effective surface temperature depends essentially only on the temperature in deep layers (below densities 10 10 g cm −3 ) and surface gravity. We use the same magnetic field configurations as in Gourgouliatos et al. (2020), so the initial magnetic field consists of several high-order multipoles with 10 ≤ l ≤ 20. We choose phases randomly for these fields. The amount of initial magnetic energy E tot,0 in this field is set at the beginning of simulations. A dipolar poloidal magnetic field is added to the stochastic magnetic field configuration with fixed B dip,0 value before the simulations are started.
We compute the same first five models as in Gourgouliatos et al. (2020). In models 6 and 7 we decrease the total magnetic energy in the crust by a factor of 4. We make this change because for larger energies we encounter certain numerical problems. We summarise the details of our simulations in Table 1.

Light curves and pulsed fraction
In order to compute the light curves we follow the same 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 NS, κ is the angle between magnetic dipole moment and rotational axis, and ∆Φ is the initial phase. These three angles uniquely describe NS orientation at the start of rotation. To simulate rotation of NS, we use simplified equations presented by Beloborodov (2002) combined with: where θ and φ are coordinates at the NS surface with respect to original dipolar magnetic poles, µ = cos ψ is the direction toward the position with coordinates θ, φ at infinity, and Φ ∈ [0, 2π] is the angular phase. Our eq. (13) reduces to eq. (5) of Beloborodov (2002) if a hot spot is located at the magnetic pole, i.e. θ = 0 = φ. We choose the photon beaming to be proportional to cos 2 θ (DeDeo et al. 2001) because of the following argument. In an atmosphere which is not affected by magnetic field, the darkening is described by the Hopf function, see e.g. Chandrasekhar (1950). The intensity emitted at angle 90 • is ≈ 1/3 of the intensity emitted in the normal direction. In a strong magnetic field (van Adelsberg & Lai 2006) the beaming is stronger than this, so at angles of 80 • the emission is 0.1 of its value at 10 • , and it goes to 0 at 90 • (see figure 15 of van Adelsberg & Lai 2006), which is roughly similar to a cos 2 θ behaviour. For CCOs the pulsed fractions are typically small, so it is better to overestimate the pulsed fraction in our models to check if the model has enough predictive power.
The flux emitted by each surface element of NS is integrated over the visible hemisphere (which is slightly more than half of the star due to the light bending in a strong gravitational field) to produce the visible flux for each rotational phase. We use 16 phases to simulate a light curve. For a few orientations we compute the pulsed fraction (PF) as: where F max and F min are maximum and minimum fluxes for a fixed NS rotational orientation. We select the maximum pulsed fraction among different orientations and show them in Table 1.
In Igoshev et al. (2020), we fitted observed light curves of magnetars in quiescence to estimate the NS rotational orientation. This was possible because magnetars seem to have regular large-scale magnetic fields. This is not the case for CCOs. If their magnetic fields are indeed formed as a result of a stochastic dynamo, the dipolar component is weak and changes location relatively quickly (see Gourgouliatos et al. 2020). Additional loops of magnetic field are located randomly; therefore, it is necessary to describe the location of each individual loop of magnetic field to reproduce a thermal map. Thus, there is no use in exact fitting of the light curve, since multiple slightly different magnetic field configurations will result in very similar light curves but slightly different NS rotational orientations.

RESULTS
We computed models 1-3 up to 50 Kyr, and model 4,5,6 and 7 up to 100 Kyr. We show the surface temperature distributions for models 2, 4 and 6 in Figures 1, 2  and 3. The surface temperature distributions in the case of model 1 and 3 are very similar to model 2. The surface temperatures of model 5 and 7 are very similar to ones produced in model 4 and 6, respectively, so we do not plot these maps. Basically the surface temperature distributions are defined only by the value of the initial magnetic energy, since the initial field topology is the same in all models. The initial strength of the dipolar component weakly affects the values of surface temperatures and the location of hot and cold regions.
In all cases the small-scale structure of the magnetic field leads to formation of complicated thermal pattern which includes alternating small hot and cold regions. Individual hot spots have simple shapes but temperatures of groups of spots differs in different parts of NS, making the pattern even more complicated. Overall, at 3.5 Kyr the surface temperature map is composed of a few relatively small hot regions and hot spot associations, and a much more extended colder continuum. We measure the linear sizes of several hot regions in the case of model 6 at 1.2 kyr. Two brightest hot regions have linear sizes of 2.2 and 2.6 km respectively, assuming that the NS radius is 10 km. At later stages of evolution (> 20 kyr) the temperature and size of hot regions seen   in all models decays. The hot regions become completely isolated and cool down.
The hot regions are much hotter in the case of model 4 and 6 with larger initial magnetic energy than in model 2. In the case of model 6, which has the largest magnetic field energy, the hot spots are reaching temperatures of 2.16 × 10 6 K, which is ≈ 0.19 keV. Gotthelf et al. (2013) mentioned that it is necessary for a successful CCO model to show variations of the temperature by a factor of two over the surface, because this temperature difference is seen in observations of RX J0822-4300. Our models 4, 5, 6 and 7 show temperature variations with a factor of 1.7 -2.2 over the surface.
Surface heating is caused by the release of magnetic energy stored in the crustal field, see Figure 4. During the first ≈ 20 kyr, the initial total magnetic energy E tot decays from 2.5 × 10 47 erg to ≈ 5 × 10 46 erg with typical power of ≈ 3×10 35 erg/s. In real NSs this power is emitted through the surface photon emission and neutrinos. In our simulations this power is partly emitted through the surface boundary condition and partly absorbed by the inner boundary condition.
We estimate the instantaneous exponent which describes the decay of total magnetic energy: in a way similar to work by Brandenburg (2020): We show the results of this evolution in Figure 5. The behaviour of p(t) is not monotonic. It reaches the first maximum with values ranging 0.39-0.58 (values for individual models can be found in Table 1). This level is reached at different physical times because the Hall timescale differs in these models. After this initial saturation, the p(t) value briefly declines and starts growing again reaching values of 1-3 at the end of our simulations, see Figure 5. It means that by the end of our simulations the energy of magnetic field decays exponentially due to the Ohmic decay. The values of 0.39-0.58 corresponds to the value 2/5 obtained for the case of helical magnetic fields by Brandenburg (2020). He noticed that even magnetic fields with small fraction of helicity evolve quickly to nearly 100% helical configurations. Therefore, it is not surprising that our simulations show an energy decay rate initially E = E 0 t −2/5 which is typical for helical configurations.
Another way to describe a complicated field using a single value is by using the root mean square of magnetic field B rms , similar to what Brandenburg (2020) obtained using the PENCIL code . We compute B rms over the whole volume of the NS crust. We show the evolution of the maximum temperature as a function of B rms in Figure 6. We see that similar values of B rms might correspond to different maximum temperatures depending on the evolution stage. On the other hand, the value of emerging dipolar, poloidal magnetic field seems to correlate with B rms , as seen in Figure 7. The initial value of the dipolar component gets erased quickly by the Hall evolution, and the field which emerges afterwards is related to B rms . In the case of model 2, the emergent field is 2 × 10 10 G. In the case of models 6 and 7, the dipolar field reaches values of ≈ 10 12 G, quite independently of the initial dipolar magnetic field.
It is interesting to note that the Hall evolution creates vortices similar to turbulence seen in the surface temperature maps, especially in Figure 1, age 9.5 kyr. Typically the Hall evolution creates power-spectra somewhat similar to turbulence, but with a different slope of l −2 , see e.g. Wareing & Hollerbach (2009) ;Goldreich & Reisenegger (1992). When the Hall evolution starts with initial condition consisting of a strong global dipole, it mostly preserves the dipole component and redistributes a part of its energy into small-scale magnetic fields. The intensity of these fields is small in comparison to the dipole, and they are not distinguishable as spatial structures, see e.g. Igoshev et al. (2020). In simulation 2 we suppress the initial dipole and choose the small-scale harmonics to be strong at the beginning of the simulation. As the result we see that the magnetic energy is redistributed among small-scale harmonics in such a way  that spatial vortices emerge in the thermal map, which shows that under certain conditions the Hall turbulence is very similar to regular turbulence. The system evolution of model 2 in scales larger than the crust scaleheight H is affected by the geometry of the system. This roughly corresponds to modes with ≈ πR N S /H, with being the spherical harmonic decomposition mode and R N S the radius of the star. For a value of H ∼ 1 km this corresponds to modes with > 30 where the behaviour of the turbulence will be more profound and the field will not be affected by the geometry of the system. Since the hot regions occur at multiple isolated locations at the surface, it means that a few of these regions are seen simultaneously, independently of the NS orien- tation. This effect is amplified further when the curved path of photons in the NS gravitational field is taken into account. We show the light curves for models 4 and 6 in Figure 8. Despite the factor two difference in temperature between hot and cold regions, the pulsed fraction is limited to 9-11%, similar to real CCOs.
Overall, the light curves have a few components which correspond to different hot regions. The shape of light curves is simple because the hot regions are compact. There is no use in discussing how much one light curve peak is higher/lower than another peak, because it is defined by the small-scale structure of the magnetic field. Slightly different configurations of the field could produce a very different light curve shape.

DISCUSSION
The 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 photon emission and (4) no magnetic field evolution in the core is assumed. Below we briefly explain how addition of these effects could change the results of our simulations.
NS cools down and the core temperature decays below 10 8 K assumed in our simulations within first 10 kyr of the NS life. Therefore, the surface temperatures will be smaller in extended simulations. Nevertheless, the magnetic field isolates certain regions from the core and heats them up, so the temperature difference between hot and cold regions might increase in realistic simulations with NS neutrino core cooling. Another factor which is not taken into account is the neutrino cooling of the crust. It might remove the excess heat from hot spots, limiting their temperatures by some threshold.
In our simulations we assume a temperature dependence between the deep crust and surface in form of eq. (12). In reality the temperature depends on magnetic fields as well, see e.g. Potekhin & Yakovlev (2001). This might change the surface temperature thermal map and could affect the light curves. We assume a strong beaming of the thermal emission proportional to cos 2 θ , which approximately follows the exact numerical calcu-lations of van Adelsberg & Lai (2006). For weak magnetic fields the beaming might be weaker which will further decrease the maximum pulsed fraction which we see in our models.
Recently, Gusakov et al. (2020) modelled the ambipolar diffusion in the NS core. They noticed that magnetic field in the core could evolve significantly on timescales of NS cooling i.e. 10 4 -10 6 years , 2020. Such a magnetic field evolution proceeding at the lower boundary condition could noticeably change both the final magnetic field configuration arising as a result of magnetic field evolution in the crust and the surface thermal pattern.

CONCLUSIONS
We perform three-dimensional magneto-thermal electron-MHD simulations of the NS crust to study the tangled magnetic field configurations. These configurations were suggested to be a possible mechanism to explain the peculiar emission properties of CCOs. We found that these magnetic field configurations lead to: • Formation of small hot regions with typical size of ≈ 2 km which are located at significant separations from each other.
• The heating correlates with initial total magnetic energy. Most of the crustal heating due to the currents is released during the first 10-20 kyr of NS evolution. The power released due to the magnetic field decay could reach values of 3 × 10 35 erg s −1 . Only a fraction of this power is emitted as thermal X-ray radiation from the NS surface.
• In our simulations with initial total magnetic energy of E tot,0 = 2.5×10 47 erg, the hot regions have temperatures 1.7 times larger than the bulk surface temperature of NS. In our simulations with E tot,0 = 10 48 erg, the hot regions have temperatures which are 2.2 times larger than the bulk surface temperature. These factors are compatible with ones seen in the X-ray observations of CCOs.
• The maximum surface temperature decays exponentially on a timescale of ≈ 15 kyr. The maximum temperature approaches the bulk surface temperature already after ≈ 20 kyr.
• The resulting light curves show modulations with maximum PFs of 9 − 11 % in the case of models 4, 5, 6 and 7. These small PFs can be explained by the fact that hot regions are located at multiple positions on the NS, and it is impossible to choose such an orientation where none are present. This is a reason why larger E tot,0 does not necessarily lead to an increase in PF. Regions become hotter, but a few of them are still seen simultaneously.
• The final poloidal magnetic field correlates with root mean square magnetic field.
Overall, the hidden magnetic field could provide enough energy to explain enhanced thermal emission of CCOs. A part of the tangled magnetic field energy is released through thermal emission from NS surface via emission of small hot spots. The small size of these spots is related to a typical size of magnetic field loops assumed as the initial condition in our simulations. These sizes are comparable to the NS crust depth and could be produced in a stochastic dynamo. Therefore, the main difference between CCOs and magnetars in this scenario is the typical size of the magnetic field. In magnetars the magnetic field is large-scale, while in CCOs it is mostly small-scale.