© The Authors 2024

1. Introduction

Neutron stars are known to sustain the strongest magnetic fields found in the universe, ranging from 10⁸ G in the case of millisecond pulsars (e.g. Possenti et al. 2003) up to 10¹⁵ G for magnetars (e.g. Nakagawa et al. 2009). The origin of these fields still remains an open question. Spruit (2008) specifically reviewed the following two scenarios. The field could be inherited by the progenitor main sequence star through flux conservation (e.g. Makarenko et al. 2021), or it might be generated by dynamo processes before the supernova explosion. The bottom line is that they only appear as plausible explanations for pulsars but not for magnetars; other interesting theories have been suggested for the latter (e.g. Braithwaite & Spruit 2006; Mösta et al. 2015; White et al. 2022; Kiuchi et al. 2024).

Whatever the source is, we have indications that magnetic fields in neutron stars undergo evolution with time. The majority of neutron stars have strong magnetic fields on their surfaces (∼10¹⁰ − 10¹³ G), while those with stronger magnetic fields host high-energy phenomena, such flares and bursts, which have been linked to changes of their magnetic field structure (e.g. Mazets et al. 1999; Turolla et al. 2015; Kaspi & Beloborodov 2017; Younes et al. 2017; Coti Zelati et al. 2018; Lin et al. 2020; Roberts et al. 2021). The observed field strengths spanning many orders of magnitude is yet another argument towards the supposed evolution of magnetization in neutron stars given the fact that the older of them are characterized by weaker fields.

Broadly speaking, the structure of a typical neutron star may be divided into its crust and core, which in turn are further split into inner and outer parts. The outer crust is physically described by a degenerate relativistic electron gas flowing on a solid crystal Coulomb lattice of protons and ions, whereas matter consisting of neutron-rich nuclei, at densities between nuclear and the neutron drip point, constitute the inner crust. At greater depth, the outer core is a superconductive superfluid of protons and neutrons with only a few free electrons. In the innermost part of a neutron star, the inner core, the extremely high dominating densities restrict us to only speculations about exotic forms of matter, such as free quarks or other unknown particles (e.g. Baym et al. 2018). The part of a neutron star between the crust and the magnetosphere (governed by plasma) is the atmosphere, containing the ocean (sometimes called envelope). Although it might be possible that the magnetic field permeates both crust and core of a neutron star, simulations typically only considered the crustal magnetic field, and only a handful of recent studies included the core as well (e.g. Ciolfi & Rezzolla 2013; Moraga et al. 2023; Skiathas & Gourgouliatos 2024). In our paper, we also assume that the magnetic flux in the core is frozen, and we study its evolution in the crust.

There are typically three processes governing the magnetic field evolution in neutron stars (Goldreich & Reisenegger 1992). The Hall drift is the transport of magnetic flux by a moving electron fluid rising electric currents and is conservative. The conversion of magnetic energy into heat because of the crust’s finite conductivity is called Ohmic dissipation. Third, the ambipolar diffusion, another dissipative process, is caused by the relative movement of the magnetic field and charged particles to neutrons (e.g. Shalybkov & Urpin 1995). The latter effect is only important in the core, where neutrons are abundant, and hence is disregarded in the context of this study. The inner battle between Hall and Ohmic terms in neutron star crusts has been thoroughly analyzed in the literature (e.g. Hollerbach & Rüdiger 2002; Cumming et al. 2004; Pons & Geppert 2007; Kojima & Kisaka 2012; Gourgouliatos & Cumming 2014; Gourgouliatos et al. 2016). However, the magnetic and thermal interactions cannot be considered isolated in neutron star physics, and because of their nonlinear interactions, it has only recently become possible to simulate their coupling in 2D (e.g. Viganò et al. 2013) and 3D (e.g. Igoshev et al. 2021; De Grandis et al. 2020; Dehman et al. 2023; Ascenzi et al. 2024). Motivated by these efforts, in this paper, we examine a yet another possible source of magnetic field production, driven the Biermann battery mechanism (Biermann 1950).

The possible generation of magnetic field through thermoelectric instabilities was initially studied in laboratory environments for laser-heated plasmas (e.g. Stamper et al. 1971; Tidman & Shanny 1974). The Biermann battery mechanism has been proposed as a mechanism for the formation of magnetic structures by weak seed fields in astrophysics contexts as well, for instance, on cosmological scales (e.g. Gnedin et al. 2000). Later works motivated by Biermann (1950) (e.g. Roxburgh 1966) suggested that differential rotation in fluid stars could produce a misalignment in isobars and equipotentials and consequently, battery effect. This effect is generally rather small compared to the dynamo action by helical motion in a fluid which is far more effective. Nonetheless, thermoelectric effects might be important in neutron stars because dynamo action cannot occur in solid crusts.

The thermoelectric mechanism was first studied in the cores of white dwarfs by Dolginov & Urpin (1980a,b). However, as noted by Blandford et al. (1983), the mechanism should not be as relevant in white dwarfs as their crystallization occurs when they are very old and little thermal energy remains. Thus, neutron stars remain as main candidates for this phenomenon.

In neutron star crusts, the battery mechanism could work the following way. For some reason, it is possible that thermal anisotropies are formed within the star. Such temperature gradients, along with density gradients that appear naturally within depth, may circulate the electron fluid, carrying both thermal and electric currents. This thermoelectrically generated electric field acts as a battery, generating a magnetic field perpendicular to the density and temperature gradients.

The first attempt of studying such effects in neutron star crusts was by Blandford et al. (1983) who extended the work of Urpin & Yakovlev (1980) about thermogalvanomagnetic effects in degenerate stars¹. In particular, they examined two mechanisms for generation of magnetic flux, in the solid crust and in the melted envelope layer above it. For the first case, they showed that heat carried by degenerate electrons could cause thermoelectric instability in the solid crust, resulting in horizontal magnetic fields that grow exponentially with time under appropriate conditions. In the second case, they found that since pressure and density gradients must be parallel, the battery mechanism cannot work in the same way as previously. However, the liquid layers contain horizontal magnetic fields as well, heat flux can drive circulation that amplifies the field strength (possibly leading to a dynamo effect), and hence, flux may be convected from the liquid to the solid layers. If these instabilities develop, the field rapidly rises to ∼10¹² G, with additional growth saturating at ∼10¹⁴ G, the point at which non-linear perturbations describing the field evolution must be included, or when the magnetic and the lattice yield stresses are comparable.

Some years later, Urpin et al. (1986) considered unstable modes in shallower depths, consisting of free degenerate electrons and fully ionized ions in a liquid or gaseous state, and unlike Blandford et al. (1983), they did not include hydrostatic equlibrium in their considerations. They worked in the linear frame as well and showed that thermoelectric modes are unstable if the star rotates. They were able to deduce toroidal field growth in quite hot neutron stars in a very narrow region close to surface. The growth rate was found to depend on the surface temperature and the typical size of the field fluctuation (with a maximum at ∼100 m). Therefore, the generated field has a strong non-dipolar character, as large-scale fields decay in linear approximation. Conversely, Blondin & Freese (1986) suggested that the accretion process in young neutron stars may remove some of their magnetic field in a short period of time.

Next, the amplification of neutron star magnetic fields through thermoelectric effects was further analyzed extensively in a series of six papers (Geppert & Wiebicke 1991, 1995; Wiebicke & Geppert 1991, 1992, 1995; Wiebicke & Geppert 1996). In the first paper (Geppert & Wiebicke 1991), the general formalism was presented, which considered an isolated non-accreting slowly rotating (thus neglecting deviations from spherical symmetry) neutron star and examining the field growth in the limit in which the Lorentz force can be safely ignored and magnetic pressure remains negligible compared to other forms of pressure (i.e. the magnetic field can alter the thermal but not the density profile of the star). Thermoelectric amplification was sought in the liquid layer, but as opposed to Urpin et al. (1986), the outer boundary condition was slightly different, as were the mathematical methods for simultaneously solving the system of heat transport and induction equations.

