Three-dimensional Modeling of the Magnetothermal Evolution of Neutron Stars: Method and Test Cases

The NS crust comprises a Coulomb lattice in which nuclei have negligible motion. Hence, the crustal currents are produced entirely by the flow of electrons, which form a highly relativistic and strongly degenerate Fermi gas. Still, their mean velocity is typically only a tiny fraction of the speed of light. We can therefore resort to the (nonrelativistic) electron magnetohydrodynamics (eMHD) approximation in treating the crustal dynamics. The evolution of the magnetic field B in the crust is described by the induction equation that, also taking into account the effects of thermal coupling, can be written in the form

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}=-c\,{\boldsymbol{\nabla }}\times \left[{{\boldsymbol{\sigma }}}^{-1}\cdot {\boldsymbol{J}}+{\boldsymbol{G}}\cdot {\boldsymbol{\nabla }}T-{\boldsymbol{\nabla }}\mu /{\rm{e}}\right]\end{eqnarray} \tag{ 1 }$$

where the term in square brackets is the electric field ${\boldsymbol{E}}$ as given by the generalized Ohm's law. Here, c is the speed of light, ${\rm{e}}$ the electron charge, ${\boldsymbol{\sigma }}$ and ${\boldsymbol{G}}$ are the electric conductivity and the thermopower tensors, and μ is the electron chemical potential. The latter, for a degenerate relativistic Fermi gas, depends only on the density, $\mu =c\hslash {(3{\pi }^{2}n)}^{1/3}$, where ℏ is the reduced Planck constant. Ignoring the displacement current, the electron current is given by ${\boldsymbol{J}}=c\,{\boldsymbol{\nabla }}\times {\boldsymbol{B}}/4\pi$.

Assuming the temperature of the crust to be well below the electron degeneracy temperature but above the ion plasma temperature, scattering of electrons can be described in terms of an energy-dependent relaxation time τ (e.g., Ziman 1972; Urpin & Yakovlev 1980), and the electron conductivity can then be approximated as

$$\begin{eqnarray}&&{\left({{\boldsymbol{\sigma }}}^{-1}\right)}_{{ij}}={\sigma }^{-1}{\delta }_{{ij}}+\displaystyle \frac{{\epsilon }_{{ijk}}{B}_{k}}{{\rm{e}}{cn}}\end{eqnarray} \tag{ 2 }$$

where the symmetric part is

$$\begin{eqnarray}&&\sigma ={{\rm{e}}}^{2}{c}^{2}\,\displaystyle \frac{n\tau (\mu )}{\mu }\end{eqnarray} \tag{ 3 }$$

and the antisymmetric part represents the Hall effect.

The thermopower can be calculated using the Mott formula and in general has an isotropic part and a part proportional to the conductivity tensor. In our model, we include only the isotropic part (the so-called Seebeck term), which is responsible for the Biermann battery effect,

$$\begin{eqnarray}&&{\boldsymbol{G}}={\rm{e}}{\left.\displaystyle \frac{\partial {{\boldsymbol{\sigma }}}^{-1}}{\partial \mu }\right|}_{T}\cdot {\boldsymbol{k}}\simeq -\displaystyle \frac{{\pi }^{2}{k}_{B}^{2}T}{{\rm{e}}\mu }{\delta }_{{ij}}\end{eqnarray} \tag{ 4 }$$

where the approximate equality is obtained by further assuming electrons to form a perfect Fermi gas. Here, ${\boldsymbol{k}}$ is the thermal conductivity tensor, which is taken to be proportional to ${\boldsymbol{\sigma }}$, according to the Wiedemann–Franz law:

$$\begin{eqnarray}&&{\boldsymbol{k}}=\displaystyle \frac{{\pi }^{2}{k}_{{\rm{B}}}^{2}T}{3{{\rm{e}}}^{2}}\,{\boldsymbol{\sigma }}\end{eqnarray} \tag{ 5 }$$

where ${k}_{{\rm{B}}}$ is Boltzmann's constant.

The evolution of temperature is, in turn, governed by the heat equation,

$$\begin{eqnarray}&&{C}_{V}\displaystyle \frac{\partial T}{\partial t}=-{\boldsymbol{\nabla }}\cdot \left(T{\boldsymbol{G}}\cdot {\boldsymbol{J}}-{\boldsymbol{k}}\cdot {\boldsymbol{\nabla }}T-\displaystyle \frac{\mu }{{\rm{e}}}{\boldsymbol{J}}\right)+{\boldsymbol{E}}\cdot {\boldsymbol{J}}+{N}_{\nu }\end{eqnarray} \tag{ 6 }$$

where C_V is the heat capacity (per unit volume) of the crust, N_ν is the neutrino emissivity due to weak processes, and the term in parentheses is the electron energy flux.

Although general-relativistic effects can be accounted for with no inherent difficulty in Equations (1) and (6) (see, e.g., Pons et al. 2009 and Viganò et al. 2013), they are not included here. The reason for this is twofold. First, given the small thickness of the crust, they are of limited importance and will not change our results qualitatively and, second, a proper general-relativistic treatment impacts the boundary conditions that are imposed on the evolution equations (see Section 2.2). While this poses no serious problem in 2D, it becomes quite troublesome in 3D. Equation (1) is analogous to the one solved in Gourgouliatos & Cumming (2014a; with the addition of thermocoupling terms) and Viganò et al. (2013), where (a different version of) Equation (6) was included as well.

