© The Authors 2025

1. Introduction

Neutron stars (NSs) are categorized into different classes (e.g., magnetars, rotation-powered pulsars, rotating radio transients, isolated thermal emitters) based on their observable properties (Ascenzi et al. 2024a). Efforts toward a “grand unification” (Kaspi 2010) of their vast observational phenomena have suggested that the magnetic field serves as a unifying factor (Harding 2013; Viganò et al. 2013). Consequently, the evolution of their magnetic field continues to be a relevant and active area of research. (For a recent review, see Igoshev et al. 2021.)

The current picture of magnetic field evolution in the fluid cores of NSs is as follows: Minutes after their formation, the charged particles in their interior (mostly protons, p; electrons, e; and muons, μ), which feel the magnetic force, are strongly coupled to each other and to neutrons, n, by collisions, and the particles essentially behave as a single fluid. However, their joint motion, which is driven by the magnetic force, is strongly hindered by fluid forces (pressure and gravity) that arise from stable stratification (Pethick 1992; Reisenegger & Goldreich 1992; Goldreich & Reisenegger 1992).

During this stage, sound and Alfvén waves propagate and are eventually damped, allowing the star to (very quickly) settle into a state of hydromagnetic quasi-equilibrium. The nature of such equilibria has been extensively studied in the literature, yielding geometries that vary from ordered axisymmetric “twisted-torus” configurations with poloidal and toroidal components that stabilize each other (Braithwaite & Spruit 2004, 2006; Braithwaite & Nordlund 2006) to disordered non-axisymmetric equilibrium magnetic fields (Braithwaite 2008; Becerra et al. 2022a), depending on the initial radial distribution of magnetic energy in the star (Braithwaite 2008) and the main energy dissipation mechanism invoked, namely, viscosity versus magnetic diffusivity (Becerra et al. 2022a).

While the star cools, a crust is formed; thus, the magnetic evolution continues in two different regions, a solid crust and a fluid core, with the transition between them occurring at approximately half of the nuclear saturation density. As discussed by Goldreich & Reisenegger (1992), the two long-term mechanisms that promote the evolution of the field in the crust are Ohmic diffusion (a dissipative effect due to electron-ion collisions) and Hall drift (non-dissipative advection of the magnetic field by the electrical current). The evolution of the magnetic field in this region, its effect on thermal evolution, and its potential impact on observations have been extensively studied in the literature (see, for instance, Hollerbach & Rüdiger 2002; Pons & Geppert 2007; Pons et al. 2009; Viganò et al. 2012; Gourgouliatos & Cumming 2014b; Viganò et al. 2013; Gourgouliatos & Cumming 2014a; Marchant et al. 2014; Lander & Gourgouliatos 2019; Viganò et al. 2021). Recent efforts have even involved full three-dimensional simulations (Wood & Hollerbach 2015; Dehman et al. 2023; Ascenzi et al. 2024b).

The core, on the other hand, follows a different picture. Young NSs have high temperatures, which in principle allow the core matter to overcome stable stratification while still moving as a single fluid strongly coupled by collisions. At such temperatures, the chemical imbalances due to the magnetic force are continuously reduced by weak interaction processes (Urca reactions). This erodes the initial quasi-equilibrium state, allowing the magnetic field and particles to move through a continuous sequence of quasi-equilibrium states as the chemical imbalances are reduced (Reisenegger et al. 2005; Reisenegger 2009; Ofengeim & Gusakov 2018). However, this so-called strong-coupling regime is short-lived. The star passively cools down due to Urca reactions that are faster than reheating from the feedback caused by the magnetic evolution (even for magnetar-strength fields), quickly rendering Urca reactions inefficient at restoring chemical equilibrium. Moreover, the timescale of magnetic field evolution increases very quickly as the temperature decreases (t_B ∝ T⁻⁶); thus, the magnetic field is found to evolve very little at this stage (Moraga et al. 2024). As noted by Moraga et al. (2024), the latter implies that the kind of barotropic equilibrium configurations often used in the literature as initial configurations for proto-NSs (e.g., Thompson & Duncan 1996; Lander et al. 2021), which are solutions of the “Grad-Shafranov” (GS) equation (see Sect. 7.1 for details), cannot be reached in this early stage of evolution. However, as shown by Castillo et al. (2017, 2020), GS equilibrium configurations are expected to appear at later times in the so-called weak-coupling regime, which we discuss next.

As the star cools further, the collisional coupling between charged particles and neutrons is reduced, allowing for a motion of charged particles relative to neutrons in a dissipative process called ambipolar diffusion (Pethick 1992; Goldreich & Reisenegger 1992). The velocity of this relative motion is proportional to the imbalance between the magnetic force and the charged particle fluid forces it induces, and it is also inversely proportional to the collision rate between charged particles and neutrons. Such motion erodes the initial hydromagnetic quasi-equilibrium, promoting the evolution of the magnetic field and fluid forces as they continually readjust and move through a continuum of successive quasi-equilibria. In Castillo et al. (2020), we report two-fluid numerical simulations of ambipolar diffusion considering a core of normal (non-Cooper-paired) npe matter in axial symmetry. We showed that as ambipolar diffusion acts, the fluid forces due to neutrons and charged particles adjust, with the former going to zero and the latter eventually fully balancing the magnetic force. Thus, an equilibrium is reached that can only be very slowly eroded by the combination of Hall drift and Ohmic dissipation. Since charged particles are fully decoupled from neutrons, this corresponds to a single-fluid barotropic equilibrium configuration, which is a solution of the GS equation.

Ambipolar diffusion is arguably the main effect promoting the magnetic field evolution in the core (in that region, it is certainly more effective than Hall drift and Ohmic decay, as conductivity is much higher). Consequently, it has been proposed as a plausible alternative mechanism explaining the weak magnetic fields of millisecond pulsars (Cruces et al. 2019). It is also sometimes invoked to explain the activity of magnetars, as it strongly depends on the magnetic field strength (Thompson & Duncan 1995, 1996; Mereghetti 2008; Beloborodov & Li 2016; Tsuruta et al. 2023), although recent evidence suggests that it is difficult for it to account for the high thermal luminosity of some of these stars, unless their internal magnetic field is extremely strong (Moraga et al. 2025). Thus, its impact on the long-term evolution of the field is a matter of interest.

The method to self-consistently solve the relevant equations for npe matter (see Sect. 2) was introduced by Gusakov et al. (2017) and later refined by Ofengeim & Gusakov (2018). The latter authors presented a detailed scheme to semi-analytically solve the relevant system of equations for a given initial magnetic field configuration and demonstrated that, in general, the bulk motion of the particle fluids can be much larger than the relative velocity between species. This finding was confirmed by the numerical simulations of the magnetic field evolution presented in Castillo et al. (2020). Taking this into account makes the evolution significantly faster than previously estimated (e.g., Goldreich & Reisenegger 1992; Thompson & Duncan 1995; Beloborodov & Li 2016; Castillo et al. 2017; Passamonti et al. 2017).

The semi-analytical scheme proposed by Gusakov et al. (2017) and Ofengeim & Gusakov (2018), while useful to calculate the velocities and chemical potential perturbations induced by an analytically given magnetic field, has limited applicability in long-term numerical simulations, as the latter involve iterating over magnetic field data stored in tables, and these data are not given as explicit analytical expressions. The authors mention that using magnetic field configurations in the form of tables in their scheme produces unreliable results; this limitation is attributed to the high number of spatial derivatives involved in the calculation (see Sect. 2), which yield numerical errors that accumulate if iterated over time (private communication).

Alternatively, in Castillo et al. (2020), we followed a different approach. Instead of directly solving the relevant equations, we introduced a fictitious friction (FF) force acting on the neutrons. The FF approach can be regarded as a numerical trick that greatly simplifies the numerical calculation of the relevant velocities and yields a very good approximation of the solution of the relevant equations (see a detailed discussion in Sect. 4). This approach allowed us to report the first two-fluid axially symmetric numerical simulations of ambipolar diffusion, correctly capturing the bulk motion of neutrons and charged particles, which cannot be done with the one-fluid approximation usually found in the literature (e.g., Goldreich & Reisenegger 1992; Castillo et al. 2017; Igoshev et al. 2023). While definitely a step forward, that work has some features that may not be entirely convincing to some readers:

1. The approach used involved two timescales that are very different for realistic magnetic fields. To capture both in the simulations, we had to consider unrealistically high values of the magnetic field (≳10¹⁷G).

2. The FF approach is conceptually not trivial, and at first glance, it is not evident that it yields the correct velocity fields.

3. The toy model equation of state (EoS) used (see Sect. 5) correctly captured the qualitative impact of stable stratification on the long-term magnetic evolution. However, the question of whether considering more realistic models would significantly affect the reported results remains.

In this paper, we compare different methods to solving the relevant equations, including the FF approach and semi-analytical approaches based on the explicit scheme of Gusakov et al. (2017) and Ofengeim & Gusakov (2018) as well as a pseudo-spectral method. Moreover, we present the results of an updated version of the code used in Castillo et al. (2020) that addresses the issues in items 1 and 3 of the above list, and we present arguments to convince the reader about the correctness and applicability of the FF approach (thus addressing item 2 in the same list).