In the linear approximation (Wiebicke & Geppert 1991), which is valid for moderate field strengths (≲10¹¹ G), the growth and decay rates of the magnetic field components for different multipolar components were calculated for several surface temperatures and gravitation intensities. The authors found that small-scale (n ∼1000) toroidal components grow exponentially in some years under certain conditions within a thin layer (∼20–100 m) below the surface creating sinusoidal variations with meridional wavelengths of ∼100 m (like the results obtained by Urpin et al. 1986). Large-scale modes such as dipole and quadruple slowly decay in linear approximation. These outputs were reinforced by the analytical approximations performed in Wiebicke & Geppert (1996). In this last paper, simple formulae for the maximum position and width of the temperature and field functions (assumed to be Gaussian) were derived. Interestingly, these two functions were found to be quite symmetric despite the strong dependence of the parameters on density. It was shown that the layer and maximum position of the functions loosely depend on n, but with larger n, the thickness decreases and the maximum position is shifted outwards. The authors were finally able to reproduce the described behavior of the growth rate in their analytical formula and reached a straightforward criterion for the occurrence of thermoelectric instability.

In the other papers, non-linear calculations were performed. Wiebicke & Geppert (1992) investigated the dominating non-linear terms on the toroidal small-scale magnetic field rapid growth, and yielded that the rate decreases with increasing field strength because of some form of self-coupling. Asymptotically, a stationary solution at 10¹¹ G was obtained. Saturation occurs when the characteristic time (given by the inverse of electron Larmor frequency) becomes comparable with the relaxation time for electron-ion collisions (in principle, when the electrons are forced to spiral around the field lines, and hence, collisions become improbable). Since the seed field is unknown a priori, Wiebicke & Geppert (1995) introduced averaged small-scale modes to obtain the structure of an individual axisymmetric mode, the relative weights of which are defined by initial conditions (assuming chaotic distributions). Regardless of the seed field geometry, the axisymmetric toroidal large-scale modes with even multipolarity in particular are strongly induced, while all other modes (including poloidal ones) are amplified at the same rate, and their (absolute) field values are significantly lower. Additionally, in contrast to the linear approximation, large-scale modes can be amplified rapidly in non-linear coupling with small-scale modes, and some selection rules of this coupling were identified. Wiebicke & Geppert (1996) described the evolution of the strongly induced axisymmetric large-scale toroidal field modes with even multipolarities (without considering their coupling with the poloidal field). Their evolution seems rather insensitive to the underlying star properties, and they may be split into a phase of exponential growth followed by inward diffusion along with gradually slowing down further enhancement.

In spite of these advancements and the somewhat promising results, thermoelectric effects on neutron stars have been quite overlooked in recent years. Some encouragements to explore this direction have indeed been formulated (for instance, we refer an interested reader to Geppert 2017 for a review of the magneto-thermal evolution of neutron stars) along with individual speculations about their application (e.g. Konenkov & Geppert 2001; De Grandis et al. 2020; Geppert et al. 2021). Nevertheless, a thorough numerical investigation of the battery term and its battle with the Hall and Ohmic effects has not been done. This is the goal of the current study.

On the grounds of a comprehensive exploration of the possible effects of a thermoelectric battery term on the magnetic field evolution within a neutron star crust, we have performed a number of numerical simulations solving the magnetic induction equation for different sets of the stellar configurations. We aim to identify under which conditions the battery mechanism can first be activated, and then prevail over the other two processes. Particularly, we are interested in the range of temperatures and initial magnetic field structures and intensities required to give some predictable results, as well as in the typical timescales for the relevant phenomena to evolve. Finally, we assess the extent to which the ensemble of parameters we examined may lead to conditions that can be applied to realistic neutron stars.

This paper is organized as follows. Section 2 contains the set of the underlying equations describing the problem and in Section 3, we discuss the numerical approach implemented to solve them. We continue with analyzing the results from the simulations (Section 4) and discuss their implications for realistic neutron stars in Section 5. A summary of this work is given in Section 6.

2. Mathematical setup

2.1. Magnetic induction equation

We worked in the one-fluid approximation of the crust, in which ions are considered to be static and to form a crystal Coulomb lattice, and in which electrons have the freedom to move. They are thus the carriers of electric current and thermal energy as the phonon contribution is considered negligible if the magnetic field is strong enough, stronger than ∼10¹² G, as deduced by Chugunov & Haensel (2007). We remark though that for such magnetic fields the phonon conductivity and the electron conductivity across the magnetic field are comparable, and this would enforce an anisotropic electron conductivity; nevertheless thermal conductivity is still assumed to be isotropic in this work. The typical electron velocities are low enough to resort to non-relativistic electron magneto-hydrodynamics (e-MHD). Thus, the electric current density j is given by:

$$\begin{array}{c}\hfill \mathit{j}=-e{n}_{\mathrm{e}}\mathit{\upsilon },\end{array}$$(1)

where −e is the elementary electron charge n_e is the electron number density, and υ is the mean electron velocity.

In the case of finite electric conductivity σ, the electric field E may be written through the generalized Ohm’s law:

$$\begin{array}{c}\hfill \mathit{E}=-\frac{\mathit{\upsilon }\times \mathit{B}}{c}+\frac{\mathit{j}}{\sigma }-\frac{1}{e}{S}_{\mathrm{e}}\mathrm{\nabla }T,\end{array}$$(2)

in which c is the speed of light, B is the magnetic field, T is the temperature and S_e represents the electron entropy, defined as:

$$\begin{array}{c}\hfill S_{\mathrm{e}}={\left(\frac{{\pi }^4}{3n_{\mathrm{e}}}\right)}^{1/3}\frac{k_{\mathrm{B}}^{2}T}{cħ},\end{array}$$(3)

where k_B is Boltzmann’s constant, and ℏ is the reduced Planck constant. In this paper, the last term of Eq. (2) is of particular importance. It expresses an additional component of the electric field, the component that is formed by a temperature gradient. In its general form, this electric field component is given by:

$$\begin{array}{c}\hfill \mathit{E}=\widehat{G}\mathrm{\nabla }T,\end{array}$$(4)

where the tensor $\widehat{G}$ is the thermopower and has two parts: an isotropic part (Seebeck term) and an anisotropic part (Ettingshausen-Nernst term). The Biermann battery is due to the isotropic term (De Grandis et al. 2020). This term has the form G_ij = −S_e/eδ_ij for a completely degenerate Fermi gas (Gourgouliatos et al. 2022), as written above.

The magnetic fields have currents as sources, according to Ampère’s law:

$$\begin{array}{c}\hfill \mathit{j}=\frac{c}{4\pi }\mathrm{\nabla }\times \mathit{B}.\end{array}$$(5)

Faraday’s law states that any change of the magnetic field induces an electric field as:

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}=-c\mathrm{\nabla }\times \mathit{E}.\end{array}$$(6)

Combining the above equations, we obtain the full magnetic induction equation:

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}=-c\mathrm{\nabla }\times (\frac{1}{4\pi e{n}_{\mathrm{e}}}(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}+\frac{c}{4\pi \sigma }\mathrm{\nabla }\times \mathit{B}-\frac{1}{e}{S}_{\mathrm{e}}\mathrm{\nabla }T).\end{array}$$(7)

In this equation, the first term expresses the Hall effect, which is dependent on B², and the second term refers to Ohmic dissipation (Jones 1988; Goldreich & Reisenegger 1992). The third term is the thermal battery, which is independent of B, and its strength is modulated by the T gradient and the thermopower.

We closely followed the formalism by Gourgouliatos & Cumming (2014). In axisymmetry, we decomposed the magnetic field into a poloidal (B_P) and a toroidal (B_T) component, and traced their evolution in terms of two scalar quantities, Ψ(r, θ) and I(r, θ):

$$\begin{array}{c}\hfill \mathit{B}={\mathit{B}}_{\mathrm{P}}+{\mathit{B}}_{\mathrm{T}}=\mathrm{\nabla }\mathrm{\Psi }\times \mathrm{\nabla }\varphi +I\mathrm{\nabla }\varphi ,\end{array}$$(8)

where we use spherical coordinates (r, θ, ϕ) and $\mathrm{\nabla }\varphi =\widehat{\varphi }/rsin\theta$. By construction, the above expression fulfills the Gauss law for the magnetic field, ∇ ⋅ B = 0, provided that Ψ and I are differentiable quantities. The function Ψ is referred to as the flux function (identical in form to the Stokes stream function), since 2πΨ is the poloidal magnetic flux passing through a spherical cap of radius r and opening angle θ. I is the poloidal current function, as cI/2 expresses the poloidal current passing through the same spherical cap. The magnetic field lines lie on the surfaces on which Ψ is constant, while the current or electron flow lines lie on surfaces of constant I.