In order to simplify the equations, all physical quantities are henceforth expressed in terms of values typical of the outer crust in a magnetar. In particular, the temperature, magnetic field, relaxation time, and chemical potential are normalized to ${T}_{0}={10}^{8}\,{\rm{K}}$, ${B}_{0}={10}^{4}{\rm{G}}$, ${\mu }_{0}=2.9\times {10}^{-5}\,\mathrm{erg}$, and ${\tau }_{0}=9.9\,\times {10}^{-19}\,{\rm{s}}$. This implies that the reference values for the number density n and the conductivity $\eta ={c}^{2}/(4\pi \sigma )$ are ${n}_{0}\,\simeq 2.6\times {10}^{34}$ cm⁻³ and ${\eta }_{0}\simeq 3.9\times {10}^{-4}$ cm²s⁻¹. It is also useful to introduce four length scales, which are relevant to the electron dynamics,

$$\begin{eqnarray}&&\lambda ={\left(\displaystyle \frac{{k}_{{\rm{B}}}{T}_{0}}{4\pi {n}_{0}{{\rm{e}}}^{2}}\right)}^{1/2}\quad \mathrm{Debye}\ \mathrm{length}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&d={\left(\displaystyle \frac{{\mu }_{0}}{4\pi {n}_{0}{{\rm{e}}}^{2}}\right)}^{1/2}\quad \mathrm{skin}-\mathrm{depth}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&L=\displaystyle \frac{{\mu }_{0}}{{\rm{e}}{B}_{0}}\quad \mathrm{Larmor}\ \mathrm{radius}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&l=c{\tau }_{0}\quad \mathrm{mean}\ \mathrm{free}\ \mathrm{path},\end{eqnarray} \tag{ 10 }$$

supplemented by the star radius ${R}_{\star }=10\,\mathrm{km}$ as the macroscopic length scale. Furthermore, the ohmic time ${\tau }_{O}={R}_{\star }^{2}/{\eta }_{0}\,\simeq 8\times {10}^{7}\,\mathrm{yr}$ is taken as the reference timescale and ${C}_{V0}={k}_{{\rm{B}}}{n}_{0}$ as the scale for the heat capacity.

The evolution equations to be solved for B and T then become

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}=\mathrm{Se}{\boldsymbol{\nabla }}\left(\displaystyle \frac{1}{\mu }\right)\times {\boldsymbol{\nabla }}{T}^{2}+\mathrm{Ha}\,{\boldsymbol{\nabla }}\times \left[\displaystyle \frac{1}{{\mu }^{3}}{\boldsymbol{B}}\times ({\boldsymbol{\nabla }}\times {\boldsymbol{B}})\right]\\ \ \ +\,-{\boldsymbol{\nabla }}\times \left[\displaystyle \frac{1}{\tau {\mu }^{2}}{\boldsymbol{\nabla }}\times {\boldsymbol{B}}\right]\end{array}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{1}{\mathrm{Ro}}\displaystyle \frac{{C}_{V}}{T}\displaystyle \frac{\partial {T}^{2}}{\partial t} & = & {\boldsymbol{\nabla }}\cdot \left(\tau {\mu }^{2}{\boldsymbol{\chi }}\cdot {\boldsymbol{\nabla }}{T}^{2}\right)+\displaystyle \frac{\mathrm{Pe}}{\mathrm{Se}}\displaystyle \frac{| {\boldsymbol{\nabla }}\times {\boldsymbol{B}}{| }^{2}}{\tau {\mu }^{2}}\\ & & +\,\mathrm{Pe}\,\mu ({\boldsymbol{\nabla }}\times {\boldsymbol{B}})\cdot {\boldsymbol{\nabla }}\left(\displaystyle \frac{{T}^{2}}{{\mu }^{2}}\right)+\displaystyle \frac{1}{\mathrm{Ro}}{N}_{\nu }\end{array}\end{eqnarray} \tag{ 12 }$$

where we defined

$$\begin{eqnarray}&&{\chi }_{{ij}}=\displaystyle \frac{{\delta }_{{ij}}+{\mathrm{Ha}}^{2}{\left(\tau /\mu \right)}^{2}{B}_{i}{B}_{j}-\mathrm{Ha}(\tau /\mu ){\epsilon }_{{ijk}}{B}_{k}}{1+{\mathrm{Ha}}^{2}{\left(\tau /\mu \right)}^{2}| {\boldsymbol{B}}{| }^{2}}\end{eqnarray} \tag{ 13 }$$

and introduced the adimensional numbers

$$\begin{eqnarray}&&\mathrm{Ha}=l/L\simeq 50\ \mathrm{Hall}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\mathrm{Se}=\displaystyle \frac{{\pi }^{2}{Ll}{\lambda }^{4}}{2{d}^{6}}\simeq 0.05\quad \mathrm{Seebeck}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\mathrm{Pe}=\displaystyle \frac{3{d}^{2}}{{Ll}}\simeq 6\times {10}^{-5}\quad \mathrm{Peclet}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{\mathrm{Ro}}=\sqrt{\displaystyle \frac{3}{2{\pi }^{2}}\displaystyle \frac{\mathrm{Pe}}{{\mathrm{SeHa}}^{2}}}\simeq 3\times {10}^{-4}\quad \mathrm{Roberts}.\end{eqnarray} \tag{ 17 }$$

Performing adimensionalization with these scales, the neutrino emissivity is expressed in units of ${N}_{\nu }^{0}=8\pi {{\rm{e}}}^{2}{c}^{2}{R}^{2}{\tau }_{0}/{\mu }_{0}{k}_{b}{T}_{0}\,=1.3\times {10}^{14}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-3}$.