This paper is organized as follows. In Sect. 2, we discuss the physical model and construct the equations to be solved in axial symmetry. In Sect. 3 and Sect. 4, we describe different approaches to obtaining the relevant velocities and chemical potential perturbations. In Sect. 5 and Sect. 6, we present the different stellar and magnetic field models used in this work. In Sect. 7, we describe the relevant equations for the magnetic field evolution, discuss the relevant timescales, and outline the mechanisms that promote the dissipation of the magnetic field. Our results are presented in Sect. 8, where we compare the different approaches to obtain the neutron velocity and use the FF approach to perform long-term time-evolving simulations. We check if the simulations are in agreement with the expected timescales and final equilibrium configurations. Finally, in Sect. 9, our results are summarized and our conclusions are outlined.

2. Physical model

Following the discussion of Castillo et al. (2020), we model the NS core as a fluid composed of normal (not superfluid or superconducting) neutrons, protons, and electrons, in axial symmetry. In the absence of a magnetic field, these are in a spherically symmetric hydrostatic and chemical equilibrium, having number densities, n_i(r), and chemical potentials, μ_i(r) (i = n, p, e), respectively. Note that charge neutrality implies n_c ≡ n_p = n_e, and chemical equilibrium implies μ ≡ μ_n = μ_p + μ_e.

For strong enough magnetic fields, the relative velocity of protons and electrons (related to the electric current density) is much smaller than their bulk velocity (Ofengeim & Gusakov 2018). Therefore, all charged particles are assumed to move together with velocity v_c, carrying along the magnetic flux according to Faraday’s law,

$$\begin{array}{c}\hfill \frac{\partial \mathit{B}}{\partial t}=\mathbf{\nabla }\times ({\mathit{v}}_c\times \mathit{B}).\end{array}$$(1)

In contrast, the neutrons are allowed to move with a different velocity, v_n. The main challenge in simulating this process is obtaining the velocity fields v_c and v_n at each time step.

In the presence of a magnetic field, B, the magnetic force causes small perturbations to the number densities of neutrons, δn_n, and charged particles, δn_c ≡ δn_p = δn_e, which evolve following the continuity equations

$$\begin{array}{cc}\hfill \frac{\partial \delta n_n}{\partial t}+\mathbf{\nabla }·\left(n_n{\mathit{v}}_n\right)& =\lambda \mathrm{\Delta }\mu ,\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \frac{\partial \delta n_c}{\partial t}+\mathbf{\nabla }·\left(n_c{\mathit{v}}_c\right)& =-\lambda \mathrm{\Delta }\mu .\hfill \end{array}$$(3)

Here, λ is a temperature-dependent parameter; Δμ ≡ δμ_c − δμ_n; δμ_c ≡ δμ_p + δμ_e; and δμ_n, δμ_p, and δμ_e are the chemical potential perturbations of the neutrons, protons, and electrons, respectively. The terms on the right-hand side of the equations represent the net rates per unit volume at which weak interactions convert protons and electrons into neutrons, and vice versa; assuming |Δμ|≪k_BT, where k_B is the Boltzmann constant and T is the absolute temperature (Goldreich & Reisenegger 1992). In what follows, we focus on the low-temperature, “weak-coupling” regime, where λ ≈ 0. Chemical potential perturbations are related to the number density perturbations, δn_n, and δn_c, by

$$\begin{array}{c}\hfill \delta {\mu }_n=K_{\mathit{nn}}\delta n_n+K_{\mathit{nc}}\delta n_c,\end{array}$$(4)

$$\begin{array}{c}\hfill \delta {\mu }_c=K_{\mathit{cn}}\delta n_n+K_{\mathit{cc}}\delta n_c,\end{array}$$(5)

where K_ij ≡ ∂μ_i/∂n_j (i, j = c, n). The off-diagonal element K ≡ K_nc = K_cn accounts for strong interactions between protons and neutrons.

For any realistic magnetic field, the evolution timescale of such perturbations will be many orders of magnitude shorter than the magnetic evolution timescale. Therefore, in time-evolving simulations, resolving the time-evolution timescale of the density perturbations and magnetic field for a realistic magnetic field strength is computationally prohibitive. However, as shown by Gusakov et al. (2017), the time derivatives in the continuity Eqs. (2), and (3) can be safely neglected in the study of the long-term evolution, as the density perturbation of charged particles and neutrons will almost instantaneously adjust to the current configuration of the magnetic field. Thus, we reduced the continuity equations to

$$\begin{array}{cc}\hfill \mathbf{\nabla }·\left(n_n{\mathit{v}}_n\right)& =0,\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill \mathbf{\nabla }·\left(n_c{\mathit{v}}_c\right)& =0.\hfill \end{array}$$(7)

Each particle species feels fluid forces due to their weight and degeneracy pressure gradient. Namely, the fluid forces on the neutrons and charged particles are, respectively, f_i = −n_i∇δμ_i − n_i(δμ_i/c²)∇Φ, for i = n, c; where Φ(r) is the gravitational potential. We note that the unperturbed star is in hydrostatic equilibrium, satisfying

$$\begin{array}{c}\hfill \mathbf{\nabla }\mu +\mu \mathbf{\nabla }\mathrm{\Phi }/c^2=0.\end{array}$$(8)

Thus, by the introduction of the relative chemical potential perturbations, χ_i = δμ_i/μ (i = n, c), fluid forces can be rewritten in a more compact form:

$$\begin{array}{cc}\hfill {\mathit{f}}_n& =-n_n\mu \mathbf{\nabla }{\chi }_n,\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill {\mathit{f}}_c& =-n_c\mu \mathbf{\nabla }{\chi }_c.\hfill \end{array}$$(10)

Charged particles also feel the Lorentz force,

$$\begin{array}{c}\hfill {\mathit{f}}_B=\frac{(\mathbf{\nabla }\times \mathit{B})\times \mathit{B}}{4\pi }.\end{array}$$(11)

An underlying assumption is that the star very quickly reaches a state of hydromagnetic quasi-equilibrium in which all the forces acting on a fluid element are close to balance, yielding

$$\begin{array}{c}\hfill {\mathit{f}}_B-n_n\mu \mathbf{\nabla }{\chi }_n-n_c\mu \mathbf{\nabla }{\chi }_c=0.\end{array}$$(12)

Consequently, in the long term, the magnetic field evolves moving through a continuum of successive hydromagnetic quasi-equilibria. In axial symmetry, the gradients, ∇χ_i, cannot have a toroidal component; therefore, Eq. (12) imposes a strong constraint on the magnetic field that must satisfy f_B, ϕ = 0.

The relative velocity of the charged particles with respect to the neutrons is the so-called ambipolar velocity, which is obtained from the force balance on the charged particles as

$$\begin{array}{c}\hfill {\mathit{v}}_{\mathit{ad}}\equiv {\mathit{v}}_c-{\mathit{v}}_n=\frac{{\mathit{f}}_B-n_c\mu \mathbf{\nabla }{\chi }_c}{{\gamma }_{\mathit{cn}}{n}_c{n}_n}.\\ \hfill \end{array}$$(13)

Here, γ_cn parametrizes the collision rate between charged particles and neutrons.

Fluid forces are written in terms of gradients of the chemical potential perturbations (see Eqs. (9) and (10)). Thus, there are additive constants in χ_n and χ_c that do not play a role in the evaluation of the velocities. If needed, these can be fixed using integral conditions corresponding to the conservation of charged particles and neutrons, i.e.,

$$\begin{array}{cc}\hfill {\displaystyle {\int }_{{\mathcal{V}}_{\text{core}}}}\delta n_{n}d\mathcal{V}& ={\displaystyle {\int }_{{\mathcal{V}}_{\text{core}}}}\mu \frac{K_{\mathit{cc}}{\chi }_n-K{\chi }_c}{K_{\mathit{nn}}{K}_{\mathit{cc}}-K^2}d\mathcal{V}=0,\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill {\displaystyle {\int }_{{\mathcal{V}}_{\text{core}}}}\delta n_{c}d\mathcal{V}& ={\displaystyle {\int }_{{\mathcal{V}}_{\text{core}}}}\mu \frac{K_{\mathit{nn}}{\chi }_c-K{\chi }_n}{K_{\mathit{nn}}{K}_{\mathit{cc}}-K^2}d\mathcal{V}=0,\hfill \end{array}$$(15)

where the density perturbations are obtained inverting Eqs. (4) and (5), and 𝒱_core is the volume of the core.

As stated before, we restricted ourselves to axial symmetry so the magnetic field could be decomposed as

$$\begin{array}{c}\hfill \mathit{B}=\mathbf{\nabla }\alpha \times \mathbf{\nabla }\varphi +\beta \mathbf{\nabla }\varphi .\end{array}$$(16)

The scalar potentials, α(r, θ), and β(r, θ), generate the poloidal and toroidal magnetic field, respectively, and $\mathbf{\nabla }\varphi =\widehat{\varphi }/(rsin\theta )$.