We now express Eq. (7) in terms of the quantities Ψ and I. As in Reisenegger et al. (2007), we define the quantity

$$\begin{array}{c}\hfill \chi =\frac{c}{4\pi \mathrm{e}{n}_{\mathrm{e}}{r}^2{sin}^2\theta },\end{array}$$(9)

the toroidal current as

$$\begin{array}{c}\hfill {\mathit{j}}_{\mathrm{T}}=\frac{c}{4\pi }\mathrm{\nabla }\times {\mathit{B}}_{\mathrm{P}}=-\frac{c}{4\pi }{\mathrm{\Delta }}^{\ast }\mathrm{\Psi }\mathrm{\nabla }\varphi ,\end{array}$$(10)

and the angular velocity of the electrons is

$$\begin{array}{c}\hfill \mathrm{\Omega }=-\frac{j_{\mathrm{T}}}{n_{\mathrm{e}}\mathrm{e}rsin\theta }=\chi {\mathrm{\Delta }}^{\ast }\mathrm{\Psi },\end{array}$$(11)

where the Grad-Shafranov operator is given by

$$\begin{array}{c}\hfill {\mathrm{\Delta }}^{\ast }=\frac{{\partial }^2}{\partial r^2}+\frac{sin\theta }{r^2}\frac{\partial }{\partial \theta }\left(\frac{1}{sin\theta }\frac{\partial }{\partial \theta }\right)·\end{array}$$(12)

After some vector calculus manipulation, the coupled differential equations for Ψ and I are as follows. The equation related to poloidal field evolution is

$$\begin{array}{c}\hfill \frac{\partial \mathrm{\Psi }}{\partial t}-r^2{sin}^2\theta \chi (\mathrm{\nabla }I\times \mathrm{\nabla }\varphi )·\mathrm{\nabla }\mathrm{\Psi }=\frac{c^2}{4\pi \sigma }{\mathrm{\Delta }}^{\ast }\mathrm{\Psi },\end{array}$$(13)

and the equation concerning toroidal field evolution is

$$\begin{array}{cc}& \frac{\partial I}{\partial t}+r^2{sin}^2\theta ((\mathrm{\nabla }\mathrm{\Omega }\times \mathrm{\nabla }\varphi )·\mathrm{\nabla }\mathrm{\Psi }+I(\mathrm{\nabla }\chi \times \mathrm{\nabla }\varphi )·\mathrm{\nabla }I)\hfill \\ \hfill & =\frac{c^2}{4\pi \sigma }({\mathrm{\Delta }}^{\ast }I-\frac{1}{\sigma }\mathrm{\nabla }I·\mathrm{\nabla }\sigma )+\frac{c}{e}rsin\theta (\mathrm{\nabla }{S}_{\mathrm{e}}\times \mathrm{\nabla }T)·\widehat{\varphi }.\hfill \end{array}$$(14)

hence, the toroidal and poloidal components of the field are strongly coupled to each other in their evolution.

The poloidal and toroidal energies were computed as follows:

$$\begin{array}{cc}\hfill E_{\mathrm{pol}}& =\frac{1}{8\pi }{\displaystyle {\int }_V}(B_r^2+B_{\theta }^2)dV,\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill E_{\mathrm{tor}}& =\frac{1}{8\pi }{\displaystyle {\int }_V}{B}_{\varphi }^{2}dV,\hfill \end{array}$$(16)

where B_r = −r⁻²∂Ψ/∂μ, B_θ = −1/(rsinθ)∂Ψ/∂r and B_ϕ = I/(rsinθ) are the components of the magnetic field in spherical coordinates.

2.2. Crust properties

We now outline our model of the microphysical parameters within the crust, used to solve Eqs. (13) and (14). First, the crust is modeled as a spherical shell, covering the outermost 5% part of the star in terms of radius. Hence, the crust extends from the crust-core boundary, r_cc, up to the base of the ocean, r_out, and r_cc/r_out = 0.95. As the radius of our modeled star is 11.3 km, we only considered the last ∼500 m of it.

We assumed that the density throughout the crust varied as a function of depth, approximating the formula deduced by Chamel & Haensel (2008) (adopted by Karageorgopoulos et al. 2019 as well) with the following expression:

$$\begin{array}{c}\hfill \rho (r)=(1+{\left(\frac{r_{\mathrm{out}}-r}{r_{\mathrm{out}}-r_{\mathrm{cc}}}\right)}^4\frac{{\rho }_{\mathrm{cc}}}{{\rho }_{\mathrm{out}}}){\rho }_{\mathrm{out}},\end{array}$$(17)

where ρ_cc = 1.3 × 10¹⁴ g cm⁻³ is the density at the base of the crust and ρ_out is the density of the outermost point considered (bottom of the ocean). We note that ρ_out ≪ ρ_cc, thus the formula is sufficiently accurate at r = r_cc (precisely ρ(r = r_cc) = ρ_cc + ρ_out).

Crystallization occurs when the electrostatic potential energy between nuclei dominates over thermal energy. Their ratio Γ is the so-called Coulomb parameter:

$$\begin{array}{c}\hfill \mathrm{\Gamma }=\frac{Z^2{e}^2}{a_I{k}_{\mathrm{B}}T}·\end{array}$$(18)

The critical value of Γ is 175, which corresponds to the melting point, i.e. crystallization occurs for larger values (Potekhin & Chabrier 2000). We assumed a pure iron crust (Z = 26, A = 56). The ion sphere radius (also referred to in literature as interionic spacing) is:

$$\begin{array}{c}\hfill a_I={\left(\frac{4\pi n_I}{3}\right)}^{-1/3}.\end{array}$$(19)

The density ρ_out at the bottom of the ocean may be found using the numerical density n_I (for a fully ionized crust):

$$\begin{array}{c}\hfill {\rho }_{\mathrm{out}}=n_{I}mA,\end{array}$$(20)

which by construction is a function of temperature, as is evident from Eq. (18). As explained in Section 3.2, in our model we used a meridionally varying temperature, however, in order to determine ρ_out corresponding to the melting point, we used the base temperature T₀ without the fluctuation. The liquid part of a neutron star (ocean) exists for lower densities, but was not included in our considerations due to its comparably small size (∼ some cm) and the fact that e-MHD fails in this region.

In general, the electron number density n_e is given by:

$$\begin{array}{c}\hfill n_{\mathrm{e}}=\frac{Z}{A}(1-x_{\mathit{fn}})\frac{\rho }{m_b},\end{array}$$(21)

in which x_fn is the fraction of free neutrons outside ions and m_b is the atomic mass unit. Lander & Gourgouliatos (2019) performed a polynomial fitting of Z, A, and x_fn in terms of ρ utilizing the equation of state by Douchin & Haensel (2001), and their result may be approximated accurately enough by the following straightforward expression:

$$\begin{array}{c}\hfill n_{\mathrm{e}}={10}^{36}(1.5{\stackrel{\sim }{\rho }}^{2/3}+1.9{\stackrel{\sim }{\rho }}^2){\mathrm{cm}}^{-3},\end{array}$$(22)

which was also used later by Gourgouliatos & Lander (2021). The scaled density is $\stackrel{\sim }{\rho }=\rho /{\rho }_{\mathrm{cc}}$.

Regarding conductivity, we assumed that it scales as $\sigma \propto n_{\mathrm{e}}^{2/3}$, based on phonon scattering or impurity scattering (Cumming et al. 2004), being σ(r = r_cc) = 10²⁴ s⁻¹ at the crust-core boundary. It was not considered a function of temperature in order to keep our code manageable and focus on the physics of the magnetic field evolution.

The typical timescale for the Ohmic evolution is:

$$\begin{array}{c}\hfill t_{\mathrm{Ohm}}\sim \frac{4\pi \sigma L^2}{c^2}=566\mathrm{Myr}{\left(\frac{L}{11.3\mathrm{km}}\right)}^2\left(\frac{\sigma }{{10}^{24}{\mathrm{s}}^{-1}}\right),\end{array}$$(23)