The large thermal conductivity of the crust, which is reflected in the fact that $\mathrm{Ro}\gg 1$, means that on the timescale of magnetic evolution, the term dependent on the heat capacity of the crust is subdominant. Rather than using a detailed microphysical model for C_V, we therefore simply take ${C}_{V}=\tau {\mu }^{2}T$, which implies a constant effective thermal diffusivity throughout the crust. Under this assumption, Equation (12) depends on temperature only through T², which proves to be an advantageous feature for numerical implementation.

Under the eMHD approximation, electrons and protons in the crust have equal, time-independent number densities n. We will assume that the crust is spherically symmetric, so that n is a function of the radius r alone. With our definition of the chemical potential (see above), this implies that μ depends on r, and we take

$$\begin{eqnarray}&&\mu (r)={\left(1+\displaystyle \frac{1-r}{0.0463}\right)}^{4/3},\end{eqnarray} \tag{ 18 }$$

following Gourgouliatos & Cumming (2014a); $\mu (r)$ increases from unity at the outer boundary (r = 1) to $\simeq 4.6$ at the inner boundary (r = 0.9). Moreover, we also assume that τ is a function of r only and, in particular, we take $\tau \equiv 1$.

The density profile corresponds to the crust model with the impurity parameter $Q\simeq 3$ of Cumming et al. (2004). The assumption of taking the relaxation time τ to be independent of temperature is adequate in the lower crust, while it is just an approximation in the upper crust, where scattering is dominated by phonons (Potekhin et al. 2015). We, nevertheless, note that taking $\tau \equiv 1$, the conductivity in the upper crust corresponds to the phonon conductivity at a realistic temperature, $T\approx {10}^{8}\,{\rm{K}}$. We assume an Fe–Ni crust, without accounting for chemical composition stratification.

The emission of neutrinos in the crust is due to a large variety of reactions. In this work, the four dominant contributions are taken into account, namely, phonon decay, neutrino pair production, neutrino bremsstrahlung, and neutrino synchrotron emission,

$$\begin{eqnarray}\begin{array}{rcl}{N}_{\nu }(n,T,{\boldsymbol{B}}) & = & {N}_{\mathrm{ph}}(n,T)+{N}_{\mathrm{pair}}(n,T)\\ & & +\,{N}_{\mathrm{bre}}(n,T)+{N}_{\mathrm{syn}}(n,T,{\boldsymbol{B}}).\end{array}\end{eqnarray} \tag{ 19 }$$

We reference Ofengeim et al. (2014) for neutrino pair bremsstrahlung decay and Kantor & Gusakov (2007) for phonon decay. A complete review can be found in Yakovlev et al. (2001). These papers provide fitting formulae for numerical evaluation, that were implemented in our code.

Solution of the evolution Equations (1) and (6) requires boundary conditions that reflect a number of physical prescriptions at the core–crust interface and at the star surface. We assume that all magnetic flux has been expelled from the superconducting core. This requires that the normal component of the magnetic field and the tangential component of the electric field must vanish at the core–crust boundary, $r={r}_{c}$. The latter results in a nonlinear boundary condition for the magnetic field due to the presence of the Hall term. Nevertheless, this contribution is negligible near the bottom of the crust due to the high electron density (see Hollerbach & Rüdiger 2004) and can hence be discarded. This allows the boundary conditions to be written in terms of the radial magnetic field and of the tangential component of the current J_t,

$$\begin{eqnarray}&&{B}_{r}({r}_{c})=0\qquad {J}_{t}({r}_{c})=0.\end{eqnarray} \tag{ 20 }$$

We assume that the electrical conductivity in the magnetosphere is negligible in comparison with that of the crust, and therefore, we match the field at the crust outer boundary to a potential one. This can be achieved in a very natural way by exploiting the spectral nature of our code, which decomposes the field using the spherical harmonics ${Y}_{{\ell }}^{m}(\theta ,\phi )$ as the basis (see Section 2.3), after introducing a poloidal–toroidal decomposition:

$$\begin{eqnarray}&&{\boldsymbol{B}}={\boldsymbol{\nabla }}\times {\boldsymbol{\nabla }}\times ({\boldsymbol{r}}{B}_{\mathrm{pol}})+{\boldsymbol{\nabla }}\times ({\boldsymbol{r}}{B}_{\mathrm{tor}}).\end{eqnarray} \tag{ 21 }$$

In such a representation, each mode of a potential field is purely poloidal, such that $({B}_{\mathrm{pol}}{)}_{{\ell }}^{m}\propto {r}^{-({\ell }+1)}$. Hence, the boundary condition can be met by requiring that

$$\begin{eqnarray}&&\begin{array}{l}{\left.\displaystyle \frac{\partial }{\partial r}{\left({B}_{\mathrm{pol}}\right)}_{{\ell }}^{m}+\displaystyle \frac{{\ell }+1}{r}{\left({B}_{\mathrm{pol}}\right)}_{{\ell }}^{m}\right|}_{r={R}_{\star }}=0;\\ \quad {\left.{\left({B}_{\mathrm{tor}}\right)}_{{\ell }}^{m}\right|}_{r={R}_{\star }}=0.\end{array}\end{eqnarray} \tag{ 22 }$$