At each time step, the Lorentz force, f_B, is known, but the velocity fields, v_i, and the scalar fields, χ_i, for i = n, c must be obtained from the self-consistent solution of Eqs. (6), (7), (12), and (13). Different approaches to this challenge are discussed in Sect. 3. To perform numerical calculations using such approaches, the equations must be normalized appropriately. The units in which the different variables are measured are given in Table 1.

3. Semi-analytical approaches

3.1. Explicit approach

The set of Eqs. (6), (7), (12), and (13) can be solved for an analytically given magnetic field configuration in axial symmetry for which f_B, ϕ = 0. A semi-analytical procedure to solve the set of equations was described by Gusakov et al. (2017) and more in depth by Ofengeim & Gusakov (2018). Their steps to obtain the poloidal component of the velocities and chemical potential perturbations are summarized here, adopting the notation and approximations used in the present paper.

The poloidal part of Eq. (12) implies scalar equations for the radial and polar components, respectively. For a given magnetic force, f_B, these equations fix χ_n and χ_c up to a radial function for each, which can be fixed by taking the appropriate boundary conditions. However, these conditions are not trivial and must be set for v_n. To proceed, we split the relative chemical imbalances into two components: χ_n(r, θ) = X_n(r, θ)+χ_n, 0(r), and χ_c(r, θ) = X_c(r, θ)+χ_c, 0(r), where ${\chi }_{n,0}={\widehat{\mathrm{P}}}_0\{{\chi }_n\}$ and ${\chi }_{c,0}={\widehat{\mathrm{P}}}_0\{{\chi }_c\}$. Here, we introduced the operator ${\widehat{\mathrm{P}}}_l$ that extracts the coefficient of the lth Legendre component of its argument,

$$\begin{array}{c}\hfill {\widehat{\mathrm{P}}}_l\{f(r,\theta )\}\equiv \frac{2l+1}{2}{\displaystyle {\int }_0^{\pi }}f(r,\theta )P_l(cos\theta )sin\theta d\theta ,\end{array}$$(17)

where f(r, θ) is a given function of r and θ, and P_l(cos θ) is a Legendre polynomial of degree l. The relation between these radial functions can be obtained by taking the radial component of Eq. (12) and projecting its 0th Legendre component

$$\begin{array}{c}\hfill {\widehat{\mathrm{P}}}_0\{f_{B,r}\}-n_c\mu {\chi }_{c,0}^{\prime }-n_n\mu {\chi }_{n,0}^{\prime }=0,\end{array}$$(18)

where prime (′) denotes derivative with respect to r. Replacing the decomposition of χ_c and χ_n into Eq. (12) together with the latter equation yields

$$\begin{array}{c}\hfill {\mathit{F}}_B-n_n\mathbf{\nabla }{X}_n-n_c\mathbf{\nabla }{X}_c=0,\\ \hfill \end{array}$$(19)

where ${\mathit{F}}_B\equiv ({\mathit{f}}_B-{\widehat{\mathrm{P}}}_0\{f_{B,r}\}\widehat{r})/\mu$, and $\widehat{r}$ is the radial unit vector. We note that by defining ϵ ≡ n_nX_n + n_cX_c and δ ≡ n_n′X_n + n_c′X_c, Eq. (19) can be rewritten as $\mathbf{\nabla }ϵ-\delta \widehat{r}={\mathit{F}}_B$. The solution for ϵ and δ can be easily obtained by first integrating its polar component with respect to θ. Thus,

$$\begin{array}{cc}\hfill ϵ(r,\theta )& ={\displaystyle {\int }_0^{\theta }}r{F}_{B,\theta }d\theta +ϵ(r,0),\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill \delta (r,\theta )& =\frac{\partial ϵ}{\partial r}-F_{B,r},\hfill \end{array}$$(21)

where the latter comes from the radial component and $ϵ(r,0)=-{\widehat{\mathrm{P}}}_0\{{\displaystyle {\int }_0^{\theta }}r{F}_{B,\theta }d\theta \}$ is chosen so ${\widehat{\mathrm{P}}}_0\{ϵ\}=0$. In terms of these variables, we can reconstruct

$$\begin{array}{cc}\hfill X_n& =\frac{n_c^{\prime }ϵ-n_c\delta }{n_n^2{(n_c/n_n)}^{\prime }},\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill X_c& =-\frac{n_n^{\prime }ϵ-n_n\delta }{n_n^2{(n_c/n_n)}^{\prime }},\hfill \end{array}$$(23)

which also satisfy ${\widehat{\mathrm{P}}}_0\{X_n\}={\widehat{\mathrm{P}}}_0\{X_c\}=0$.

To fully determine χ_n and χ_c, we also needed to evaluate χ_n, 0 and χ_c, 0. We note that by integrating Eq. (6) over the volume, 𝒱(r), of any centered sphere of radius, r, we obtained

$$\begin{array}{c}\hfill 0={\displaystyle {\int }_{\mathcal{V}(r)}}\mathbf{\nabla }·\left(n_n{\mathit{v}}_n\right)d\mathcal{V}=2\pi r^2{n}_n{\displaystyle {\int }_0^{\pi }}{v}_{n,r}(r,\theta )sin\theta d\theta ,\end{array}$$(24)

implying ${\widehat{\mathrm{P}}}_0\{v_{n,r}\}=0$. Analogously, from Eq. (7) we obtained ${\widehat{\mathrm{P}}}_0\{v_{c,r}\}=0$, and thus also ${\widehat{\mathrm{P}}}_0\{v_{ad,r}\}=0$. On the other hand, replacing Eq. (12) and χ_n = X_n + χ_n, 0 into Eq. (13) would yield

$$\begin{array}{c}\hfill n_c{\mathit{v}}_{\mathit{ad}}=\frac{\mu }{{\gamma }_{\mathit{cn}}}(\mathbf{\nabla }{X}_n+{\chi }_{n,0}^{\prime }\widehat{r}).\\ \hfill \end{array}$$(25)

Taking the radial component of the latter and projecting its 0th Legendre component implies

$$\begin{array}{c}\hfill {\chi }_{n,0}^{\prime }=0.\end{array}$$(26)

This, together with Eq. (18), implies

$$\begin{array}{c}\hfill {\chi }_{c,0}^{\prime }=\frac{{\widehat{\mathrm{P}}}_0\{f_{B,r}\}}{n_c\mu }.\end{array}$$(27)

These conditions define χ_n, 0 and χ_c, 0 from f_B, r up to additive constants that do not affect the calculation of the velocities. If needed, they can be fixed using Eqs. (14) and (15).

Finally, the velocities can be obtained as follows: v_ad, which is purely poloidal, can be obtained from Eqs. (25) and (26),

$$\begin{array}{c}\hfill {\mathit{v}}_{\mathit{ad}}=\frac{\mu \mathbf{\nabla }{X}_n}{n_c{\gamma }_{\mathit{cn}}},\\ \hfill \end{array}$$(28)

whereas v_n follows by writing v_c = v_n + v_ad and replacing Eq. (28) into Eq. (7), which implies

$$\begin{array}{c}\hfill \mathbf{\nabla }·\left(n_c{\mathit{v}}_n\right)+\mathbf{\nabla }·\left(\frac{\mu }{{\gamma }_{\mathit{cn}}}\mathbf{\nabla }{X}_n\right)=0.\end{array}$$(29)

By multiplying Eq. (6) by n_c, Eq. (29) by n_n, and then subtracting them, we obtained, after some algebra,

$$\begin{array}{c}\hfill v_{n,r}=-\frac{1}{n_n{(n_c/n_n)}^{\prime }}\mathbf{\nabla }·\left(\frac{\mu }{{\gamma }_{\mathit{cn}}}\mathbf{\nabla }{X}_n\right).\\ \hfill \end{array}$$(30)

The polar component can be obtained from Eq. (6)

$$\begin{array}{c}\hfill v_{n,\theta }=-\frac{1}{rsin\theta n_n}\frac{\partial }{\partial r}(r^2{n}_n{\displaystyle {\int }_0^{\theta }}{v}_{n,r}sin\theta d\theta ),\end{array}$$(31)

where the limits of the integral are chosen to be consistent with conditions arising from axial symmetry, i.e., v_n, θ(r, 0) = v_n, θ(r, π) = 0.

On the other hand, the toroidal magnetic force must remain null (∂f_B, ϕ/∂t = 0), which imposes a quite complicated condition that can, in principle, be used to obtain v_n, ϕ (see Ofengeim & Gusakov 2018). So far, no feasible numerical implementation to obtain v_n, ϕ for an arbitrary magnetic field configuration has been proposed.

In summary, the steps to obtain the poloidal components of v_n and v_c for an analytically given magnetic field satisfying f_B, ϕ = 0 are

1. Evaluate ϵ and δ from the magnetic force using Eqs. (20) and (21).

2. Evaluate X_n from Eq. (22).

3. Evaluate v_ad from Eq. (28).

4. Evaluate v_n, r and v_n, θ from Eqs. (30) and (31), respectively.

5. Evaluate v_{c, Pol} = v_{n, Pol} + v_{ad, Pol},