whereas the Hall timescale is

$$\begin{array}{c}\hfill t_{\mathrm{Hall}}\sim \frac{4\pi \mathrm{e}{L}^2{n}_{\mathrm{e}}}{\mathit{cB}}=\frac{81.2\mathrm{Myr}}{B_{14}}{\left(\frac{L}{11.3\mathrm{km}}\right)}^2\left(\frac{n_{\mathrm{e}}}{{10}^{36}{\mathrm{cm}}^{-3}}\right),\end{array}$$(24)

and the timescale of the battery is

$$\begin{array}{c}\hfill t_{\mathrm{battery}}\sim \frac{e{L}^{2}B}{cT{S}_{\mathrm{e}}}=33.6\mathrm{Myr}\frac{B_{14}}{T_9^2}{\left(\frac{L}{11.3\mathrm{km}}\right)}^2{\left(\frac{n_{\mathrm{e}}}{{10}^{36}{\mathrm{cm}}^{-3}}\right)}^{1/3},\end{array}$$(25)

where B₁₄ is the strength of the magnetic field expressed in units of 10¹⁴ G and T₉ is the temperature in 10⁹ K. The relative importance of the Hall and Ohm terms was gauged by the ratio of these timescales, which is called magnetic Reynolds number (or Hall parameter):

$$\begin{array}{c}\hfill R_M=\frac{t_{\mathrm{Ohm}}}{t_{\mathrm{Hall}}}=\frac{\sigma |B|}{ce{n}_{\mathrm{e}}}·\end{array}$$(26)

R_M does vary through the crust and time, and in our simulations we keep it a few times of 10² at most, because when it reaches some 10³ simulations become unstable. The corresponding dimensionless parameter for the battery term is expressed as:

$$\begin{array}{c}\hfill R_T=\frac{t_{\mathrm{battery}}}{t_{\mathrm{Hall}}}=\frac{3^{1/3}cħ{|B|}^2}{4{\pi }^{7/3}{k}_{\mathrm{B}}^2{n}_{\mathrm{e}}^{2/3}{T}^2},\end{array}$$(27)

and ranged at values between ∼10⁻² and 10³ in our simulations.

3. Numerical setup

3.1. Numerical implementation

Our simulations were performed in a Fortran 90 environment. We considered the magnetic field of the star as axisymmetric, and our calculations were therefore restricted to 2D. We discretized the numerical domain into r and μ = cos θ (in order to have a better control of the grid). The resolution was typically 100 × 100 but we experimented in several cases with a double analysis to confirm the numerical convergence of our outputs.

Derivatives were numerically calculated using finite differences; a central difference scheme was used to calculate spatial derivatives, and a three- and five-point stencil for second and third derivatives, respectively. The numerical method adopted for the time intergations of Eqs. (13) and (14) was the second order Adams-Bashford method (AB2). The timestep is adaptable, according to the Courant condition (Courant et al. 1952):

$$\begin{array}{c}\hfill \mathrm{\Delta }t\le \frac{\mathrm{\Delta }r}{{\upsilon }_{max}},\end{array}$$(28)

in which Δr is the grid spacing and we set as υ_max the maximum electron plasma velocity, as we worked in e-MHD.

Three sets of boundary conditions were adopted, enclosing the numerical domain of our simulation world. Regarding the surface, we assumed a vacuum magnetic field outside of star using a multipole fit. This implies that B_ϕ = 0 at the surface (r = r_out), giving a multipolar current-free poloidal magnetic field, I(r_out, θ) = 0, so that no current flowed from the star to magnetosphere. The continuity of B_θ comes as follows. In the absence of current outside of the star, we may write ∇ × B = 0, and hence, the vacuum field can be expressed as B = ∇V, where V is a scalar function. Along with ∇ ⋅ B = 0, we get Laplace equation, ∇²V = 0. Working in axisymmetry, we decomposed the external field into a sum of multipolar expansion and obtained (Marchant et al. 2011):

$$\begin{array}{c}\hfill V(r,\theta )={\displaystyle \sum _{l=1}^{\infty }}\frac{a_l}{r^{l+1}}{P}_l(cos\theta ),\end{array}$$(29)

where P_l(cos θ) are the Legendre polynomials. We have:

$$\begin{array}{c}\hfill B_r(r,\theta )=\frac{\partial V}{\partial r}=-{\displaystyle \sum _{l=1}^{\infty }}(l+1)\frac{a_l}{r^{l+2}}{P}_l(cos\theta ),\end{array}$$(30)

and hence the coefficients are calculated as:

$$\begin{array}{c}\hfill a_l=-r^{l+2}\frac{2l+1}{2l+2}{\displaystyle {\int }_0^{\pi }}{P}_l(cos\theta )B_r(r,\theta )sin\theta d\theta .\end{array}$$(31)

Thus, B_θ is given by:

$$\begin{array}{c}\hfill B_{\theta }(r,\theta )=\frac{1}{r}\frac{\partial V}{\partial \theta }={\displaystyle \sum _{l=1}^{\infty }}\frac{a_l}{r^{l+2}}\frac{d{P}_l(cos\theta )}{d\theta }=-\frac{1}{rsin\theta }\frac{\partial \mathrm{\Psi }}{\partial r},\end{array}$$(32)

giving the boundary condition:

$$\begin{array}{c}\hfill \frac{\partial \mathrm{\Psi }(r_{\mathrm{out}},\theta )}{\partial r}=-sin\theta {\displaystyle \sum _{l=1}^{\infty }}\frac{a_l}{r_{\mathrm{out}}^{l+1}}\frac{d{P}_l(cos\theta )}{d\theta }·\end{array}$$(33)

Field lines were not allowed to cross the symmetry axis, that is the functions Ψ and I functions were set to zero on the axis of the star, i.e. Ψ(r, 0) = Ψ(r, π) = I(r, 0) = I(r, π) = 0. Concerning the inner boundary, ∂I/∂r = 0, whereas we considered two cases for the boundary condition for Ψ. For the magnetic field confined in the crust, it is Ψ(r_cc, θ) = 0 or equal to a time-independent function f, Ψ(r_cc, θ) = f(θ), when allowed to permeate to the core. The crust-confined case was motivated by a magnetic field expelled out of the core faster than thermal or magnetic timescales due to the Meissner effect assumption (e.g. Lander 2014; Passamonti et al. 2017). In the core-threading case, no magnetic field evolution was allowed in the core itself or at the crust-core boundary, because the evolution timescales of the magnetic field in the core are typically much longer than those of the crust.

3.2. Initial conditions

We examined two scenarios for the initial poloidal magnetic field. In the option in which the field is restricted in the crust, the initial poloidal flux is given by the expression:

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_1(t=0)={\mathrm{\Psi }}_0(1.05-\frac{r}{r_{\mathrm{out}}})(\frac{r}{r_{\mathrm{out}}}-\frac{r_{\mathrm{cc}}}{r_{\mathrm{out}}}){\sin }^2\theta ,\end{array}$$(34)

while it is slightly different in the alternative in which the field is assumed to thread the core (without evolving within it, however, as explained in Section 3.1):

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_2(t=0)={\mathrm{\Psi }}_0(1.05-\frac{r}{r_{\mathrm{out}}})\frac{r}{r_{\mathrm{out}}}{\sin }^2\theta .\end{array}$$(35)

In both situations, the intensity of the field is controlled by the parameter Ψ₀ (in units of 6.4  ×  10²⁴ G cm²). The poloidal current function was fixed at zero at the beginning of the simulations (I₀ = 0).

The temperature profile of the neutron star is of crucial importance for our calculations too. We chose to adopt a Gaussian-type function for the temperature distribution:

$$\begin{array}{c}\hfill T(r,\theta )=T_0+\lambda T_0\exp (-\frac{{(r-r_0)}^2}{{\sigma }_r^2})\exp (-\frac{{(\mu -{\mu }_0)}^2}{{\sigma }_{\mu }^2}),\end{array}$$(36)