The core conducts heat even more efficiently than the crust and is therefore approximately isothermal. It cools by neutrino emission according to the equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial {T}_{c}}{\partial t}=-\displaystyle \frac{{N}_{c}(T)}{{C}_{c}}\end{eqnarray} \tag{ 23 }$$

where C_c is the core specific heat and ${N}_{c}(T)={N}_{0}{T}^{k}$ the neutrino emissivity of the core. The star long-term thermal evolution is governed by Equation (23) once the neutrino emissivity is specified. Our model uses a standard slow-cooling scenario with k = 8, ${C}_{c}={10}^{20}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{{\rm{K}}}^{-2}$, and ${N}_{0}={10}^{21}\,\mathrm{erg}\,{{\rm{c}}}^{-3}\,{\rm{m}}\,{{\rm{K}}}^{-8}$ (Page et al. 2004).

The surface temperature is controlled by the properties of the thermal blanketing envelope. This layer is geometrically very thin but hosts a large temperature gradient. Thus, the widespread approach is to treat it separately, using a plane-parallel approximation to obtain a relation between the temperature at the bottom of the envelope T_b (that is, the temperature of the top of the crust) and the surface temperature T_s (Gudmundsson et al. 1983). Assuming that no energy gains or losses occur in the envelope, the temperature gradient at the top of the crust is given by (Tsuruta & Cameron 1966)

$$\begin{eqnarray}&&-({\boldsymbol{k}}\cdot {\boldsymbol{\nabla }}T)\cdot \hat{{\boldsymbol{r}}}={\sigma }_{{\rm\small{}}\mathrm{sb}}\,{T}_{s}^{4}({T}_{b},{\boldsymbol{B}}),\end{eqnarray} \tag{ 24 }$$

where the left-hand side is evaluated at the top of the crust and $\hat{{\boldsymbol{r}}}$ is the radial unit vector. We have chosen the form

$$\begin{eqnarray}&&{T}_{s}({T}_{b},g,{\boldsymbol{B}})={T}_{s}^{(0)}({T}_{b},g){ \mathcal X }({T}_{b},{\boldsymbol{B}})\end{eqnarray} \tag{ 25 }$$

where g is the gravitational acceleration at the surface. We used the expressions for ${T}_{s}^{(0)}$ as calculated in Gudmundsson et al. (1983) for an iron envelope ignoring magnetic fields and the magnetic correction ${ \mathcal X }({T}_{b},{\boldsymbol{B}})$ obtained in Potekhin & Yakovlev (2001).

Equations (11) and (12) were solved in three dimensions using a suitably modified version of the code parody, which was originally developed by Dormy et al. (1998) and Aubert et al. (2008). A version of the same code, which did not include the thermomagnetic coupling, was first used to investigate the magnetic field evolution in NSs in Wood & Hollerbach (2015). The code is pseudo-spectral: it uses a finite grid in the radial direction and an expansion in spherical harmonics ${Y}_{{\ell }}^{m}(\theta ,\phi )$ for the angular part. The NS crust is assumed to be a perfect spherical shell. The time-stepping algorithm is Crank–Nicholson for the ohmic diffusion, backward Euler for the isotropic part of the thermal diffusion, and Adams–Bashforth for all other terms.

We typically use 128 radial grid points and spherical harmonics up to degree ${\ell }\approx 100$, obtaining a typical resolution of ≲100 m on the surface. Parallelization is implemented using mpi, and the code is run on a cluster of CPUs. Work is distributed in such a way that each thread takes care of a spherical shell containing ${N}_{r}\mathrm{mod}{N}_{\mathrm{cores}}$ points, where N_r is the radial grid size. In order to compute space derivatives within our finite difference scheme in each thread, a single shell should contain at least four grid points; hence, to achieve the desired resolution, the code is typically run on $\lesssim 32$ cores.

In order to set the initial⁷ magnetic configuration for our simulations, we followed the widespread approach of confining the field in simple, large-scale structures (Rüdiger et al. 2013). In particular, we selected a force-free B-field matching our boundary conditions, with both nonzero poloidal (ℓ = 1, m = 0) and toroidal (ℓ = 2, m = 0) components. For such a field, the components of Equation (21) take the form ${B}_{{\ell }}^{m}\propto {Y}_{{\ell }}^{m}(\theta ,\phi ){\zeta }_{{\ell }}(r)/r$, where ${\zeta }_{{\ell }}$ is a linear combination of spherical Bessel functions of degree $\pm {\ell }$, constructed in such a way as to obey the boundary conditions (see Chandrasekhar & Kendall 1957 for a full derivation). We stress that the evolution of the poloidal/toroidal components is strongly coupled by the action of the Hall term, which can transfer energy both ways between them (Pons & Geppert 2007). The initial temperature profile is assumed to be a constant $T(r,t=0)\equiv {10}^{8}\,{\rm{K}}$, but we note that the overall evolution is virtually independent of this choice. This is due to the fact that the term $\propto \partial T/\partial t$ in Equation (12) is suppressed by a factor ${\mathrm{Ro}}^{-1}\approx {10}^{-4}$, so that the temperature rapidly achieves a quasi-steady state.