where the subscript “Pol” refers to the poloidal part of a vector field, ${\mathit{v}}_{n\text{, Pol}}=v_{n,r}\widehat{r}+v_{n,\theta }\widehat{\theta }$.

We note that this procedure completely defines the velocity fields everywhere in the NS core, leaving no room for imposing boundary conditions at the crust-core interface. Thus, any such conditions, for example no flux of particles through the interface (v_n, r(r = R) = v_c, r(r = R) = 0), would require constraining the magnetic field (Ofengeim & Gusakov 2018), and it is unclear whether these constraints will be maintained as the magnetic field evolves in time.

In Sect. 8.1, we show plots of the velocities, which are computed using the procedure described above. All steps are evaluated using Wolfram Mathematica¹ software, which produces involved symbolic expressions for the velocities and chemical potential perturbations in terms of the coordinates r, θ, the magnetic force, the functions n_c(r),n_n(r),μ(r),γ_cn(r), and their radial derivatives. The latter background functions are evaluated from numerical data computed from the EoS.

3.2. Pseudo-spectral approach

In addition to the semi-analytical procedure outlined in Sect. 3.1, the system of Eqs. (6), (7), (12), and (13) can be solved using a pseudo-spectral approach. An optimal implementation of such an algorithm requires rewriting some of the previous equations in a spectral-friendly manner.

Spectral codes expand scalar and vector fields using polynomial bases for the radial and angular domains within a spherical region. The angular component is typically expanded in spherical harmonics as they satisfy the regularity conditions along the axis. In axial symmetry, these reduce to the Legendre polynomials P_l(cos θ). Consequently, the chemical perturbationsread as

$$\begin{array}{c}\hfill {\chi }_i(r,\theta )={\displaystyle \sum _{l=0}^{n_{\theta }}}{\chi }_{i,l}(r)P_l(cos\theta ),(i=n,c),\end{array}$$(32)

where n_θ represents the truncation order of the polynomial basis in the angular domain, and ${\chi }_{i,l}(r)={\widehat{\mathrm{P}}}_l\{{\chi }_i\}$; thus X_i = ∑_l ≥ 1χ_i, l(r)P_l(cos θ). For simplicity, we have omitted the explicit expansion of χ_i, l(r) in radial polynomials. The solution procedure is the following:

1. First, we computed the chemical potential perturbations. Instead of using Eqs. (22) and Eqs. (23) to obtain the {χ_i, l}_l ≥ 1, we proceeded in the following manner. To obtain χ_c, we divided Eq. (12) by μn_n and then took the curl

   $$\begin{array}{c}\hfill \mathbf{\nabla }\times \left[{\left(\frac{n_c}{n_n}\right)}^{\prime }{\chi }_c\widehat{r}\right]=-\mathbf{\nabla }\times \left(\frac{{\mathit{f}}_B}{\mu n_n}\right).\end{array}$$(33)

   Taking again the curl of Eq. (33) and computing only its radial component, we obtained

   $$\begin{array}{c}\hfill {\left(\frac{n_c}{n_n}\right)}^{\prime }{\mathrm{\Delta }}_H{\chi }_c=-\widehat{r}·\mathbf{\nabla }\times \mathbf{\nabla }\times \left(\frac{{\mathit{f}}_B}{\mu n_n}\right).\end{array}$$(34)

   Here, Δ_H represents the horizontal Laplacian, which in axial symmetry can be expressed as

   $$\begin{array}{c}\hfill {\mathrm{\Delta }}_H\equiv \frac{1}{r^{2}sin\theta }\frac{\partial }{\partial \theta }(sin\theta \frac{\partial }{\partial \theta }),\end{array}$$(35)

   and its eigenfunctions are the Legendre polynomials:

   $$\begin{array}{c}\hfill {\mathrm{\Delta }}_H{P}_l(cos\theta )=-\frac{l(l+1)}{r^2}{P}_l(cos\theta ).\end{array}$$(36)

   Thus, with the aid of equations (34)–(36), we obtained χ_c, l, while χ_n, l was determined by repeating the same procedure but dividing Eq. (12) by μn_c instead. Thus, for l ≥ 1, the functions χ_i, l(r) read as

   $$\begin{array}{cc}\hfill {\chi }_{n,l}& =-\frac{r^2}{l(l+1){(n_n/n_c)}^{\prime }}{\widehat{\mathrm{P}}}_l\{\widehat{r}·\mathbf{\nabla }\times \mathbf{\nabla }\times \left(\frac{{\mathit{f}}_B}{\mu n_c}\right)\},\hfill \end{array}$$(37)

   $$\begin{array}{cc}\hfill {\chi }_{c,l}& =-\frac{r^2}{l(l+1){(n_c/n_n)}^{\prime }}{\widehat{\mathrm{P}}}_l\{\widehat{r}·\mathbf{\nabla }\times \mathbf{\nabla }\times \left(\frac{{\mathit{f}}_B}{\mu n_n}\right)\}.\hfill \end{array}$$(38)

2. Next, we obtained the velocity fields. The ambipolar velocity is the simplest to acquire, and this can be done once X_n = ∑_l ≥ 1χ_n, l(r)P_l(cos θ) has been determined (see Eq. (28)). Since n_nv_n is purely solenoidal, its poloidal component can be written as

   $$\begin{array}{c}\hfill {\mathit{v}}_{n\text{, Pol}}=\frac{1}{n_n}\mathbf{\nabla }\times \mathbf{\nabla }\times (\mathrm{\Psi }\widehat{r}).\end{array}$$(39)

   Here, Ψ(r, θ) is a scalar field. Taking the radial component of the latter and projecting its lth Legendre component yields

   $$\begin{array}{c}\hfill {\widehat{\mathrm{P}}}_l\left\{v_{n,r}\right\}=\frac{l(l+1){\mathrm{\Psi }}_l}{r^2{n}_n},\end{array}$$(40)

   where ${\mathrm{\Psi }}_l={\widehat{\mathrm{P}}}_l\left\{\mathrm{\Psi }\right\}$. This, together with Eq. (30), implies

   $$\begin{array}{c}\hfill {\mathrm{\Psi }}_l=-\frac{r^2}{l(l+1){(n_c/n_n)}^{\prime }}{\widehat{\mathrm{P}}}_l\{\mathbf{\nabla }·\left(\frac{\mu \mathbf{\nabla }{X}_n}{{\gamma }_{\mathit{nc}}}\right)\},\end{array}$$(41)

   for l ≥ 1, which can be used to compute v_n, θ from Eq. (39), whereas the term Ψ₀ is not relevant, as it vanishes from the right-hand side of the latter equation.

These equations are solved using the pseudo-spectral code Dedalus²: a flexible framework for spectrally solving differential equations, encompassing boundary, initial, and eigenvalue problems (Burns et al. 2020). It expands in terms of Jacobi polynomials and spherical harmonics in the radial and angular domains, respectively (for more details, see, e.g., Vasil et al. 2019).

To compute the velocity field, we implemented the procedure described above. We computed χ_i, l (i = n, c) for l ≥ 1 by straightforward substitution in Eqs. (37) and (38). After this, obtaining v_ad and v_n was achieved through a simple substitution in Dedalus (Eqs. (28), (39), and (41)). If needed, we could also determine χ_n, 0 and χ_c, 0. To do so, a boundary value problem is solved consisting of Eq. (18) together with conditions (26)-(27) and the integral conditions (14)-(15). In the implementation of this approach, we set the truncation order to n_r = 100 and n_θ = 71 radial and angular polynomials, respectively.

We note that the pseudo-spectral approach just described has the same limitations as the explicit approach. Namely, it does not provide a way to compute the toroidal component of v_n and v_c, which makes the approach unsuitable for time-evolving simulations, and it does not allow one to impose boundary conditions at the crust-core interface.

4. Fictitious friction approach

In Sect. 3, we showed how the set of Eqs. (6), (7), (12), and (13) can, in principle, be solved for any given magnetic field configuration. However, the numerical implementation is cumbersome, as it requires knowledge of high-order derivatives of the magnetic field (Ofengeim & Gusakov 2018). In Castillo et al. (2020), we instead introduced a FF force in Eq. (12) (see Yang et al. 1986; Roumeliotis et al. 1994; Hoyos et al. 2008, 2010; Viganò et al. 2011), defined as

$$\begin{array}{c}\hfill {\mathit{f}}_{\zeta }=-\zeta n_n{\mathit{v}}_n,\end{array}$$(42)

where ζ is a (small) adjustable parameter, which is a strategy that, to our knowledge, was first described by Chodura & Schlüter (1981). With this additional force, the force balance equation takes the form

$$\begin{array}{c}\hfill {\mathit{f}}_B+{\mathit{f}}_n+{\mathit{f}}_c-\zeta n_n{\mathit{v}}_n=0,\end{array}$$(43)

from which the neutron velocity can be determined:

$$\begin{array}{c}\hfill {\mathit{v}}_n=\frac{1}{\zeta n_n}({\mathit{f}}_B+{\mathit{f}}_n+{\mathit{f}}_c).\end{array}$$(44)