which spread over both r and μ (with regard to the induction equation, the meridional gradient is mostly relevant). In this formula, a hot region centered at the point (r₀, μ₀) with a spread of 1σ (drop to a factor of 1/e) controlled by the parameters σ_r and σ_μ in radial and meridional directions, respectively, was assumed to be superimposed on a background of constant temperature T₀ (r₀ and σ_r were normalized over the simulated radius, i.e. they are expressed in units of r_out). The magnitude of the additional temperature at the point (r₀, μ₀) is λ times T₀, where λ can be any positive number determining the steepness of the temperature gradient.

We are of course aware that in most types of stars, their temperature should rise moving radially inward to their core. However, since our model considered only the outermost 5% of the neutron star (its crust), it is rather safe to assume that a uniform background temperature distribution governs this region (though it is true that any temperature gradient is rather confined precisely to the outermost layers). The additional heat region was added to account for any temperature anisotropy that might rise in the outer layers of the neutron star (the reasons for this behavior are explained more thoroughly in Section 5), which was subsequently supposed to spread throughout the star. As for the choice of the temperature function (Gaussian profile), lacking any essential evidence of its mathematical formulation, we argue that exponential decays appear rather naturally in such environments. Besides, in any tests we performed with other functions with similar profiles (e.g., quadratic) we found no signs of significant changes in our results.

We are mostly interested in the influence of the temperature gradient on the magnetic field architecture, hence the point of interest is therefore mainly limited to a region near the pole, in which the temperature is highest. For most of our simulations, we set the central point equal to (r₀ = 0.99, μ₀ = 0.99) and the standard deviations σ_r = 0.1 and σ_μ = 0.1, while unless otherwise stated we kept λ = 1 (i.e. the highest temperature is twice the background temperature). Sections 4.1 and 4.2 are dedicated to exploring the field evolution with different combinations of T₀ (high and low, accordingly) and Ψ₀. A representative temperature profile is presented in Fig. 1 (the star’s crust was scaled up by a factor of eight for visualization purposes in all such figures in this paper). On the other hand, alternative values for the remaining parameters in Eq. (36) (λ, r₀, μ₀, σ_r, σ_μ) are examined in Section 4.3.

Fig. 1. Temperature profile for λ = 1, r0 = 0.99, μ0 = 0.99, σr = 0.1 and σμ = 0.1.

Lastly, it should be remarked that throughout the simulation time-span, the temperature is kept constant in every run. In other words, we did not account for the evolution of temperature with time and assumed that some form of heat source maintained it for some reason at the same levels. While this is apparently an over-optimistic assumption (especially for the timescales considered), it is a necessary approximation we needed to make at the current time, in order to examine battery effects in neutron stars and to keep our calculations manageable. Future studies should include the thermal and magnetic interactions in a more realistic context.

4. Results

Below, we present the outputs from different models of the temperature profile. All cases presented here were simulated for a maximum time of ones million years, in order to be able to observe the saturation of the system each time.

4.1. Model A. High temperature

As stated in Section 3.2, our basic layout hd the following parameter values: r₀ = 0.99, μ₀ = 0.99, σ_r = 0.1, σ_μ = 0.1, and λ = 1. In this subsection, we examine a rather high main crust temperature, of T₀ = 10⁹ K.

We began with an initial poloidal strength of Ψ₀ = 10.0 and the case of Ψ₁, Eq. (34). Three snapshots of the evolution of the magnetic field, at 10, 100 and 1000 kyr are given in Fig. 2. The colorbar on the left refers to the value of the toroidal magnetic field, and the density of the contour lines gives an estimation of the poloidal magnetic field.

Fig. 2. Evolution of Model A through time in the case of Ψ1 and Ψ0 = 10.0. In each snapshot, the same number of poloidal lines is drawn, so that direct comparisons may be made.

Soon after the start of the simulation, a region of strong toroidal field is formed near (r₀, μ₀), which gradually overpassed in strength the structures formed by the Hall effect (yellow and light blue regions covering the two hemispheres in Fig. 2a). As time passes, the region of strong field due to the battery shifted deeper and toward the axis of the star, as a result of velocity, resistivity, and density differences between different parts of the crust. The addition of the thermoelectric current, and its resulting toroidal field creation, cause the poloidal lines to change their initial shape, and to start to twist around the battery region. These lines eventually tended to encircle the region of the strong toroidal field, creating interesting arcades and multipoles near the pole, but were kept almost intact at distant regions.

We now turn on our attention on the evolution of the poloidal and toroidal energies through time. Their resulting values are plotted in Fig. 3 as a function of time. We notice that although the poloidal energy descreased monotonically until reaching zero, the toroidal energy started from zero (since no initial toroidal field is assumed) and increased until about 2.5 × 10⁴³ erg, at which it stabilized.

Fig. 3. Evolution of toroidal and poloidal energies through time of Model A in the case of Ψ1, Ψ0 = 10.0.

The transition to constant values for the toroidal energies was located around 400 kyr (consistent with the analytical calculation of Eq. (25), with the appropriate scale), at which we assumed there is the saturation point. In other words, saturation occured when the energy and field provided by the thermoelectric battery were compensated by the Ohmic losses. After saturation was achieved, the structure and values of the toroidal field remained constant, while the poloidal field decayed with time.

We explored the evolution of the above configuration for other values of Ψ₀ as well, ranging from 0.01 up to 100.0. The final snapshots (at simulated time of 1 Myr) are presented in Fig. 4 for the cases of Ψ₀ = 0.1, 1.0, 10.0, and 100.0. Simulations with even larger Ψ₀ were not successful. Interestingly, the larger the initial poloidal field, the sparser the final poloidal lines (near the battery region), but typically, the final architecture of the lines was rather similar, especially for the low-Ψ₀ cases. Additionally, all poloidal energies gradually dropped to zero, despite starting at different values depending on Ψ₀.

Fig. 4. Final snapshots for Model A in the case of Ψ1 for different values of Ψ0 (0.1 – a, 1.0 – b, 10.0 – c, and 100.0 – d). In each snapshot, the same number of poloidal lines is drawn, so that direct comparisons may be made.

Regarding the toroidal field, it appears that the battery term dominated immediately the Hall term for lower Ψ₀, whereas it took more time to achieve that for higher values. Still, we note that saturation is observed at roughly the same time for all of them, while the final value of the maximum toroidal field in the star is also the same, 3.16  ×  10¹⁴ G, as shown in Fig. 4. This value is representative for magnetars, while the remaining star is described by values around −5  ×  10¹² G for B_T, which is two orders of magnitude lower. This small negative toroidal field is a result of the Hall term, which caused oscillations of the poloidal lines and diverts some field along the ϕ-axis; the battery term does not contribute far away from the main thermal gradient area. In addition, for all choices of Ψ₀, the toroidal energy also reached the same steady value, but for Ψ₀ ≳ 10 , the toroidal energy reaches a peak before dropping at this value.

We also experimented with even higher temperatures. For instance, setting T₀ = 5 × 10⁹ K resulted in a maximum of B_T = 2.35 × 10¹⁵ G, noticing even sparser structure of the poloidal lines. However, since even higher temperatures are not physically acceptable for the crusts of even highly magnetized magnetars, we did not consider this scenario further.

The case of Ψ₂, Eq. (35) is different. The simulations were only successful for very low initial values of Ψ₀, as in the opposite case features of an extremely large field were created, rendering the code unstable. Hence, we simulated temperatures of T₀ = 10⁹ K and Ψ₀ = 0.01 (given in Fig. 5), or even lower. The respective evolution of the energies is shown in Fig. 6. The high-B_T region is observed to shift to inner depths and smaller polar angles. At roughly the same time as in the previous configurations, saturation is again achieved (though at a slightly higher toroidal energy value, while the poloidal energy also maintains a constant value after ∼300 kyr because the poloidal flux is non-zero at the crust-core interface), apparently stabilizing the geometry of poloidal lines after ∼700 kyr. The maximum value of the toroidal field was now slighlty higher (3.27 × 10¹⁴ G) and the minimum ones around −2.50 × 10¹³ G.

Fig. 5. Evolution of Model A through time in the case of Ψ2 and Ψ0 = 0.01. In each snapshot, the same number of poloidal lines is drawn, so that direct comparisons may be made.

Fig. 6. Evolution of toroidal and poloidal energies through time of Model A in the case of Ψ2, Ψ0 = 0.01.