The evolution of the B-field over a few Hall timescales is shown in Figure 1 for three different initial magnetic configurations: a purely dipolar field (${B}_{\mathrm{pol}}(0)\approx {10}^{14}\,{\rm{G}}$, ${B}_{\mathrm{tor}}(0)=0$) and a field with poloidal and toroidal components of the same order but opposite relative polarities (${B}_{\mathrm{pol}}(0)\,\approx \pm {B}_{\mathrm{tor}}(0)\approx {10}^{14}\,{\rm{G}}$). Our simulations confirm the previous finding that the magnetic field evolves toward the so-called Hall attractor (Gourgouliatos & Cumming 2014b), where the magnetic field tends to reach a configuration dominated by the modes ℓ = 1, 2, 3, 5, 7 (see again Figure 1). The dominance of odd modes with respect to the nearby even ones is a general feature of the Hall attractor. We remark that for a better comparison with previous works (Viganò et al. 2013; Turolla et al. 2011), we chose initial conditions that are essentially axisymmetric. Our results show that initially axisymmetric configurations tend to maintain their symmetry as they evolve.

Zoom In Zoom Out Reset image size

The components of the magnetic field in spherical coordinates, B_r, B_θ, and B_ϕ, at the beginning of the simulation (t = 0) and at $t=3\times {10}^{4}\,\mathrm{yr}\simeq {\tau }_{H}$ are shown in Figure 2 for the case with ${B}_{\mathrm{pol}}(0)\approx +{B}_{\mathrm{tor}}(0)$. A general feature of the magnetic field is the appearance of an equatorial structure in which the field is stronger (Gourgouliatos et al. 2016) and of small-scale structures due to the Hall term.

Zoom In Zoom Out Reset image size

As already mentioned, the temperature distribution tends to follow the magnetic field. The structure of the Hall attractor, in which an equatorial current ring forms, is reflected in a hotter equatorial region. Moreover, formula (5) implies that heat tends to be transported preferentially along the field lines. Hence, the equatorial region is hotter not only because of higher dissipation, but also because heat is trapped by the closed field lines appearing in that region. Figure 3 shows a typical case, which is representative—at least qualitatively—of all our nearly axisymmetric runs. We note that, owing to the dependence of the properties of the heat-blanketing envelope on the geometry of the magnetic field, the observable surface map can be quite different from the one on the top of the crust. As an example, the last panel of Figure 3 shows the surface temperature for the very same case: the overall topology is quite different, as the equatorial belt is not just hotter, but instead shows a colder ring at the very equator (see e.g., the recent results in Kondratyev et al. 2019, in which a similar behavior is discussed in a 3D stationary framework). Even though the various features are on a large scale, they exhibit a smaller scale—yet well-resolved—structure, due to the Hall term.

Zoom In Zoom Out Reset image size

As already noted in Gourgouliatos & Pons (2020), the presence of a strong toroidal field can trigger a resistive tearing eMHD instability (Wood et al. 2014). This instability, even when starting from an initial condition that is essentially symmetric, produces non-axisymmetric small-scale magnetic structures that, due to Joule dissipation, translate into localized heat deposition. A strong toroidal component in the star crust passes undetected and is invoked to explain the observed activity in the so-called low-B magnetars, i.e., sources with a dipole field comparable to that of the radio-pulsar population (see Section 5.1 for further details).

To explore this issue better, we ran a simulation assuming an ℓ = 1, m = 0 initial magnetic field with a poloidal field ${B}_{\mathrm{pol}}(0)\approx {10}^{14}\,{\rm{G}}$ and a toroidal one ${B}_{\mathrm{tor}}\approx 4\times {10}^{15}\,{\rm{G}}$. Given the nature of the solution we are looking for, the resolution for this case was improved to ${{\ell }}_{\max }=250$, corresponding to cells of a few tenths of meters on the surface. Indeed, an instability is triggered after about a Hall time ${\tau }_{H}$. The spectrum of all the ℓ modes at $t\simeq {10}^{4}\,\mathrm{yr}$ (Figure 4) exhibits the characteristic features of the Hall attractor: even modes are suppressed with respect to the nearby odd ones up to ${\ell }\lesssim 100$, and this produces the typical wavy profile. However, the onset of an instability is marked by the appearance of well-resolved structures that form up to ${\ell }\simeq 100$, with a complex structure of secondary peaks on top of the Hall structure. The flatness of the spectrum at high ℓ guarantees that the instability is of physical and not numerical origin. The slight increase at very high ℓ is due to numerical aliasing. The spectrum of m modes, on the other hand, is sharply peaked toward zero and hence is not shown.

Zoom In Zoom Out Reset image size

In the small structures, the magnetic field can reach values up to $\sim 2\times {10}^{16}\,{\rm{G}}$, and this drives a local temperature increase, as shown in Figure 5. Such strong fields generate high magnetic stresses in the crust. As a gauge to determine whether such stresses are strong enough to lead to crustal failure, we compared them to the maximum mechanical yield of the crust through the von Mises criterion (see, e.g., Pons & Rea 2012; Lander & Gourgouliatos 2019),

$$\begin{eqnarray}&&\sqrt{{\bar{M}}_{{ij}}{\bar{M}}^{{ij}}}\gtrsim {\tau }_{\max }(n,T)\end{eqnarray} \tag{ 28 }$$

where ${\bar{M}}_{{ij}}$ is the traceless part of the magnetic stress tensor ${M}_{{ij}}={B}_{i}{B}_{j}/4\pi$. Chugunov & Horowitz (2010) derived estimates for ${\tau }_{\max }$ by means of molecular dynamic simulations and elucidated the strong dependence of the breaking stress on temperature. In our calculation, we used the fit they provide for the maximum crustal yield,

$$\begin{eqnarray}&&{\tau }_{\max }=\left(0.0195-\displaystyle \frac{1.27}{{\rm{\Gamma }}-71}\right){n}_{i}\displaystyle \frac{{Z}^{2}{{\rm{e}}}^{2}}{a}\end{eqnarray} \tag{ 29 }$$