Within the framework of the FF approach, instead of solving the original system of Eqs. (6), (7), (12), and (13), it is proposed to solve a modified system consisting of equations Eqs. (6), (7), (13), and (44). This system of equations is hereafter referred to as the FF equations or the FF system, and the corresponding calculation method is called the FF approach. The main advantage of this approach is the simplicity of calculating all three components of the neutron velocity v_n using Eq. (44). We note that in the limiting case ζ → ∞, we recover the one-fluid approximation, in which neutrons are modeled as a motionless background (see, e.g., Castillo et al. 2017).

However, we wanted to know how valid the FF approach is. To answer this question, we assumed that the solution of Eqs. (6), (7), (12), and (13) obtained using the rigorous method of Sect. 3.1 is known. We denote the quantities describing this solution with the superscript ⁽⁰⁾. For example, ${\chi }_n^{(0)}$ and ${\mathit{v}}_n^{(0)}$ represent the relative neutron chemical potential perturbation and neutron velocity in the quasi-stationary approximation, respectively, for a given B.

Then, we supposed that for a sufficiently small dimensionless parameter ζ (in the units specified in Table 1), the solution of the FF system can be expressed as a series expansion in ζ. For instance, for the quantities χ_n and v_n, we can write

$$\begin{array}{cc}\hfill {\chi }_n& ={\chi }_n^{(0)}+\zeta {\chi }_n^{(1)}+{\zeta }^2{\chi }_n^{(2)}+\cdots ,\hfill \end{array}$$(45)

$$\begin{array}{cc}\hfill {\mathit{v}}_n& ={\mathit{v}}_n^{(0)}+\zeta {\mathit{v}}_n^{(1)}+{\zeta }^2{\mathit{v}}_n^{(2)}+\cdots .\hfill \end{array}$$(46)

(Similar expansions are assumed to hold for other problem parameters, e.g., f_c, f_n.) Substituting these expansions into Eq. (44) gives

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\gamma }_{\mathit{cn}}{n}_c{n}_n[\zeta {\mathit{v}}_n^{(0)}+{\zeta }^2{\mathit{v}}_n^{(1)}+\cdots ]=& [{\mathit{f}}_B+{\mathit{f}}_n^{(0)}+{\mathit{f}}_c^{(0)}]+\hfill \\ \hfill & \zeta [{\mathit{f}}_n^{(1)}+{\mathit{f}}_c^{(1)}]+\cdots .\hfill \end{array}\end{array}$$(47)

Equating terms of the same order in ζ yields an infinite set of equations relating quantities of order k and k + 1 (k = 0, 1, …). The first two of these equations are

$$\begin{array}{cc}& {\mathit{f}}_B+{\mathit{f}}_n^{(0)}+{\mathit{f}}_c^{(0)}=0,\hfill \end{array}$$(48)

$$\begin{array}{cc}& {\mathit{v}}_n^{(0)}=\frac{1}{{\gamma }_{\mathit{cn}}{n}_c{n}_n}[{\mathit{f}}_n^{(1)}+{\mathit{f}}_c^{(1)}].\hfill \end{array}$$(49)

The first equation is satisfied automatically, as it coincides with the force balance equation (12) of our exact problem. The second equation, in turn, provides a relationship between the neutron velocity ${\mathit{v}}_n^{(0)}$ of the exact solution and the first-order corrections ${\mathit{f}}_n^{(1)}$ and ${\mathit{f}}_c^{(1)}$ to the forces ${\mathit{f}}_n^{(0)}$ and ${\mathit{f}}_c^{(0)}$. It can be shown that the remaining Eqs. (6), (7), and (13) allow, in principle, for the complete calculation of all first-order corrections in terms of the zeroth-order quantities (computed using the rigorous approach). The higher-order terms can be found following the same strategy.

From this example, the essence of the FF approach becomes clear. Within this approach, the solution approximately corresponds to the exact solution of Eqs. (6), (7), (12), and (13) with an accuracy of 𝒪(ζ). For instance: ${\mathit{v}}_n-{\mathit{v}}_n^{(0)}\sim \mathcal{O}(\zeta )$. Thus, for sufficiently small ζ, the FF approach should yield a solution that deviates only slightly from the exact calculation.

We have just discussed how, given the exact solution of the system Eqs. (6), (7), (12), and (13), the solution of the FF system can be derived. However, in practice, the FF approach was specifically proposed to avoid solving the exact system of Eqs. (6)–(13). Supposing we chose a sufficiently small parameter, ζ, and found the solution of the FF equations, we wanted to determine whether this solution would always be close to the exact solution of the system (6)–(13). Unfortunately, this is not the case. As discussed below, the FF system requires more boundary conditions than the original system (6)–(13). Specifically, to uniquely solve the FF equations, some boundary conditions must be imposed, for example, by specifying the r components of the velocities v_n and v_c at the crust-core interface. Meanwhile, in the exact system, these components are determined automatically during the solution process, as they depend on the magnetic field structure. Thus, the FF system of equations describes a broader class of solutions, which only becomes close [to an accuracy of ∼𝒪(ζ)] to the exact solution if the boundary conditions for the FF system (values of the r components of the velocities v_n and v_c at the boundary) are chosen to be sufficiently close to the corresponding boundary values obtained in the exact calculation. This interpretation has been verified numerically; see Sect. 8.1.2 for details. With these caveats, the FF approach can be successfully used to determine velocity fields for a given magnetic field configuration, as described below.

To solve the FF system for a given initial condition for the magnetic field, we can replace the velocities from Eqs. (13) and (44) into the continuity Eqs. (6) and (7), together with Eqs. (9) and (10). After rearranging terms, leaving only the terms explicitly dependent on the magnetic field on the right-hand side of the equations, we obtained

$$\begin{array}{cc}\hfill \mathbf{\nabla }·(n_n\mu \mathbf{\nabla }{\chi }_n+n_c\mu \mathbf{\nabla }{\chi }_c)& =\mathbf{\nabla }·{\mathit{f}}_B,\hfill \end{array}$$(50)

$$\begin{array}{cc}\hfill \mathbf{\nabla }·[n_c\mu \mathbf{\nabla }{\chi }_n+g(\zeta )n_c\mu \mathbf{\nabla }{\chi }_c]& =\mathbf{\nabla }·\left[g(\zeta ){\mathit{f}}_B\right].\hfill \end{array}$$(51)

Here, we introduced a dimensionless radial function,

$$\begin{array}{c}\hfill g(\zeta ,r)=\frac{\zeta }{{\gamma }_{\mathit{cn}}{n}_n}+\frac{n_c}{n_n}.\end{array}$$(52)

The differential equations (50), (51), and integral conditions (14) and (15) fully define the functions χ_c and χ_n when combined with appropriate boundary conditions such as the non-penetration conditions at the crust-core interface for the velocity fields, i.e., v_n, r = v_ad, r = 0, which translates into

$$\begin{array}{cc}\hfill {\frac{\partial {\chi }_c}{\partial r}|}_{r=R}& ={\frac{f_{B,r}}{n_c\mu }|}_{r=R},\hfill \end{array}$$(53)

$$\begin{array}{cc}\hfill {\frac{\partial {\chi }_n}{\partial r}|}_{r=R}& =0.\hfill \end{array}$$(54)

We chose these boundary conditions because they can be easily implemented numerically. However, it should be noted that, as shown by Ofengeim & Gusakov (2018), for an arbitrary magnetic field, the exact solution to the set of Eqs. (6), (7), (12), and (13) will typically yield non-null radial components of the velocities at the boundary.

4.1. Numerical implementation using finite difference and finite volume

In order to numerically solve our equations using a finite difference and finite volume implementation, we discretize the relevant variables over a staggered polar grid composed of N_r points that can be homogeneously or inhomogeneously distributed in the radial direction and N_θ points equally spaced in the polar direction, following the rule

$$\begin{array}{cc}\hfill r_i& ={\left(\frac{i-1}{N_r-1}\right)}^{u}i=1,\cdots ,N_r,\hfill \end{array}$$(55)

$$\begin{array}{cc}\hfill {\theta }_j& =\pi \frac{j-1}{N_{\theta }-1}=(j-1)\mathrm{\Delta }\theta j=1,\cdots ,N_{\theta }.\hfill \end{array}$$(56)

The exponent, u, is specified later. We discretized the relevant variables on this grid in the following way, where half-integer indices denote averages (e.g., r_i + 1/2 ≡ (r_i + r_i + 1)/2):

1. The variable α was evaluated at the corner of each cell: α_i, j = α(r_i, θ_j).

2. The variable β, χ_n, and χ_c and the ϕ component of the forces and velocities were evaluated at the center of each cell: β_i, j = β(r_i + 1/2, θ_j + 1/2).

3. The radial components of vector fields, such as the magnetic field, the forces, and the velocities were evaluated at the centers of cell faces defined by a constant r = r_i: B_r i, j = B_r(r_i, θ_j + 1/2).

4. The polar components of vector fields, such as the magnetic field, the forces, and the velocities were evaluated at the centers of cell faces defined by a constant θ = θ_j: B_θ i, j = B_θ(r_i + 1/2, θ_j).