4.2. Model B. Low temperature

Keeping the same set of parameters as in Section 4.1 (r₀ = 0.99, μ₀ = 0.99, σ_r = 0.1, σ_μ = 0.1, λ = 1), we now considered lower crust temperatures, namely T₀ being between 10⁷ and 10⁸ K. Starting with the case of Ψ₁, Fig. 7a presents the final snapshot (1 Myr) when T₀ = 10⁸ K and Ψ₀ = 1.0. Having examined initial poloidal fields Ψ₀ in the range 0.01 to 10.0 (for lower temperatures the battery term in the resulting toroidal field can only prevail over weak poloidal fields, the final outputs were insensitive to this choice, all reaching a steady state of 2.92  ×  10¹² G at most for the toroidal field. At this low value, the toroidal field is not strong enough to twist the existing poloidal lines, hence no such feature is observed now, apart from minor oscillations after about 700 kyr. As a result, the toroidal energy rises to values of ∼10³⁹ erg, while the poloidal energy follows a similar behavior as in Model A, diminishing to zero as time passes. The remaining star (far from the battery region) was maintained at fields ∼ − 10¹⁰ G, caused by the interplay between Hall drift and Ohmic dissipation alone.

Fig. 7. Final snapshots for Model B in the case of Ψ1 and Ψ0 = 1.0 for T0 = 108 K (a) and T0 = 5 × 107 K (b).

As expected from Eq. (7), lower temperatures lead to even lower thermally induced magnetic fields. We simulated a temperature T₀ = 5 × 10⁷ K, and the final result in the case of Ψ₁ and Ψ₀ = 1.0 is shown Fig. 7b. B_T has a value of 7.2 × 10¹¹ G near the pole, which is weaker by almost an order of magnitude than the field in the case of Fig. 7a, as is the toroidal energy after saturation.

We notice that for the temperatures examined here (10⁷ − 10⁸ K), the toroidal fields are in the range of ∼10¹¹ − 10¹² G, which is representative for typical pulsars. At even lower temperatures, the mechanism becomes rather ineffective. In particular, for T₀ ≲ 10⁷ K no battery feature is apparent and the resulting magnetic field evolution is mainly described by the combination of Hall and Ohm terms.

In the field Ψ₂ scenario, the Hall term outweighs battery even at 10⁸ K because no signs of the latter are observed. We considered Ψ₀ as low as 0.01, since as explained in Section 4.1, the code cannot sustain higher values. The case of radial field assumed to penetrate the core was therefore not considered further.

4.3. Experiments with the temperature function

This subsection is dedicated to modifying the parameters kept constant (r₀, μ₀, σ_r, σ_μ, λ) and to exploring their effect on the results. Unless otherwise stated, we will adopt as our main simulation setup the one presented in Figs. 2, 3 (Model A, Ψ₁, Ψ₀ = 10.0, T₀ = 10⁹ K) and made any changes (one at a time) regarding the temperature function parameters to this one.

We first experimented with the position of the point with the highest temperature throughout the crust (referred to as the heat point). We assumed so far that the maximization of temperature occurs in r₀ = 0.99. The reason behind this option lies in the fact that magnetospheric currents (which contribute to the augmentation of temperature near the polar cap) are attenuated within the first ∼100 m of the star (at a depth corresponding to approximately 1% of the stellar radius). While this seems a realistic choice, heat might accumulate at other points inside of this. Consequently, we ran two tests with r₀ = 0.95 and r₀ = 0.90. In the former case, the results appear to be almost identical compared to when r₀ = 0.99 in terms of the final values and the architecture of the toroidal and poloidal field, their energies, and saturation times. Setting r₀ = 0.90 (the heat point is located at the base of the crust) some changes occured, however. The final snapshot, at 1 Myr, is given in Fig. 8, in which the same number of poloidal lines as in Fig. 2 are drawn for comparison. The highest B_T values appear to have dropped to almost half (∼1.68 × 10¹⁴ G), as does the toroidal energy (about 1.0 × 10⁴³ erg). By examining the time evolution, we additionally conclude that saturation occurs approximately 100 kyr later.

Fig. 8. Final snapshot for Model A in the case of Ψ1 and Ψ0 = 10.0 for r0 = 0.90.

We simulated the symmetric configuration for the position of the heat point in the meridional (symmetric with respect to the equator), that is, μ₀ = −0.99. Fig. 9 shows the final output (again the same number of poloidal lines as in Fig. 2 are drawn). Clearly, the heat point located at the south pole yields the opposite values (i.e., the battery region is surrounded by strong negative magnetic fields), but the absolute values remain the same. Moreover, the same occurs with the other results (saturation time and magnetic energies). Despite that, moving the heat point to other latitudes (e.g., to the equator or to mid-latitudes) was not possible. The development of a strong toroidal field (10¹⁴ G) in regions far from the poles twists the existing poloidal lines in narrow areas very strongly, which prevents the code from working properly.

Fig. 9. Final snapshot for Model A in the case of Ψ1 and Ψ0 = 10.0 for μ0 = −0.99.

We then tested the spread of the Gaussian function in r and μ. By simulating an idealized neutron star, we have the freedom to vary the values σ_r and σ_μ, as the actual temperature profile of a neutron star is unknown. We chose to examine three additional scenarios. The first one is the most conservative, at which σ_r and σ_μ are reduced ten times, at 0.01. We then doubled these parameters, at 0.2, while an even more optimistic scenario, in which these parameters were as large as 0.5 each, was examined as well.

These three cases, along with the default case, are presented in Figs. 10a–d at a simulated time of 1 Myr (all of the figures have the same number of poloidal lines). In principle, the small yellow region of strong toroidal field is shifted, as expected, because different spreads also change the point of the highest temperature gradient. Low standard deviations only restrict the battery region near the pole, but with multipoles evident farther from it. The highest B_T values look similar for σ_r, σ_μ equal to 0.01 and 0.2. The evolution in the latter case reveals that around 400 kyr, the poloidal lines start to increase in density forming wavy structures, whereas the toroidal field reaches a maximum value (1.83 × 10¹⁴ G) that has more than one region with a strong field. These fields are later unified as the field drops and saturates at values 1.58 × 10¹⁴ G at roughly 600 kyr. The large spread of the last case, Fig. 10d, along with the fact that the point of the highest temperature gradient now lies in the core, have as a consequence that the toroidal field is even lower, 7.0 × 10¹³ G, and more time (about 900 kyr) is required for saturation. Lastly, we notice that for all of these configurations, field lines become almost radial and parallel around it, because of twisting in narrow region near the initial heat point.

Fig. 10. Final snapshots for Model A in the case of Ψ1 and Ψ0 = 10.0 for 4 different values of σr and σμ: (0.01 – a, 0.1 – b, 0.2 – c, and 0.5 – d).

Moving on, we briefly discuss the parameter λ. We considered the following situation: the typical temperature is set at T₀ = 10⁸ K (Model B), which is a more realistic temperature for a neutron star, but we assumed that the temperature rose 11 times at the heat point, that is λ = 10. The final figure, with the remaining model choices being Ψ₁ and Ψ₀ = 10.0, is given in Fig. 11. It is noticed that the field values reached are almost two orders of magnitude higher (1.26 × 10¹⁴ G) than in the case of λ = 1 (Fig. 7). Hence, the battery is now capable of twisting the poloidal field. We argue that the scenario in this simulation is quite realistic, as it does not require a high temperature along the crust, but only a region that may be greatly enhanced.

Fig. 11. Final snapshot for Model B in the case of Ψ1, Ψ0 = 10.0 and T0 = 108 K for λ = 10.

While it may be tempting to experiment more with variations of these parameters, or combinations of them, and how they may influence the architecture and evolution of the thermoelectrically induced magnetic field, we restricted ourselves to the experiments performed above. This is because the temperature distribution of neutron star has not been conclusively determined so far, and hence, any deep investigation of the chosen model parameters might be too theoretical, not resembling any particular relevance for similar future studies. The implications, practicability, and robustness of our model and results are discussed in the next section.

5. Discussion

5.1. Analysis

To further understand the thermoelectric mechanism in neutron star crusts, we now focus on simplified cases of our model and attempt to specify any general behavior of how the system evolves. We initially examined the case of a zero magnetic field as an initial condition, that is apart from I₀ = 0, we set Ψ₀ = 0 as well. Hence, we isolated the influence of the poloidal field in the growth and resulting architecture of the toroidal one, as an initial purely toroidal field does not generate any poloidal field, whereas the poloidal field can generate a toroidal field. The final outcome (1 Myr) is given in Fig. 12a, for which all parameters other than Ψ₀ were kept the same as in Model A (Fig. 2). In both cases, the toroidal field grows linearly in time and saturates at the same values and timescales, as expected by Eq. (25). However, in the absence of poloidal field, there is no field that can been twisted along the ϕ-direction and generate Hall waves. Because of this, any toroidal field needs to be advected through the crust, and as this is a rather slow process, B_T remains zero far from the battery region instead of getting small negative values as occurs even with the lowest non-zero poloidal field.

Fig. 12. Final snapshots for Model A in the case of Ψ1 and Ψ0 = 0 with the Hall term (a) and without the Hall term (b).

We further simplified this scenario by omitting the Hall effects as well (Fig. 12b). The maximum strength of the toroidal field achieved was slightly lower in this case, 3.02 × 10¹⁴ G, though its structure was modified. Compared to the case of including the Hall term in the induction equation (Fig. 12a or Fig. 2c), the region of the strong field appears to be more compact and shifted toward the equator. This behavior is of course expected, since, for a purely toroidal field, the Hall drift tends to move electrons along the lines of constant χ (Reisenegger et al. 2007, Fig. 1 therein), depending on the field’s polarity. Consequently, when Hall effects are considered, the crust regions governed by strong fields move toward the pole and spread. We confidently determined that the high-B_T is transported toward the axis because of Hall drift, while the higher conductivity at larger depths causes the currents to move near the core, as an out-turn of Ohmic losses. Moving the gradient symmetrically in the southern hemisphere of the star, we obtained equivalent outcomes.

Another useful conclusion drawn by this simplified simulation is related to the comparative strengths of the involved terms in the induction equation. Specifically, although the Hall effect does play a role in determining the final structure and strength of the field, its contribution is rather small compared to the other two factors, and may be ignored as a first approximation. It does not affect the energy injected in the system and can only moderately modify the field configuration around the gradient region. As expected, the Hall effect is related with the transient behavior, which at sufficiently long time will fade away, whereas the final magnetic field geometry and intensity is primarily determined by the interplay between the battery term, which acts as a source, and the Ohmic diffusion.

The battle between the battery and Ohm dissipation reached a steady-state of equilibrium in all our simulations. The resulting field strengths and energy are dependent on the temperature and the temperature gradient. For instance, comparing Figs. 7a and 7b, one can see that the final toroidal field is four times higher when the maximum temperature is two times higher, while their profiles are the same. In other words, the final configuration of the crust’s magnetic field is not insensitive to this input in particular. To elucidate, we have shown that for the gradients assumed in this work, the process becomes ineffective below ∼10⁷ K. Apart from that, the poloidal energies (when including poloidal fields) dwindle exponentially in the course of time. This implies that the toroidal field is not capable of sustaining some of the poloidal field around it; the contribution of the poloidal field to the phenomenon is almost insignificant. All these remarks lead to the realization that the final toroidal field culminates a saturation equilibrium (under the proper conditions of activating the battery).

This saturation toroidal field can be approximately deduced semi-analytically. Based on Eq. (14), we derived a scaling relation at saturation, ∂I/∂t = 0. Ignoring the term involving the poloidal field (Ψ does not play a particular role compared to I), the final expression writes as:

$$\begin{array}{cc}& I(\mathrm{\nabla }\chi \times \mathrm{\nabla }\varphi )·\mathrm{\nabla }I\hfill \\ \hfill & =\frac{c^2}{4\pi \sigma }({\mathrm{\Delta }}^{\ast }I-\frac{1}{\sigma }\mathrm{\nabla }I·\mathrm{\nabla }\sigma )+\frac{c}{e}rsin\theta (\mathrm{\nabla }{S}_{\mathrm{e}}\times \mathrm{\nabla }T)·\widehat{\varphi }.\hfill \end{array}$$(37)

The term on the left-hand-side corresponds to the Hall non-linear advection of the toroidal field (Vainshtein et al. 2000), whereas the other two terms correspond to Ohmic losses and the battery. Figs. 13a-c show the strength of each term (Hall, Ohmic, and battery) throughout the crust at the final configuration step (1 Myr) for the model of Fig. 12a (Model A, Ψ₁, Ψ₀ = 0). In the high-B_T region, the battery term is dominant, reaching values around 1.3 × 10⁸, and the Ohmic term is about −9.5 × 10⁷. They sum up to 3.6 × 10⁷ representing the Hall term. The latter is three and two times weaker than the former terms, respectively (absolute values).

Fig. 13. Strength of the Hall (a), Ohmic (b) and battery (c) terms as deduced in Eq. (37) for Model A in the case of Ψ1 and Ψ0 = 0 at 1 Myr (final step). The outermost radial grid layer of b panel (Ohmic term) was masked out for visualization purposes, due to the large values in this region (∼6000).

We can use dimensional analysis to estimate the scale of the saturation field, by keeping only the Ohmic and battery terms. The toroidal magnetic field obtained from I writes as:

$$\begin{array}{cc}\hfill B_{\varphi }& =\frac{I}{rsin\theta }\sim \frac{S_{\mathrm{e}}T{L}_{S}4\pi \sigma }{\mathit{ec}}=\frac{k_B^2{T}^{2}4\pi \sigma }{c^{2}eħ}{\left(\frac{{\pi }^4}{3n_{\mathrm{e}}}\right)}^{1/3}\hfill \\ \hfill & =17B_{14}{T}_9^2{\left(\frac{n_{\mathrm{e}}}{{10}^{36}{\mathrm{cm}}^{-3}}\right)}^{-1/3}\left(\frac{\sigma }{{10}^{24}{\mathrm{s}}^{-1}}\right)\hfill \end{array}$$(38)

where L_S is the scale at which the temperature and the magnetic field change, assuming that the battery term will activate currents within the region of the temperature anisotropy. The estimate of the magnetic field obtained through dimensional analysis is higher by a factor of 5 than the one typically found through simulations, and this can be attributed to the slight offset of the spatial structure of the magnetic field from that of the temperature. To quantify the impact of the Hall non-linear advection term we can use the results of the runs shown in Figs. 12a and b, where the impact of the Hall term affects the structure of the configuration, but not the maximum value of the toroidal field. In this sense, there is an analog of the Hall attractor, where the final outcome of the run converges to a magnetic field structure, where the toroidal field is determined by the interplay of the battery term, Ohmic dissipation and the Hall term. We obtain very similar results regardless of the initial field structure, and quite remarkably, even when no magnetic field is present in the initial conditions.

5.2. Model practicability and limitations

It is possible that temperature anisotropies are formed near the surface of a neutron star. Some possible energy sources that create a hotter region in a specific part of the crust are the following. First, it can be a result of the thermo-magnetic evolution of the star. For instance, the magnetic field itself might form hotspots through phenomena like reconnection. Additionally, the anisotropic heat conductivity due to the magnetic field can create horizontal temperature gradients, which is particularly interesting as it provides another coupling between the temperature and the magnetic field, so that there will be feedback in both directions. Second, it can be generated by backflowing magnetospheric currents. The current returning from the light cylinder forms a closed loop, and it may heat the region near the rim of the polar cap as it enters through the outer crust layers. Third, another possible mechanism only relevant to non-isolated neutron stars is accretion. Any accretion (e.g., through a disk) of different material (e.g., hydrogen) other than the iron-nickel onto the crust could lead to thermonuclear reactions, which increases the local temperature as well.

Whatever the mechanism is, it is believed that temperature differences may appear rather naturally on neutron stars. Numerous sources of strong magnetic fields and ages in the range of 10 kyr have hotspots whose temperatures exceed 5 × 10⁶ K (Viganò et al. 2013). For instance, in the magnificent seven, thermal pulses from kilometer-sized areas are observed (Haberl 2007; Turolla 2009; Potekhin et al. 2015; De Grandis et al. 2022). These pulses correspond to temperatures of 10⁶ K while the characteristic ages of the sources are in the order of Myrs. Moreover, J0030+0451 has been thermally mapped in X-rays (Riley et al. 2019; Miller et al. 2019), and both circular and linear thermal regions of ∼10⁶ K were identified near its south pole. While these temperatures are much smaller than those considered in this work, this source is rather old. Thus, surface temperature anisotropies can be maintained for rather long timescales.

Thermal blanketing (Gudmundsson et al. 1983; Beznogov et al. 2021) is another important point. The temperature in the interior of the crust can be much higher than the surface temperature, as also assumed by thermal evolution calculations (Viganò et al. 2013; De Grandis et al. 2020). There, the surface temperature is scaled to the crustal one by application of the following formula (Igoshev et al. 2021):

$$\begin{array}{c}\hfill \frac{T_{\mathrm{c}}}{{10}^8\mathrm{K}}={\left(\frac{T_{\mathrm{s}}}{{10}^6\mathrm{K}}\right)}^2,\end{array}$$(39)

where T_c is the temperature at the top of the crust, and T_s is the observed temperature at the surface. Thus, an observed temperature 5 × 10⁶ K, which is the case in several sources, corresponds to a crustal temperature of 2.5 × 10⁹ K, which lies well within the range of our calculations. We note though that the above relation of determining the inner crust temperature does not hold true when the heat comes from the surface hotspot, and it is only relevant for internal heat origin.

In our model, we assumed that temperature anisotropies appear in an idealized neutron star crust, ignoring their origin. Nevertheless, the main limitation of our considerations lies in the timescale on which this temperature gradient is maintained. We implicitly assumed that the thermal energy of the star rises in some way and remains constant over a period of hundreds of thousands of years. This implies that the temperatures we consider are somehow conserved by an external mechanism, although no such term exists in our equations (including the heat conduction in our set of equations is beyond the scopes of this paper). Realistically, temperature anisotropies drop, e.g. by neutrino emission (Yakovlev et al. 2001), over the course of time (e.g., outburts at any point of the crust typically last for 5–10 years, or accretion may last for some thousand years). While we acknowledge this limitation, we think our study will be crucial first step in integrating battery terms in equations regarding the magnetic field evolution in neutron stars. Having calculated the maximum possible impact of a battery term, we show that at least for an idealized environment, thermoelectric effects are rather important and can potentially provide solutions to open problems regarding neutron star physics, hence thermo-magnetic interactions should not be omitted.

The suitability of the temperature function used was thoroughly reviewed in Section 3.2. We now discuss the repercussions of the choices of its relevant parameters. We assumed that the hot region of the star is merely connected to the polar cap. Its position in the meridional is given as a first approximation (for dipolar field):

$$\begin{array}{c}\hfill {\theta }_0={\sin }^{-1}\left(\frac{r_{\mathrm{ns}}\mathrm{\Omega }}{c}\right),\end{array}$$(40)

implying that it is related to the rotation of star (r_ns stands for the neutron star radius and Ω is the rotational angular velocity). As we do not model a specific type of neutron star and restrict ourselves to an idealized situation, we did not underwent into much detail into the exact position of the central point, and placing it near the pole(s) will suffice. Other than that, the temperature gradient may be formed by outbursts as well, as explained, giving us the freedom to experiment with the geometry of the temperature profile. In this sense, that is why we examined a wide range for the extent of the temperature anisotropy.

Lastly, other minor limitations of our considerations are related to the models adopted for density and conductivity. Particularly, conductivity is a function of density, which in turn is dependent on the ocean temperature; thus, conductivity is that way only indirectly modulated by temperature (it also depends on temperature, at least as affected by phonon scattering). We performed a test simulation assuming the conductivity linearly dependent on temperature, without noticing any significant changes in our outputs. Later work may adopt more accurate microphysics, including both electrical and thermal conductivities regulated by the respective temperature at each point.

6. Conclusions

Under appropriate conditions, magnetic field may be generated and enhanced through battery mechanisms, one of which is by temperature anisotropies causing thermal flux. High temperature gradients are required to achieve that, and neutron star crusts are a potential candidate for such a phenomenon. In this paper, we considered an idealized neutron star and showed that magnetic field enhancement might be possible through thermoelectric effects.

In our model, we made the unrealistic approximation that a temperature profile was kept constant throughout the simulated time in order to assess the maximum effect the battery term has on the magnetic field evolution. Regions of localized strong toroidal field were formed (ranging from 10¹¹ − 10¹⁵ G, depending on the temperature assumed), around of which poloidal field twists creating multipoles. For high temperatures (≳10⁸ K), the thermoelectric effect dominates the Hall drift, and the system requires some thousand years for Ohmic losses to be in equilibrium with the battery, with only minor oscillations observed afterward. Hall advection may only influence the final configuration of the toroidal field, rather than determine its value.

Later studies should shed more light into the actual role of thermoelectric effects in neutron stars. More realistic models need to be adopted, that is considering temperature decay and spread of heat through time or taking into account the thermo-magnetic interactions. Follow-up experiments with microphysical parameters and the temperature profile can be performed in the future as well, using the most up-to-date data deduced by observations. Lastly, we expect that 3-D simulations instead of axisymmetry might be worth investigating.

Acknowledgments

We are grateful to an anonymous referee for constructive criticism that significantly improved the quality of this work. We thank Dimitris Ntotsikas for his assistance in the optimization of the code used for this work and Davide De Grandis for his insightful comments on the manuscript. KNG acknowledges discussions during the annual meeting of Simons Collaboration on Extreme Electrodynamics of Compact Sources, New York, 29 February to March 1st 2024. This work was supported by computational time granted from the National Infrastructures for Research and Technology S.A. (GRNET S.A.) in the National HPC facility – ARIS – under project ID pr015026/simnstar”.

References

All Figures

	Fig. 1. Temperature profile for λ = 1, r0 = 0.99, μ0 = 0.99, σr = 0.1 and σμ = 0.1.
In the text

	Fig. 2. Evolution of Model A through time in the case of Ψ1 and Ψ0 = 10.0. In each snapshot, the same number of poloidal lines is drawn, so that direct comparisons may be made.
In the text

	Fig. 3. Evolution of toroidal and poloidal energies through time of Model A in the case of Ψ1, Ψ0 = 10.0.
In the text

	Fig. 4. Final snapshots for Model A in the case of Ψ1 for different values of Ψ0 (0.1 – a, 1.0 – b, 10.0 – c, and 100.0 – d). In each snapshot, the same number of poloidal lines is drawn, so that direct comparisons may be made.
In the text

	Fig. 5. Evolution of Model A through time in the case of Ψ2 and Ψ0 = 0.01. In each snapshot, the same number of poloidal lines is drawn, so that direct comparisons may be made.
In the text

	Fig. 6. Evolution of toroidal and poloidal energies through time of Model A in the case of Ψ2, Ψ0 = 0.01.
In the text

	Fig. 7. Final snapshots for Model B in the case of Ψ1 and Ψ0 = 1.0 for T0 = 108 K (a) and T0 = 5 × 107 K (b).
In the text

	Fig. 8. Final snapshot for Model A in the case of Ψ1 and Ψ0 = 10.0 for r0 = 0.90.
In the text

	Fig. 9. Final snapshot for Model A in the case of Ψ1 and Ψ0 = 10.0 for μ0 = −0.99.
In the text

	Fig. 10. Final snapshots for Model A in the case of Ψ1 and Ψ0 = 10.0 for 4 different values of σr and σμ: (0.01 – a, 0.1 – b, 0.2 – c, and 0.5 – d).
In the text

	Fig. 11. Final snapshot for Model B in the case of Ψ1, Ψ0 = 10.0 and T0 = 108 K for λ = 10.
In the text

	Fig. 12. Final snapshots for Model A in the case of Ψ1 and Ψ0 = 0 with the Hall term (a) and without the Hall term (b).
In the text

	Fig. 13. Strength of the Hall (a), Ohmic (b) and battery (c) terms as deduced in Eq. (37) for Model A in the case of Ψ1 and Ψ0 = 0 at 1 Myr (final step). The outermost radial grid layer of b panel (Ohmic term) was masked out for visualization purposes, due to the large values in this region (∼6000).
In the text