where ${\rm{\Gamma }}={Z}^{2}{{\rm{e}}}^{2}/{{ak}}_{B}T$ is the classical Coulomb coupling parameter, n_i is the ion density, $a={(4\pi {n}_{i}/3)}^{-1/3}$ is the ion sphere radius and, following Horowitz et al. (2007), we took Z = 29.4 as the mean ion charge in the Fe–Ni crust. Because ${\tau }_{\max }$ decreases at higher T, a 3D, coupled magnetothermal code provides the most accurate way to investigate the onset of crustal failures in magnetars.

Zoom In Zoom Out Reset image size

Figure 6 shows the ratio between the magnetic and the breaking stress in our simulation after a time $6.2\,\times {10}^{3}\,\mathrm{yr}\lesssim {\tau }_{H}$. The map refers to the region of the crust where the ratio is maximum, at about half the crust depth. The magnetic stress reaches values up to $\sim 50 \%$ of the maximum yield in our simulation so that the von Mises criterion for crustal yielding is likely to be fulfilled. Crustal failures can therefore be triggered in our simulation because of the large magnetic stress coupled with the heating produced by magnetic dissipation, which significantly increases the temperature in the equatorial region, thus lowering the breaking stress. The instability lasts for some 1000 yr before it is damped by dissipation. This directly concerns magnetar activity, because magnetically induced crustal failures ("starquakes") are thought to be responsible for magnetar bursts and outbursts (Pons & Rea 2012).

Zoom In Zoom Out Reset image size

We conclude this section by noticing that the use of the von Mises criterion as expressed by Equation (28), albeit widespread in the literature, should be taken with some care. In fact, it does not take into account the effects of the enormous gravity, which tends to inhibit any radial displacement (see, e.g., Haskell 2008). However, because the resistive tearing instability arises as a consequence of the presence of a strong toroidal field, our result is not much affected even when setting to zero all radial shear terms (the maximum stress-to-yield ratio decreases from ∼50% to 45%). Nevertheless, only a consistent, nonlocal calculation that takes into account the global hydrostatic structure of the crust could unambiguously solve the issue.

As already mentioned, magnetar activity is believed to be associated with crustal failures (see Section 3.2). However, the crustal dynamics in such events is little explored as yet, owing to its inherent complexity (e.g., the crust may flow plastically; Lander 2016). Such a study is beyond the capability of our code, which does not incorporate a description of the motion of crustal matter.

As a minimal model to address the physics of crustal failures, and in particular the way in which heat is transported to the surface, we therefore performed a simulation in which energy is injected during a short time interval in a localized region of the crust, much in the same way as in Pons & Rea (2012), but exploiting our fully 3D approach. As the background state, we take an NS with an initial field ${B}_{\mathrm{pol}}\approx {10}^{12}\,{\rm{G}}$ and ${B}_{\mathrm{tor}}\approx {10}^{13}\,{\rm{G}}$ that has been consistently evolved for a Hall time. This high toroidal field configuration was chosen in the spirit of the results of Section 3.2 and mimics a low-B magnetar.

In our test model, heat has been released in the northern hemisphere and in the innermost half of the crust, assuming a Gaussian profile along the three spatial dimensions with ${\sigma }_{r}\simeq 100\,{\rm{m}}$, ${\sigma }_{\theta }\simeq {\sigma }_{\phi }\simeq \pi /5\,\mathrm{rad}$. The additional heating term in Equation (12) is $\dot{H}\simeq 5\times {10}^{37}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, modulated by the Gaussian profile. Heating is assumed to be quasi-instantaneous (the $\dot{H}$ term is active for ${\rm{\Delta }}{t}_{\mathrm{inj}}\simeq 3\,{\rm{s}}$). The NS luminosity has then been calculated assuming blackbody emission at the local temperature, after deriving T_s from Equation (25). The time evolution follows a typical FRED (fast-rise-exponential-decay) pattern, as shown in Figure 7. The two curves in Figure 7 illustrate the role played by neutrino losses in the crust. The temperature increase produced by heat deposition, in fact, is large enough in this case to make neutrino emission sizable (contrary to what occurs when the crust is not heated), and this results in a photon luminosity lower by a factor of ∼2 with respect to the case in which neutrino losses are turned off. In the present case, the peak luminosity is $\sim {10}^{33}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, with an increase of a factor ≈10 above the quiescent level.

Zoom In Zoom Out Reset image size

The hot structure that develops onto the surface exhibits a somehow peculiar evolution. In fact, its shape is determined by heat diffusion, which is not isotropic but depends on the magnetic field direction according to Equations (2) and (5). Hence, heat tends to flow along field lines. Figure 8 shows how during the luminosity rise time, heat is not just flowing radially to reach the surface, but does so following the magnetic field. Moreover, once it is formed, the hotter region tends to drift as it cools down, both in latitude, toward the equator, and in longitude. Such behavior is clearly visible in Figure 9, which shows four snapshots of the heated surface patch evolution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The duration of this event is of some thousands of years, with a rise time of about a century (however, timescales in this case are affected by code limitations; see Section 5.2). Still, on a qualitative basis, it can be taken as a representation of the observed flux variation during magnetar outbursts, which happen on shorter timescales (Coti Zelati et al. 2018).