The value of the divergence of a vector field F at the center of each cell (∇ ⋅ F)_i, j = ∇⋅F|_{ri + 1/2, θj + 1/2} was given by

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {(\mathbf{\nabla }·\mathit{F})}_{i,j}=\frac{1}{{\mathcal{V}}_{i,j}}(& S_{ri+1,j}{F}_{ri+1,j}-S_{ri,j}{F}_{ri,j}\hfill \\ \hfill & +S_{\theta i,j+1}{F}_{\theta i,j+1}-S_{\theta i,j}{F}_{\theta i,j}),\hfill \end{array}\end{array}$$(57)

where F_r i, j = F_r(r_i, θ_j + 1/2), F_θ i, j = F_θ(r_i + 1/2, θ_j); S_r i, j and S_θ i, j are the surface areas of the faces of each cell at constant r = r_i and θ = θ_j, respectively; and 𝒱_i, j is the volume of the cell. Therefore,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {(\mathbf{\nabla }·\mathit{F})}_{i,j}=& 3\frac{r_{i+1}^2{F}_{ri+1,j}-r_i^2{F}_{ri,j}}{r_{i+1}^3-r_i^3}\hfill \\ \hfill & +\frac{3(r_{i+1}^2-r_i^2)}{2(r_{i+1}^3-r_i^3)}\frac{sin{\theta }_{j+1}{F}_{\theta i,j+1}-sin{\theta }_j{F}_{\theta i,j}}{cos{\theta }_j-cos{\theta }_{j+1}}.\hfill \end{array}\end{array}$$(58)

Assuming boundary conditions v_n, r(r = R) = v_c, r(r = R) = 0, Eqs. (50) and (51) can be solved for the relative chemical perturbations, χ_n, and χ_c. Considering all grid cells, we have two sets of variables, {χ_n i, j}_{i = 1, j = 1}^{N_r − 1, N_θ − 1} and {χ_c i, j}_{i = 1, j = 1}^{N_r − 1, N_θ − 1}, giving a total of 2(N_r − 1)(N_θ − 1) unknowns. Eqs. (50) and (51) are evaluated at the center of each grid cell. Thus, each of them provides (N_r − 1)(N_θ − 1) equations, making a total of 2(N_r − 1)(N_θ − 1). We note that these equations only involve spatial derivatives of the relative chemical perturbations. Thus, they determine χ_n and χ_c up to additive constants. Fixing these constants requires extra conditions provided by Eqs. (14) and (15). Therefore, for the system of equations to be consistent, we must replace two redundant equations from the discretization of Eqs. (50) and (51). We choose to replace the equations corresponding to the grid cell (i, j) = (1, 1).

In the evaluation of each divergence, the boundary conditions on the velocities must be taken into account. For each cell touching the symmetry axis, we impose the fluxes through the corresponding face of the cell to be null (i.e., θ = 0 or θ = π). In the cells touching the center of the star, we impose fluxes through the face defined by r = 0 to be null, and at the faces defined by r = 1 we impose v_n, r = 0 and v_ad, r = 0, thus fluxes through those faces are null as well. In other words, radial fluxes at r = 1 of the cells touching the boundary are equal on both sides of equations (50) and (51).

Schematically, the linear system described before can be represented by a vector of unknown values of the relative chemical potential perturbations χ = {{χ_n i, j}_{i = 1, j = 1}^{N_r − 1, N_θ − 1},{χ_n i, j}_{i = 1, j = 1}^{N_r − 1, N_θ − 1}}; a vector b containing the right-hand side of Eqs. (50), (51), (14), and (15); and a sparse matrix, Λ, generated from the coefficients extracted from the discretization of the left-hand side of the previous equation:

$$\begin{array}{c}\hfill \mathrm{\Lambda }\left\{\begin{array}{c}{\{{\chi }_{ni,j}\}}_{i=1,j=1}^{N_r-1,N_{\theta }-1}\\ {\{{\chi }_{ci,j}\}}_{i=1,j=1}^{N_r-1,N_{\theta }-1}\end{array}\right\}=\mathit{b}\equiv \left\{\begin{array}{c}{\{\text{rhs}(50)\}}_{(i,j)\ne (1,1)}^{N_r-1,N_{\theta }-1}\\ {\{\text{rhs}(51)\}}_{(i,j)\ne (1,1)}^{N_r-1,N_{\theta }-1}\\ 0\\ 0\end{array}\right\}.\end{array}$$(59)

We note that all the dependence on the magnetic field is in b. Thus, Λ needs to be inverted only once in a time-evolving simulation. It is also worth noting that this approach has some limitations. Setting a very small value of ζ, is not only limited by the integration time. As ζ decreases, the matrix Λ becomes progressively more singular.

In summary, this solution consists of the following steps:

1. An initial seed for the magnetic field is chosen.

2. From the magnetic field, b is computed and used to solve for the chemical potential perturbations:

   $$\begin{array}{c}\hfill \mathit{\chi }={\mathrm{\Lambda }}^{-1}\mathit{b}.\end{array}$$(60)

3. The values of χ_n and χ_c are replaced into Eqs. (13) and (44), yielding the relevant velocities.

In the case we are also interested in following the time evolution:

1. The time derivative of the magnetic field from Eqs. (1) is computed and used to evolve the magnetic field one time step (see more details in Sect. 7).

2. This procedure is iterated, going back to step 1 but using the evolved magnetic field.

4.2. Numerical implementation with the Dedalus pseudo-spectral code

As an alternative numerical implementation of the FF approach, we solved Eqs. (50) and (51) with the normalization described in Table 1 using Dedalus. The truncation orders were set to n_r = 92 for the radial polynomials and n_θ = 65 for the angular polynomials. As explained in Sect. 3.2, Dedalus expands the radial and angular domains using Jacobi polynomials and spherical harmonics, respectively, leveraging their properties to solve the system in Eq. (60) efficiently.

We remark that all regularity conditions imposed by axial symmetry are naturally satisfied, as the polynomial basis used in the expansion is inherently regular at the axis and the origin. Additionally, the boundary conditions, Eqs. (53) and (54), are enforced in Dedalus using the τ-method (see, e.g., Boyd 2001).

5. Background neutron star model

We performed our analysis considering five different EoSs. One is the same toy model described in Castillo et al. (2020). There, neutrons and protons are treated as a non-relativistic Fermi gas, and electrons are treated as an extremely relativistic Fermi gas. Both neutrons and protons are considered to have the same mass m. Therefore, the chemical potentials are related to the particle densities by

$$\begin{array}{cc}\hfill {\mu }_n(r)& =m{c}^2+\frac{p_{\mathit{Fn}}{(r)}^2}{2m},\hfill \end{array}$$(61)

$$\begin{array}{cc}\hfill {\mu }_c(r)& =m{c}^2+\frac{p_{\mathit{Fc}}{(r)}^2}{2m}+p_{\mathit{Fc}}(r)c,\hfill \end{array}$$(62)

where p_Fi(r) = ℏ[3π²n_i(r)]^1/3  (i = n, c) are the Fermi momenta of neutrons and charged particles, respectively. Thus, the density profiles can be written as

$$\begin{array}{cc}\hfill n_n(r)& ={[\frac{2mc}{ħ{(3{\pi }^2)}^{1/3}}{n}_c{(r)}^{1/3}+n_c{(r)}^{2/3}]}^{3/2},\hfill \end{array}$$(63)

$$\begin{array}{cc}\hfill n_c(r)& =n_0{\text{sinc}}^6(0.7\pi r/R),\hfill \end{array}$$(64)

where sinc(x) = sin(x)/x. Equation (63) comes from the chemical equilibrium condition, μ_n(r) = μ_c(r), and Eq. (64) is an analytic approximation to the radial profile obtained by numerically integrating the (non-relativistic) stellar structure equations.

The other EoSs used are HHJ (Heiselberg & Hjorth-Jensen 1999), BSk24, BSk25, and BSk26 (Goriely et al. 2013), which have been adapted for npe matter. For all realistic models, the central density is adjusted so that the mass of the star is 1.4M_⊙. The basic parameters of the stellar models are summarized in Table 2. Unlike the toy model, these realistic EoSs consider strong interactions by means of a coefficient K ≠ 0 (see Sect. 2). A comparison between n_n, n_c, and n_c/n_n as functions of r for the different models can be seen in Figs. 1(a) and (b). The collision coefficient between charged particles and neutrons, γ_cn, for the different models is computed from J_pn(r, T) = n_c(r)n_n(r)γ_cn(r, T) given by Yakovlev & Shalybkov (1991), where T is the temperature.

Fig. 1. Comparison between the five different EoSs studied. (a) nn and nc, in units of nc(r = 0); (b) the ratio nc/nn; and (c) ℓc/R, where ℓc ≡ −[dln(nc/nn)/dr]−1 (see Sect. 7.2) and R is the core radius. All these variables are functions of the normalized radial coordinate, r/R.

6. Analytical magnetic field models

In our simulations, initial conditions are given by analytical expressions for α and β as functions of r and θ. These functions are normalized so that the rms value of the magnetic field in the core is 1:

$$\begin{array}{c}\hfill ⟨\mathit{B}⟩\equiv {\left(\frac{1}{{\mathcal{V}}_c}{\displaystyle {\int }_{{\mathcal{V}}_c}}{\mathit{B}}^{2}d\mathcal{V}\right)}^{1/2}=1.\end{array}$$(65)

We focus on the following models:

Model I: Arguably, the simplest expressions for a purely poloidal magnetic configuration are of the form

$$\begin{array}{c}\hfill \alpha (r,\theta )=f(r){sin}^2\theta ,\end{array}$$(66)

where f(r) is a polynomial on r. In particular, the following has been widely used in the literature (Akgün et al. 2013; Armaza et al. 2015; Passamonti et al. 2017; Ofengeim & Gusakov 2018; Castillo et al. 2020):

$$\begin{array}{c}\hfill {\alpha }_{\text{I}}(r,\theta )=\frac{35}{8}\sqrt{\frac{11}{118}}(r^2-\frac{6}{5}{r}^4+\frac{3}{7}{r}^6){sin}^2\theta .\end{array}$$(67)

The numerical coefficients on the polynomial come from imposing continuity of B on the surface of the core with a current-free dipole solution outside [f(r)∝1/r], continuity of the toroidal component of the current density J_ϕ (although this is not a physical requirement), and the normalization of the magnetic field (in addition to a regularity condition at the axis, which does not allow for terms proportional tor⁰ and r¹).

Model II, a model with null radial velocities at the crust-core interface: Following Ofengeim & Gusakov (2018), we can construct poloidal magnetic field configurations for which the radial components of both v_n and v_ad vanish at the crust-core interface. First, we take α as a multipolar expansion ofthe form

$$\begin{array}{c}\hfill \alpha (r,\theta )={\displaystyle \sum _{l\ge 1}}{C}_l{\alpha }_l(r,\theta ),\end{array}$$(68)

where

$$\begin{array}{c}\hfill {\alpha }_l(r,\theta )=f_l(r)P_l^1(cos\theta )sin\theta .\end{array}$$(69)

Here, C_l are constants that determine the relative weights of each multipole, f_l(r) = ∑_ia_l, irⁱ are polynomials of r, and P_l¹ are the associated Legendre functions. The potentials α_l are normalized so their generated magnetic field B_l ≡ ∇α_l × ∇ϕ satisfies ⟨B_l⟩ = 1, implying ∑_lC_l² = 1. For this model, we only considered a dipole; however, equations are written for a more general case, as they also apply to model III. The order of the polynomial(s) depends on the number of conditions to be imposed. We enforce the same three conditions as for model I, namely

$$\begin{array}{cc}& f_l^{\prime }(1)+l{f}_l(1)=0,\hfill \end{array}$$(70)

$$\begin{array}{cc}& f_l^{″}(1)-l(l+1)f_l(1)=0,\hfill \end{array}$$(71)

$$\begin{array}{cc}& ⟨{\mathit{B}}_l⟩=1.\hfill \end{array}$$(72)

These are, respectively, continuity of B at the crust-core boundary considering a current-free field everywhere outside the core; continuity of J_ϕ, which is null outside; and a normalization.

Additional conditions follow from setting v_n, r(1, θ) = v_ad, r(1, θ) = 0. If we expand v_ad, r(1, θ), using Eq. (28), into Legendre polynomials, all Legendre components should be null, which implies ${\widehat{\mathrm{P}}}_l\{\partial X_n/\partial r{|}_{r=1}\}=0$, l = 1, 2, …. Thus, the first extra condition is

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^{\pi }}{\frac{\partial X_n}{\partial r}|}_{r=1}{P}_l(cos\theta )sin\theta =0,l=1,2,..\end{array}$$(73)

Analogously, Eq. (30) implies ${\widehat{\mathrm{P}}}_l\{\mathbf{\nabla }·[(\mu /{\gamma }_{\mathit{cn}})\mathbf{\nabla }{X}_n]{|}_{r=1}\}=0$, l = 1, 2, …, which, together with Eq. (73), yields the second extra condition

$$\begin{array}{c}\hfill {\displaystyle {\int }_0^{\pi }}{[\frac{{\partial }^2{X}_n}{\partial r^2}-\frac{l(l+1)}{r^2}{X}_n]}_{r=1}{P}_l(cos\theta )sin\theta =0,l=1,2,...\end{array}$$(74)

To apply these two conditions, we first need to obtain X_n in terms of the set of coefficients {a_l, i} in f_l(r). This is done first by computing the magnetic force in terms of the a_l, i coefficients, then using (20), (21), and (22). The number of non-zero Legendre components in X_n, and consequently in Eqs. (73) and (74), increases nonlinearly with the number of non-zero Legendre components in B (as X_n is nonlinear in B). For a dipolar magnetic field configuration (corresponding to l = 1), the only non-zero Legendre component in X_n is l = 2, which imposes non-linear relations for the coefficients {a_{l = 1, i}} in Eqs. (73) and (74). We have a total of five conditions; thus, we used a polynomial of five coefficients of the form

$$\begin{array}{c}\hfill f_1(r)=a_{1,2}{r}^2+a_{1,4}{r}^4+a_{1,6}{r}^6+a_{1,8}{r}^8+a_{1,10}{r}^{10}.\end{array}$$(75)

Solving for the coefficients, using the conditions previously described, yields eight solutions, four of which, evaluated with the HHJ EoS, are summarized in Table 3 and displayed in Fig. 2. The other four are the same, just with a global minus sign. We focus on “Solution 1” (hereafter model II, α_II); as in the others, the magnetic flux is fully enclosed inside the core.

Fig. 2. Radial functions f1(r) for model II. The coefficients for each solution are given in Table 3.

Model III, a non-equatorially symmetric model with null radial velocities at the crust-core interface: Using the conditions from the previous model, we can build a non-equatorially symmetric poloidal magnetic configuration with null radial velocities at the crust-core interface. For this, we consider a combination of a dipolar and quadrupolar magnetic field configuration of the form ${\alpha }_{\text{III}}(r,\theta )=\sqrt{0.5}{\alpha }_1(r,\theta )+\sqrt{0.5}{\alpha }_2(r,\theta )$, where

$$\begin{array}{cc}\hfill {\alpha }_1(r,\theta )& ={\displaystyle \sum _{i=1}^6}{a}_{1,2i}{r}^{2i}{sin}^2\theta ,\hfill \end{array}$$(76)

$$\begin{array}{cc}\hfill {\alpha }_2(r,\theta )& ={\displaystyle \sum _{i=1}^6}{a}_{2,2i+1}{r}^{2i+1}{sin}^2\theta cos\theta .\hfill \end{array}$$(77)

The number of a_l, i coefficients considered comes from the number of conditions needed to fix the configuration. For each of the multipoles (l = 1 and l = 2), we have to satisfy the conditions from Eqs. (70) and (72); a total of four. We note that we relaxed condition (71), as it is not strictly needed. Also, α_III yields a function, X_n(r, θ), that includes non-zero multipoles l = 1, 2, 3, and 4. Thus, Eqs. (73) and (74) yield eight equations. So, we have set enough coefficients to satisfy the 12 conditions needed for α_III. There is a very large number of solutions for this system. The coefficients for the solution we use in this work are presented in Table 4.

Model IV, a non-equatorially symmetric model with a toroidal component: We use the same non-equatorially symmetric magnetic field introduced in Castillo et al. (2020), with poloidal and toroidal components produced by the potentials

$$\begin{array}{cc}\hfill {\alpha }_{\text{IV}}(r,\theta )& =\sqrt{0.18}{\alpha }_{\text{I}}(r,\theta )+\sqrt{0.42}{\alpha }_{\text{aux}}(r,\theta ),\hfill \end{array}$$(78)

$$\begin{array}{cc}\hfill {\beta }_{\text{IV}}(r,\theta )& =\sqrt{0.4}{\beta }_{\text{aux}}(r,\theta ),\hfill \end{array}$$(79)

where

$$\begin{array}{cc}\hfill {\alpha }_{\text{aux}}(r,\theta )& =3.718(r^3-\frac{10}{7}{r}^5+\frac{5}{9}{r}^7){sin}^2\theta cos\theta ,\hfill \end{array}$$(80)

$$\begin{array}{cc}\hfill {\beta }_{\text{aux}}(r,\theta )& =112.546r^5{(1-r)}^2{sin}^2\theta sin(\theta -\pi /5).\hfill \end{array}$$(81)

The quadrupolar poloidal potential, α_aux, (which satisfies the same conditions as α_I), and the toroidal potential, β_aux, are both normalized so ⟨B_Pol^aux⟩ = 1 and ⟨B_Tor^aux⟩ = 1, respectively³; so α_IV and β_IV will produce a magnetic configuration in which the toroidal magnetic field has 40% of the internal magnetic energy.

The magnetic field models I, II, III, and IV are shown in the first row of Fig. 3.