The framework discussed in Section 4.1 can be easily adapted to study the shape of pulsar hot spots. In fact, the physical ingredients remain the same as long as it can be assumed that no other effects apart from heating come into play. The major difference with respect to the model presented in the previous section is that now energy is deposited in the outermost crustal layers, as it is the case, e.g., for the heat deposited by backflowing currents on the surface of radio pulsars.

As a background state, we take an NS with initial field ${B}_{{pol}}\approx {B}_{{tor}}\simeq {10}^{12}\,{\rm{G}}$ and temperature $T\simeq 5\times {10}^{7}\,{\rm{K}}$, which was evolved for some Hall times, $t\sim {10}^{5}\,\mathrm{yr}$. Then, a heating source is placed in a small region of size ∼0.5 km close to the magnetic pole. We chose this spot position because even though the external magnetic field in our simulations is not a pure dipole, its qualitative shape is similar to it, as displayed in Figure 10, and heating of the polar regions is to be expected.⁸

Zoom In Zoom Out Reset image size

Backflowing currents in pulsars can reach a depth ranging from about a tenth to the entire width of the crust (Karageorgopoulos et al. 2019). In our simulations, we choose to insert the heat source uniformly from the surface down to a quarter of the crust width. We at any rate checked that different depth values provide quite similar results, possibly because this length is anyway much smaller than the other relevant lengths in the problem.

Results show that if heat injection is steady, the hot spot reaches a state of quasi-equilibrium in a few years. Starting with a heated patch of size ∼1 km, the spot tends to assume a quasi-circular shape, staying on top of the injection region. The evolution is shown in Figure 11 (top row), where the initial shape and the equilibrium configuration of the hot spot are compared for steady heat injection $\dot{H}=5\times {10}^{25}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. We then followed the evolution of the same spot after the heating term is turned off. In a time ≈1000 yr, the spot cools down in such a way that only a ring corresponding to the region rim is left. In the final cooling phases (Figure 11, bottom row), a crescent-like structure drifting toward the equator becomes visible. It is a somehow subdominant feature, as the main hot spot is still in correspondence with the initially heated region, but we observed that it can eventually become hotter than the central spot in the very late cooling stages (when temperature differences from the background are small). We note that some asymmetry in the position/shape of the initially heated region with respect to the magnetic pole is necessary for the formation of crescent-like structures during the evolution. Heat injection in a circular patch exactly below the magnetic pole results in a nearly circular cooling spot. However, simulations show that such deviations need to be indeed small, and in real sources, they are expected to be produced, e.g., by the effect of the coupling of crustal heating currents with the rotation of the star (Karageorgopoulos et al. 2019) or the presence of subdominant nondipolar components.

Zoom In Zoom Out Reset image size

In the case of a higher heat injection, $\dot{H}=5\times {10}^{26}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, we observe a similar phenomenology, but a turbulent-like pattern emerges (see Figure 12), which is nonetheless well resolved by our grid. This results in a more complex evolution of the shape as the spot cools down. Its relic, in fact, gets fragmented into many smaller structures that do not exhibit the ordered, ring-like shape of the previous case. A drifting crescent-shaped subdominant structure is again formed in the final phases.

Zoom In Zoom Out Reset image size

Moreover, the backreaction on the magnetic field becomes important: in fact, in this situation, a temperature gradient perpendicular to the (radial) density one develops, hence the Biermann battery effect, which provides negligible feedback in the long-term evolution of isolated NSs, can give rise to a substantial local enhancement of the magnetic field. Such behavior is displayed in Figure 13. When the stationary state is reached, some small magnetic structures appears on top of the ≈10¹² G large-scale (quasi) dipolar field, where the field strength can reach values up to $6\times {10}^{14}\,{\rm{G}}$. Thus, localized heating may also account for small-scale magnetic structures in the crust, which are therefore not originated by dynamo-like processes.

Zoom In Zoom Out Reset image size

The appearance of these crustal magnetic features is reflected in the creation of local magnetic structures in the magnetosphere, even though the overall B-field remains very close to dipolar. Figure 14 shows the external field lines for the same case as in Figure 13, as derived by solving the magnetospheric structure (see footnote 2). After subtracting the contribution of the m = 0 modes, which are dominated by the dipolar field, a small magnetic field loop is clearly visible above the heated region, extending outwards to a distance $\lesssim {R}_{\star }$ with a typical strength $\approx {10}^{9}\,{\rm{G}}$. This shows that magnetic structures are not necessarily confined to the crust but can extend in the inner magnetosphere.

Zoom In Zoom Out Reset image size

The cooling phase lasts some thousands of years. It is therefore possible that the aftermath of powerful heating events can produce a long-lasting thermal structure on an NS crust, evolving in complex patterns along field lines.

Comparison with low-B magnetars—The presence of strong toroidal fields in magnetars has long been invoked to explain their distinguishing activity compared to radio pulsars with similar spin-down magnetic fields, the high-B pulsars with ${B}_{\mathrm{dip}}\approx {10}^{13}-{10}^{14}\,{\rm{G}}$ (see, e.g., Turolla et al. 2015). On the other hand, some sources with ${B}_{\mathrm{dip}}$ as low as $\approx {10}^{12}\,{\rm{G}}$ can show magnetar-like activity (the low-B magnetars; see, e.g., Turolla et al. 2011 and references therein). According to our simulations, the resistive tearing instability appears on a timescale $\lesssim {\tau }_{H}\approx {10}^{4}\,\mathrm{yr}$ and lasts for about $\approx 1000\,\mathrm{yr}$. This mechanism can hence provide a viable explanation for the activity (bursts and outbursts) detected in young sources (age $\lesssim {10}^{4}$ yr), which are the vast majority of the magnetar population.⁹ Whether such an instability can be triggered under the conditions typical of older objects, like the low-B sources SGR 0418+5729 and Swift J1822.3–1606 (age $\approx {10}^{5}$–10⁶ yr; Turolla et al. 2011; Rea et al. 2012), or the onset of outbursts is produced by a different, possibly related, mechanism is an open question.

Crescent-shaped features and observations—Nonpolar, crescent-like hot spots have been recently detected in NICER X-ray observations of the millisecond pulsar PSR J0030+0451 and interpreted as due to heating from backflowing currents in a nondipolar magnetic field (Miller et al. 2019). Our results show that such features can actually form as thermal relics of past events of heat deposition even in the presence of a dipole-dominated field, provided that the crustal transport properties are properly accounted for. Even though the evolutionary history of PSR J0030+0451 is likely quite different from that of a passively cooling NS and its (dipolar) field is lower than the one used in our model, our results show that a qualitatively similar behavior of the crust may be responsible for the observed pattern even without invoking strong multipolar field components.

Battery effects and magnetar magnetospheres—In Section 4, we showed how the magnetic field created through battery effects by an external heating source can reach strong local values in a turbulent-like pattern. The existence of small-scale magnetic structures, in which the field strength is orders of magnitude higher than in the surrounding dipole, has been invoked to explain the (relatively) large energy (≈1–10 keV) of the absorption features detected in the (quiescent) emission of some magnetars, if these are interpreted to be due to cyclotron absorption/scattering onto protons, ${E}_{\mathrm{cp}}\simeq 0.6(B/{10}^{14}\,{\rm{G}})$ keV. The prototypical source is the low-field (${B}_{\mathrm{dip}}\sim 6\times {10}^{12}\,{\rm{G}}$) magnetar SGR 0418+5729, where a phase-dependent absorption feature at ∼2–10 keV was discovered in the XMM-Newton data by Tiengo et al. (2013). According to their interpretation, the line arises as radiation from a cooling spot on the star surface crosses a baryon-loaded, small (≈100 m) magnetic loop where absorption occurs. Although associating this kind of magnetic structure with those produced by the battery effect in our simulations is tempting, we warn that thermocoupling effects turn out to be less important in the case of deep heating (see Section 4.1), where the local enhancement of the magnetic field is modest.

In this work, we highlighted the perspectives that a novel three-dimensional approach can open up in the study of NS magnetothermal evolution. There are, nevertheless, some limitations that must be taken into account when interpreting our numerical results.

In fact, we had to reduce the microphysical input to a realistic yet simplified model for the computing time to be manageable. This concerns in particular the use of a simplified form for the hydrostatic equilibrium density profile of Equation (18), which was also assumed to be independent of the temperature and magnetic field, and the use of a constant τ throughout the crust. Moreover, we have chosen some strong prescriptions on thermal conductivity and heat capacity. In particular, the assumption that C_V is linearly dependent on the temperature is valid only for the electron contribution and does not take into account the contribution of the lattice. This implies that Equation (6) depends on T² only, which is a key point for the efficiency of the numerical scheme. However, such an assumption becomes questionable when the term $\propto \partial T/\partial t$ starts to dominate, as in the case of impulsive heating in Section 4.1. In particular, this affects our estimates for the duration of thermal relaxation events. In fact, the heat diffusion timescale across a length L can be estimated as (Chaikin et al. 2018)

$$\begin{eqnarray}&&{\tau }_{\mathrm{diff}}\sim \displaystyle \frac{1}{4}{\left[{\int }_{L}{\rm{d}}l{\left(\displaystyle \frac{{C}_{V}}{k}\right)}^{1/2}\right]}^{2},\end{eqnarray} \tag{ 30 }$$

hence it is regulated by the specific heat-to-thermal conductivity ratio. According to the estimates of Chaikin et al. (2018), for the typical conditions of an NS, the timescale of heat transport from an internal heater to the surface is ≲1 yr, whereas in our model, the value of the characteristic diffusion time across the crust turns out to be much longer, ${\tau }_{\mathrm{diff}}\approx 50\,\mathrm{yr}$. This may well be related to our assumptions that make the ratio ${C}_{V}/\kappa$ be independent of the temperature, while it is expected to depend on the temperature as well as on the properties of the crustal superfluidity (Potekhin et al. 2015 and references therein). Hence, our results for this case should simply be regarded as indicative of the general evolution of such events. For the model discussed above, the evolution timescale is longer than what is expected under more realistic conditions by a factor of $\approx 100$, although extending this to other cases is haphazard.

Another strong prescription is that the whole physics of the core is embodied in the boundary conditions (23) and (20). Addressing the complex microphysics of the core and the description of the crust–core transition is beyond the scope of this paper (and is in general a problem best suited for one-dimensional studies). However, a direct implication of Equation (23) is that heat conduction from the crust to the core is inhibited. This is not a problem for our models but could become an issue when dealing with extremely high heat injections in the deep crust. Moreover, Equation (20) implies that the core is assumed to be in a Type I superconducting phase, so that no magnetic field is allowed to enter it, and that the magnetic flux has been completely expelled during the phase transition. However, current models of pulsar glitches (see, e.g., Baym et al. 1969) suggest that the state is of a Type II superconductor, or in any case, that at least some field is present in the core. Dealing with the modeling of this more complicated transition is again beyond the scope of this work.