Fig. 3. Magnetic field configurations, relative chemical potential perturbations, and velocity fields produced by the models I, II, III, and IV. For B, lines represent the poloidal magnetic field lines labeled by the corresponding value of α. For models I–III, color represents the poloidal potential α, while in model IV, it represents the toroidal magnetic field Bϕ. For χn and χc, color represents their local value. For vn and vad, lines represent the poloidal component and color its magnitude; these are computed through the explicit semi-analytic approach described in Sect. 3.1 using the background NS model with the HHJ EoS.

7. Time evolution

The time-evolving equations for the magnetic potentials can be obtained from Eq. (1) together with Eq. (16)

$$\begin{array}{cc}\hfill \frac{\partial \alpha }{\partial t}& =rsin\theta [({\mathit{v}}_n+{\mathit{v}}_{\mathit{ad}})\times \mathit{B}]·\widehat{\varphi },\hfill \end{array}$$(82)

$$\begin{array}{cc}\hfill \frac{\partial \beta }{\partial t}& =r^2{sin}^2\theta \mathbf{\nabla }·\left\{\frac{[({\mathit{v}}_n+{\mathit{v}}_{\mathit{ad}})\times \mathit{B}]\times \widehat{\varphi }}{rsin\theta }\right\}.\hfill \end{array}$$(83)

The evolution in the core and crust of the star are coupled. For simplicity, in this work we assume that the crust acts as a vacuum or perfect resistor (∇ × B = 0). Therefore, the magnetic field outside of the core can be written as B = ∇φ, where $\phi (r,\theta )={\displaystyle \sum _{l=1}^{\infty }}(b_l/r^{l+1})P_l(cos\theta )$ and b_l are coefficients. The magnetic field in the crust and outside of the star is entirely determined by the radial component of the magnetic field at the crust-core interface, where continuity implies $b_l=-R^{l+2}{\widehat{\mathrm{P}}}_l\{B_r(R,\theta )\}/(l+1)$ (see Marchant et al. 2011). Of course, in practice, we are only considering a finite number of external multipoles, N_Exp, satisfying the criterion N_Exp ≥ N_θ/4. It should be noted that the assumption of a current-free crust is properly justified only if currents therein dissipate much faster than in the core. The latter, however, is very uncertain (Jones 2001; Cumming et al. 2004), and simulations coupling the evolution of the core and crust may be more appropriate (see discussion in Sect. 9, numeral 5).

In our time-evolving simulations, the numerical computation of the time derivative of the toroidal potential β is performed conservatively using a finite-volume scheme, and the time derivative of the poloidal potential α is calculated using finite differences. The decomposition of the magnetic field in terms of α and β ensures that the condition ∇ ⋅ B = 0 is maintained at all times with machine precision. The time evolution is computed using the midpoint method (a second-order Runge–Kutta method), which allows the system to be evolved with second-order accuracy in time. For more details on the numerical code, please refer to Castillo et al. (2017).

7.1. Hydromagnetic quasi-equilibrium

As discussed in Sect. 1, after a few Alfvén times, the magnetic field will evolve toward a hydromagnetic quasi-equilibrium state in which, under axial symmetry, the toroidal magnetic force must be null, since the fluid forces, f_i = −n_iμ∇χ_i (i = n, c), which are purely poloidal, cannot balance it. This imposes a strong condition on the magnetic field; Eq. (16) implies f_B, ϕ ∝ |∇α × ∇β| = 0; which can only be achieved if there is a relation between the magnetic potentials (e.g., β = β(α)). Since in our model there is no toroidal field outside of the core, the latter implies that the toroidal field must be confined to the regions in which the poloidal magnetic field lines close inside of the core, forming a “twisted torus” configuration (Chandrasekhar & Prendergast 1956; Braithwaite & Spruit 2004). It also implies that the magnetic force takes the form

$$\begin{array}{c}\hfill {\mathit{f}}_B=-\frac{{\mathrm{\Delta }}^{\ast }\alpha +\beta d\beta /d\alpha }{4\pi r^2{sin}^2\theta }\mathbf{\nabla }\alpha ,\end{array}$$(84)

where

$$\begin{array}{c}\hfill {\mathrm{\Delta }}^{\ast }\equiv \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}$$(85)

is the so-called GS operator.

As time passes, this quasi-equilibrium is slowly eroded as the magnetic field lines, which are frozen into the charged particle fluid, keep moving relative to the neutron fluid with velocity v_ad. Eventually, a final barotropic equilibrium state is reached on a timescale ∼t_B (see Sect. 7.2), in which the magnetic force is mainly balanced by the pressure gradient and gravitational force of the charged particle fluid, while neutrons stay in diffusive equilibrium,

$$\begin{array}{cc}& {\mathit{f}}_B-n_c\mu \mathbf{\nabla }{\chi }_c=0,\hfill \end{array}$$(86)

$$\begin{array}{cc}& \mathbf{\nabla }{\chi }_n=0.\hfill \end{array}$$(87)

The latter formula indeed means diffusive equilibrium, since it is equivalent to the condition ∇δμ_n + (δμ_n/c²)∇Φ = 0 (see Eq. (8) and the definition of the function χ_n).

Eqs. (84) and (86) imply that, in this kind of equilibrium, not only β, but also χ_c, must be a function of α, and these variables must satisfy the so-called GS equation

$$\begin{array}{c}\hfill {\mathrm{\Delta }}^{\ast }\alpha +\beta \frac{d\beta }{d\alpha }+4\pi r^2{sin}^2\theta n_c\mu \frac{d{\chi }_c}{d\alpha }=0.\end{array}$$(88)

Solutions to this equation are called “GS equilibria” (Grad & Rubin 1958; Shafranov 1966) and have been widely studied in the literature (e.g., Armaza et al. 2015 and references therein). We note that in this equilibrium χ_n must be uniform, but χ_c is non-uniform. Together with Eq. (4), this implies that, unlike what one might have expected, δn_n cannot be uniform.

7.2. Time-scale estimates

From Eq. (1), we can estimate the timescale for magnetic evolution in the core of the star as

$$\begin{array}{c}\hfill t_B\sim \frac{{\ell }_B}{|{\mathit{v}}_c|}=\frac{{\ell }_B}{|{\mathit{v}}_n+{\mathit{v}}_{\mathit{ad}}|},\end{array}$$(89)

where ℓ_B is the length scale over which the magnetic field varies. The usual assumptions found in the literature (Goldreich & Reisenegger 1992; Thompson & Duncan 1996; Castillo et al. 2017; Passamonti et al. 2017) are that the neutron velocity can be neglected and that the magnitude of the poloidal component of the ambipolar velocity can be estimated from Eq. (13) as v_{ad, Pol} ∼ f_B/(γ_cnn_cn_n) (v_ad, ϕ = 0 in the quasi-stationary state). From Eq. (11), we see that the magnitude of the magnetic force can be estimated as f_B ∼ B²/(4πℓ_B). Under these assumptions, we obtain the usual estimate for magnetic evolution under ambipolar diffusion as t_B ∼ t_ad, where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill t_{\mathit{ad}}\equiv \frac{{\ell }_B}{|{\mathit{v}}_{ad\text{, Pol}}|}& \sim \frac{4\pi {\gamma }_{\mathit{cn}}{n}_c{n}_n{\ell }_B^2}{B^2}\hfill \\ \hfill & \sim 8\times {10}^5{\left(\frac{{10}^{14}\text{G}}{B}\right)}^2{\left(\frac{T}{{10}^8\text{K}}\right)}^2{\left(\frac{{\ell }_B}{1\text{km}}\right)}^2\text{yr},\hfill \end{array}\end{array}$$(90)

for BSk24. This derivation, however, has a few shortcomings. First, from Eqs. (12) and (13), we see that a better estimate of the magnitude of the poloidal component of the ambipolar velocity is v_ad ∼ f_n/(γ_cnn_cn_n); however, as demonstrated by Moraga et al. (2024), an arbitrary magnetic field can induce fluid forces larger than the Lorentz force, i.e., f_c ∼ f_n ∼ (ℓ_c/ℓ_B)f_B, where we introduced a length scale

$$\begin{array}{c}\hfill {\ell }_c\equiv -{[\frac{d}{\mathit{dr}}ln(n_c/n_n)]}^{-1},\end{array}$$(91)

which is a measure of the stable stratification due to the composition gradient in the NS core [more strongly stratified for smaller ℓ_c, see Fig. 1(c)]. Thus,

$$\begin{array}{c}\hfill v_{ad\text{, Pol}}\sim \frac{{\ell }_c}{{\ell }_B}\frac{f_B}{{\gamma }_{\mathit{cn}}{n}_c{n}_n}.\end{array}$$(92)

More importantly, the implicit assumption of fixed neutrons is not justified. As shown in Ofengeim & Gusakov (2018) and numerically checked by Castillo et al. (2020), the neutron velocity can be much larger than the ambipolar diffusion velocity. Its magnitude can be estimated from Eqs. (9) and (30) as

$$\begin{array}{c}\hfill v_{n\text{, Pol}}\sim \frac{{\ell }_c}{{\ell }_B}\frac{f_n}{{\gamma }_{\mathit{cn}}{n}_c{n}_n}\sim \frac{{\ell }_c}{{\ell }_B}{v}_{ad\text{, Pol}}.\end{array}$$(93)