Magneto-thermal evolution in the cores of adolescent neutron stars: The Grad-Shafranov equilibrium is never reached in the ‘strong-coupling’ regime

At the high temperatures present inside recently formed neutron stars ( 𝑇 ≳ 5 × 10 8 K), the particles in their cores are in the "strong-coupling" regime, in which collisional forces make them behave as a single, stably stratified, and thus non-barotropic fluid. In this regime, axially symmetric hydromagnetic quasi-equilibrium states are possible, which are only constrained to have a vanishing azimuthal Lorentz force. In such equilibria, the particle species are not in chemical ( 𝛽 ) equilibrium, so 𝛽 decays (Urca reactions) tend to restore the chemical equilibrium, inducing fluid motions that change the magnetic field configuration. If the stars remained hot for a sufficiently long time, this evolution would eventually lead to a chemical equilibrium state, in which the fluid is barotropic and the magnetic field, if axially-symmetric, satisfies the non-linear Grad-Shafranov equation. In this work, we present a numerical scheme that decouples the magnetic and thermal evolution, enabling to efficiently perform, for the first time, long-term magneto-thermal simulations in this regime for different magnetic field strengths and geometries. Our results demonstrate that, even for magnetar-strength fields ≳ 10 16 G, the feedback from the magnetic evolution on the thermal evolution is negligible. Thus, as the core passively cools, the Urca reactions quickly become inefficient at restoring chemical equilibrium, so the magnetic field evolves very little, and the Grad-Shafranov equilibrium is not attained in this regime. Therefore, any substantial evolution of the core magnetic field must occur later, in the cooler "weak-coupling" regime ( 𝑇 ≲ 5 × 10 8 K), in which Urca reactions are effectively frozen and ambipolar diffusion becomes relevant.


INTRODUCTION
Neutron stars (NSs) are grouped into several classes according to their very different surface magnetic field strengths, rotational periods, and the presence or absence of accretion from a binary companion.These properties vary noticeably; for example, millisecond pulsars (MSPs) have surface magnetic fields that typically range around ∼ 10 8−9 G, whilst in magnetars these values can reach ∼ 10 14−15 G at their surface (Kaspi & Beloborodov 2017), with possibly even stronger internal fields.Therefore, achieving a detailed comprehension of the different processes that lead to the evolution of the magnetic field of NSs is a key step to understand their observed phenomenology and how some of the different NS classes could be related to each other.
The problem of how the NS magnetic field evolves has been deeply studied in the literature throughout the last decades, from the seminal paper by Goldreich & Reisenegger (1992) to a substantial amount of recent work that approached the problem either theoretically or numerically, given its complicated and non-linear characteristics (see Pons & Viganò 2019 for a recent review).Among all the mechanisms that promote the magnetic evolution, by far the best understood are ★ E-mail: nicolas.moraga@ug.uchile.clthose occurring in the NS crust, where ions have very restricted mobility (unless the magnetic field is strong enough to break the crust), so the evolution depends only on the motion of the electrons.Therefore, the long-term mechanisms that control the evolution of the field in this region are Ohmic diffusion, i. e., current dissipation by electric resistivity; and Hall drift, which corresponds to advection of the magnetic field lines by the electron fluid motion.Though Hall drift by itself does not dissipate magnetic energy, it changes the magnetic field configuration, creating small-scale structures that dissipate more quickly, and, together with Ohmic dissipation, it may lead to a special, quasi-steady magnetic field configuration called a "Hall attractor" (Cumming et al. 2004;Gourgouliatos & Cumming 2014).
In the NS core, the physics is more complex and not yet well understood.The NS core is a mixture of neutrons (mostly), protons, and electrons, joined by muons and other, more exotic species at increasing densities.Also, it is generally accepted that neutrons and protons can become superfluid and superconducting respectively, at temperatures  ∼ 10 8 − 10 10 K (Migdal 1959;Page et al. 2013), leading to a complex dynamics (Glampedakis et al. 2011;Gusakov & Dommes 2016;Gusakov et al. 2020;Dommes & Gusakov 2021).In the present paper, we will restrict ourselves to the case of a "normal" core, i. e., non-superconducting and non-superfluid, which is likely to be realistic for the high temperatures and strong magnetic fields that our work will mostly focus on.
After the supernova explosion, a proto-NS is born in an extremely hot, liquid state opaque to neutrinos.However, within a minute, the star becomes neutrino-transparent and the young NS is formed (Burrows & Lattimer 1986;Keil & Janka 1995;Pons et al. 1999;Lander et al. 2021).At nearly the same time, the magnetic field reaches an equilibrium configuration that is likely to fill the entire volume of the NS, as in the simulations of Braithwaite & Spruit (2004) and Becerra et al. (2022).This configuration is stabilized by the composition gradient of matter (radially decreasing proton-to-neutron ratio) within the NS core, which causes radial buoyancy forces that oppose convective motions (Pethick 1992;Reisenegger & Goldreich 1992;Goldreich & Reisenegger 1992).On the other hand, since the electric conductivity of the NS core is very high, the magnetic field is effectively frozen into the charged particles.Thus, any process that promotes the magnetic evolution has to move the charged particles in such a way that the stable stratification is overcome.This can be achieved in different ways in the so-called "strong-coupling" and "weak-coupling" regimes, which are effective in different temperature ranges (Goldreich & Reisenegger 1992;Hoyos et al. 2008;Reisenegger 2009;Gusakov et al. 2017;Castillo et al. 2020): (i) At high temperatures, all particle species in the core are strongly coupled to each other due to frequent collisions.As a result, the charged particles and neutrons move together as a single fluid that is strongly constrained by buoyancy forces as a consequence of stable stratification.However, by the effect of -decays (so-called "Urca-reactions", e. g., Shapiro & Teukolsky 1983), the fluid elements can gradually adjust their chemical composition, so that the stable stratification is overcome, and the fluid can transport the magnetic flux.Therefore, the dynamics of this "strong-coupling" regime is characterized by bulk motions where matter behaves as a single, stably stratified, non-barotropic fluid (i.e., the pressure depends on the non-uniform chemical composition in addition to density), where the interplay between the magnetic field dynamics and Urca reactions determines the time scale over which the magnetic field evolves.
(ii) As the star cools, the NS core transits into the "weak-coupling" regime, in which the collisional coupling has decreased enough to make ambipolar diffusion (the relative motion of charged particles with respect to neutrons) possible, while the rates of the Urca reactions quickly drop, preventing further conversions between the particle species.In this regime, the finite relative velocity of neutrons and charged particles, although smaller than the bulk flow velocity (Ofengeim & Gusakov 2018;Castillo et al. 2020), is crucial to allow the magnetic field to evolve and determines the rate at which magnetic energy is dissipated through collisions between neutrons and charged particles.This process heats the NS core and has been invoked to explain the activity of magnetars due to its strong dependence on the magnetic field intensity (Thompson & Duncan 1995, 1996;Beloborodov & Li 2016;Tsuruta et al. 2023).It has also been suggested to explain the low magnetic fields of MSPs, as the collisional coupling decreases at low temperatures, so the time scale of ambipolar diffusion may become short enough to cause substantial magnetic field decay in old NSs before their spin-up by accretion (Cruces et al. 2019).However, it is important to note that this scenario ignores the effects of superfluidity and superconductivity, which may have a significant impact on the dissipation of magnetic energy in the core of old NSs (see e. g.Kantor & Gusakov 2018;Gusakov et al. 2020).
We note that this kind of dynamics (especially for the strongcoupling regime) was not correctly described by Goldreich & Reisenegger (1992) and Thompson & Duncan (1996), who treated the neutrons as a fixed background and thus required the motion of the charged particles to overcome the friction force due to protonneutron collisions, which is very strong at high temperatures, but which does not act if all particles move together.A similar approach was also adopted by Passamonti et al. (2017); Viganò et al. (2021) and later by Igoshev & Hollerbach (2023).
The formal model with moving neutrons has been discussed by several authors, for instance, Shalybkov & Urpin (1995); Reisenegger et al. (2005); Reisenegger (2009).In these studies, the velocity of the neutrons was not explicitly calculated.The original method for calculating all velocities of different particle species for a given stellar magnetic field configuration was proposed in Gusakov et al. (2017).That work demonstrated that a magnetic field inevitably leads to the emergence of large-scale flows of stellar matter, which influence the evolution of the magnetic field.Further, in the study by Ofengeim & Gusakov (2018), the velocity of these flows was explicitly calculated for a range of magnetic field models.It was shown that the velocity of matter as a whole can significantly exceed the relative (diffusion) particle velocities, which were calculated in previous studies.A conclusion was drawn about the substantial acceleration of magnetic evolution when accounting for this (previously overlooked) effect.These findings were later confirmed in the detailed self-consistent numerical simulations of the magnetic field evolution by Castillo et al. (2020) (see below for details).Castillo et al. (2017) performed the first simulations that evolved the magnetic field and the small density perturbations it induces on the charged particles inside a NS core in axial symmetry (2 dimensions).This model did not include Urca reactions and assumed constant friction coefficients (thus implicitly constant temperature) and motionless neutrons, which represents a simplified version of the weak-coupling regime.Later, Castillo et al. (2020) included the motion of the neutrons as well as the different density profiles of neutrons and charged particles (with a toy-model equation of state), thus explicitly accounting for the stable stratification, while still focusing on the weak-coupling regime at constant temperature.
Despite the improvements that have gradually been implemented in these papers, an important limitation persists within their numerical framework as the quotient between the time scales of interest is proportional to the ratio of magnetic to fluid pressure, which is ≲ 10 −6 for realistic magnetic fields.Thus, unrealistically large magnetic field strengths are needed to run the simulations for a long enough time.In addition, the crucial influence of the thermal evolution on the dynamics of the NS core has not been considered in the contexts of the strong-coupling regime, or in the weak-coupling regime with independent velocity fields for neutrons and charged particles.This work can be considered as the continuation of the previous papers, addressing these issue.
The aim of this work is to report the first long-term simulations of the magneto-thermal evolution of a NS core in the strong-coupling regime.We present a new numerical scheme where, by neglecting the time derivatives in the continuity equations, the quotient between the time scales of interest becomes independent of the magnetic field strength, allowing to perform simulations with realistic parameters and thus speeding up the simulation code considerably.Also, we apply a strategy that allows us to separate the magnetic and thermal evolution in the equations, in which we first run the evolution of the magnetic field at constant temperature and then introduce the thermal evolution through a change in the time variable (see § 4 for more details).
This paper is structured as follows: In § 2, we present the main equations and assumptions governing the physical model, which is to be numerically solved in axial symmetry.We discuss the time scales associated with the problem and describe the application of the novel numerical scheme to the strong-coupling regime.In § 3, we present the outcomes of our simulations initially under constant temperature conditions.Subsequently, in § 4, we describe the strategy to include the thermal evolution and present the results of the magneto-thermal evolution.Lastly, in § 5, we summarize our results and outline the main conclusions.

General equations
Following the model of Castillo et al. (2020), we consider a spherically symmetric, non-rotating, non-magnetic background NS core composed of neutrons (), protons (), and electrons () in chemical and hydrostatic equilibrium, i. e.
where the radial functions () and   () ( = , ) are the gravitational potential and the chemical potentials, respectively.The subscript  denotes charged particles, i. e.,   ≡   +   .The presence of the magnetic field perturbs this equilibrium.However, since the ratio between the magnetic pressure and degeneracy pressure  in NS interiors is  2 /8 ≲ 10 −6 (Reisenegger 2009), the perturbations induced by the magnetic field are expected to be of the same order.Thus, the chemical potentials and particle densities can be split into a time-independent, spherically symmetric background and a much smaller, time-dependent perturbation: where the condition of charge neutrality allows to write   (, ) =   (, ) ≡   (, ).The number density perturbations are related to the chemical potential perturbations by and satisfy baryon number conservation: where V is the core volume.For simplicity, we consider an axially symmetric magnetic field, which can be decomposed as where the scalar potentials (, , ) and (, , ) generate the poloidal (meridional) and toroidal (azimuthal) magnetic field components, respectively.Here,  is the radial coordinate, and  and  are the polar and azimuthal angles, respectively, so ∇ = φ/( sin ), where φ is the unit vector in the  direction.The functions (, , ) and (, , ) are known as the poloidal flux and poloidal current functions, respectively, because 2(, , ) is the magnetic flux and (, , )/2 is the electric current enclosed by an azimuthal circle at given  and .The curves  = const are the poloidal magnetic field lines.
In order to study the magnetic field evolution, we need to integrate the Faraday induction equation where   =   +   is the charged-particle velocity, obtained as the sum of the neutron velocity   and the ambipolar velocity   (relative motion of charged particles with respect to neutrons).Using the poloidal-toroidal decomposition (equation 7), equation ( 8) can be written as the following coupled equations for  and : Therefore, in order to evolve the magnetic field, we need to compute   , which is obtained from the sum of and The first term in the parenthesis on the right-hand side of equation ( 11) is the Lorentz force, whose poloidal and toroidal components can be written as and respectively, where is the "Grad-Shafranov (GS) operator".The other two terms include the forces due to degeneracy pressure and gravity on each of the particle species (neutrons and charged particles) and can be written as with   ≡   / ( = , ).The latter forces are purely poloidal since we are assuming axial symmetry.Throughout this work, we use the Cowling approximation, i. e., we do not consider the perturbations of the gravitational potential.The parameter  in equation ( 11) is introduced through the artificial friction force   ≡ −     acting on the neutrons in order to follow the long-term quasi-stationary evolution (Hoyos et al. 2008;Castillo et al. 2020).This force replaces and mimics the short-term dynamics associated to the inertial terms in the Euler equations, which evolve the NS core towards a hydromagnetic quasi-equilibrium state, in which all the physical (i.e., non-fictitious) forces acting on the fluid elements are close to balancing each other:   +   +   ≈ 0. This approximation is justified because the fluid in the NS core is expected to reach the quasi-equilibrium state with the magnetic field within a few Alfvén crossing times ∼ 1 (10 14 G/) s, and we are interested in the NS evolution on much longer time scales.
In equation ( 12),   is a rate coefficient that parametrizes the frictional drag forces due to collisions between charged particles and neutrons, and can be written as (Yakovlev & Shalybkov 1991) where  is the mass density and  nuc = 0.16 fm −3 and  nuc = 2.8 × 10 14 g cm −3 are the nuclear saturation number density and mass density, respectively.Here,  is the core temperature.
The velocity fields   and   are also further constrained to satisfy the continuity equations at all times, namely, where the time derivatives of the number density perturbations,   ( = , ), have been neglected since a simple estimate gives Here, ΔΓ is the net conversion rate per unit volume of charged particles to neutrons by nonequilibrium Urca processes (direct or modified), which in general has a non-linear dependence on the variable Δ/   (Reisenegger 1995), where Δ ≡   −   is the chemical imbalance and   the Boltzmann constant.However, we will see that in the present problem the departure from chemical equilibrium remains modest, |Δ| ≲   , in the so-called "sub-thermal" regime or approximation (Haensel et al. 2002), in which ΔΓ is proportional to the chemical imbalance where  is a function that depends strongly on temperature and on the specific Urca processes involved.
For this EoS, the direct Urca mechanism is allowed in the core only for stellar masses  > 1.83  ⊙ .Hence, we focus exclusively on the modified Urca process, for which  is given by (Gusakov et al. 2017)

"Strong-coupling" regime
Since the coefficients  and   depend on temperature, it is necessary to also follow the evolution of the temperature in order to obtain the evolution of the magnetic field.At any given time, and thus for a given magnetic field configuration and temperature, equations ( 11)-( 12) and ( 19)-( 20) would need to be used to obtain the chemical potential perturbations and the velocities, which would then be used to take a step in time to calculate the new temperature and magnetic field configuration.
In order to simplify this process, we note that Urca reactions and ambipolar diffusion have strong (and opposite) dependencies on the temperature, namely  ∝  6 and   ∝  −1  ∝  −2 (see equations 18 and 22).Due to this strong dependence, the moment when both mechanisms operate with equal efficiency implicitly defines a temperature,    , that marks the transition from the strong-coupling to the weak-coupling regime, and whose derivation will be given in § 2.4.2.Here, we focus on the strong-coupling regime, which corresponds to high temperatures  >    and can be described by setting the ambipolar velocity to zero (  → ∞).Thus,   =   ≡ , and the full set of equations to solve in axial symmetry becomes Equations ( 25)-( 28) must be solved first to obtain the chemical potential perturbations and the velocity for a given magnetic field configuration, and then the obtained velocity must be used to evolve the magnetic field through equations ( 23) and ( 24).For the first step, we replace the expression for the velocity (equation 27) into the continuity equations ( 25)-( 26), obtaining where Δ ≡   −   .For a given magnetic field configuration (which determines the right-hand side of both equations) and with the additional constraint of baryon number conservation (equation 6), these can be solved numerically for   and   by a finite-element method.
Then, it is easy to compute the velocity  through equation ( 27) and forward-step the magnetic field components through equations ( 23) and ( 24), as illustrated in Fig. 1.

Boundary conditions
In the work of Castillo et al. (2017Castillo et al. ( , 2020)), it was assumed that the currents in the crust decay much faster than the typical evolution time scales in the core, so the crust was treated as a vacuum (a perfect resistor) whose magnetic field at any time is fully determined by the field in the core.This is unlikely to be satisfied in the strong-coupling regime because of the short time scales involved and because the crust is still partly liquid and only progressively freezes during this period of time (e.g., Aguilera et al. 2008).Nevertheless, and for simplicity, we proceed in the same way as Castillo et al. (2017Castillo et al. ( , 2020)), considering the external magnetic configuration as a current-free, poloidal magnetic field that is computed at all time-steps as a multipolar expansion and imposing continuity of the poloidal component at the crust-core interface,  =    .We note that this assumption should allow a faster evolution than the opposite extreme of a magnetic field frozen into the crust, and even so we will see that the observed evolution is very limited.At the crust-core interface, we also assume that the radial component of the fluid velocity  vanishes.Thus, from equation ( 27), we obtain On the other hand, in the strong-coupling regime, the radial velocity component can be determined solely from the continuity equations ( 25)-( 26), leading to the expression where we introduced a length scale with the baryon density given by   =   +   .As a consequence, at the crust-core interface, the chemical equilibrium condition is satisfied, i. e.
Further boundary conditions are related to the axial symmetry constraint.In order for the velocity field and the magnetic field to be continuous at the axis, one must impose   =   = 0 and   =   = 0 at  = 0 and  = .The latter conditions imply / ( = 0, ) = ( = 0, ) = 0. Hence,  must be a constant along the axis, whose value we set to zero, so that the magnetic flux enclosed by any circle around the axis at fixed  and  can be obtained as 2(, ).

Short-term evolution through artificial friction
In Castillo et al. (2020) and in the present work, in order to focus on the long-term evolution, the inertial terms in the equations of motion for neutrons and charged particles are neglected, and the short-term dynamics associated to the excitation and relaxation of sound, gravity, and Alfvén waves is mimicked by introducing an artificial friction force,   = −     , acting on the neutron fluid.
In Castillo et al. (2020), which considered the weak-coupling regime, the aim was to start with an arbitrary initial magnetic field configuration, and then, by properly choosing , allow the bulk motions controlled by   to reduce the net force imbalance on a fluid element (in parenthesis on the right-hand side of equation 11) much more quickly than ambipolar diffusion (with relative velocity   ) reduced the force imbalance on the charged-particle component (in parenthesis on the right-hand side of equation 12).This implied that  had to fulfill  ≪     , while at the same time being large enough for a feasible numerical simulation.Thus, Castillo et al. (2020) presented a detailed hierarchy of time scales that mimic the relaxation due to sound, gravity, and Alvén waves propagating through the NS core with artificial friction.The first two represent the typical time scales in which the density perturbations grow and the poloidal hydromagnetic quasi-equilibrium is reached.
In the new scheme presented here, however, the time derivative terms in the continuity equations are neglected, which is valid for times longer than the gravity or buoyancy time scale (   , in the notation of Castillo et al. 2020).This approach speeds up the simulation code and allows to consider realistic magnetic field strengths since the simulations now start after this buoyancy-like time scale, when the (non-fictitious) poloidal forces are already close to a hydromagnetic equilibrium (  Pol  +   +   ≈ 0).Nevertheless, from this new starting point, the toroidal Lorentz force,  Tor  , is still finite and balanced by fictitious friction.It must decay to zero by toroidal fluid motions, with velocity  Tor ∼  2 /(4   ℓ  ) (see equation 27), generated by the magnetic field itself, which occurs on a longer, Alfvén-like time scale where ℓ  is the typical length-scale of the magnetic field.

Long-term evolution through Urca reactions
As we have previously discussed, the short-term dynamics, mimicked by the fictitious friction, evolves the NS core towards an almost complete balance of forces, leaving it out of chemical equilibrium.The excess of chemical energy can only be dissipated by Urca reactions, which in turn locally change the degeneracy pressure forces, so that the magnetic field must rearrange itself, moving particles in such a way that a new hydrostatic equilibrium is reached.As we shall explain in detail in this section, this process changes the typical time scale at which the core reaches chemical equilibrium, and the time scale for the decay of Δ can be much longer when considering the effect of the magnetic field dynamics than in the absence of a magnetic field (Reisenegger 2009).We emphasize this because various authors (e. g., Thompson & Duncan 1996;Lander et al. 2021) have assumed that, since the initial NS temperature is so high, the stellar matter can be regarded to be in chemical equilibrium and thus to be described as a barotropic fluid.As we shall show in § 3.3, this assumption imposes strong, but unrealistic, constraints on the initial magnetic field configuration.
By studying the time derivative of Δ, the typical time scale of the long-term magneto-chemical evolution can be derived.Using the continuity equations ( 25) and (26) (but now retaining the timederivatives of the number densities on their left-hand sides), as well as equation ( 5), one obtains This expression shows that, if the magnetic field suddenly disappeared and there were no fluid motions, Δ would decrease expo-  nentially, reaching chemical equilibrium on a time scale where we used   ≳   , |  |, and evaluated   and  at the center ( = 0) to give the numerical estimate.However, when the magnetic and velocity fields are present, there are fluid motions that try to keep the force balance and largely compensate the effect of the Urca reactions.The time scale of this coupled magnetic and chemical (magneto-chemical) evolution follows from the continuity equations, which yield a characteristic fluid velocity that can be estimated from equation (32) as  ∼ |Δ|ℓ  /  .On the other hand, considering zero artificial friction in equation ( 11), the degeneracy pressure forces plus gravity balance the poloidal Lorentz force, (as we will see in § 3.2, this assumption is not strictly correct, but here we only aim to provide a rough estimate of   ), so that   |Δ| ∼  2 /4, the characteristic time scale of this process is given by yr (see also Thompson & Duncan 1996;Reisenegger et al. 2005;Reisenegger 2009, but note that the first of these references incorrectly assigns this time scale to the low-temperature, "weak-coupling" regime).
We note that the ratio between the latter two time scales is roughly where we used  2    ∼   /  and again evaluated at the center ( = 0) to give the numerical estimate.It can be seen that, even for magnetar-strength magnetic fields, the time scale for the combined magneto-chemical evolution is expected to be much longer than the naive estimate   .
Here, we have consistently distinguished between the length scales ℓ  and ℓ  in equation ( 38) and (39) because, for an arbitrary initial magnetic field configuration, ℓ  is typically larger than ℓ  .This stems from the fact that the radial profiles of the densities,   () and   () are not identical, but quite similar, smoothly decreasing over a lengthscale ℓ  ∼    (see equation 33), while ℓ  is generally only a fraction of    .
In order to identify the transition temperature to the weak-coupling regime, we set   , evaluated at  = 0, equal to the ambipolar diffusion time, as given by Castillo et al. (2020), and also evaluated at where the correction factor 0.1 was included in consistency with the findings of Ofengeim & Gusakov (2018) and Castillo et al. (2020).Estimating ℓ  ≈    /4 = 2.8 km and ℓ  ≈ 4ℓ  ≈    , this yields the transition temperature Lastly, as in the weak-coupling regime, the value of  is not arbitrary.It must be set so as to separate the short (unphysical) and long (physical) time scales,    ≪   , which requires  2. The lines represent the poloidal magnetic field, labeled by the magnitude of   , and in colors the poloidal current function   .These variables are plotted in code units, see Table .1. Table 2. Parameters for the initial magnetic field configurations used in our simulations.The normalization constants  0 and  0 are fixed by the condition ⟨ Pol ⟩ = 1, and ⟨ Tor ⟩ = 1, respectively, where ⟨.⟩ denotes an rms average over the stellar core.Here, we choose the poloidal potentials as , where   ℓ (cos  ) is the associated Legendre polynomial of order ℓ with azimuthal index ; and the coefficients,  1 and  2 , are determined so that the azimuthal current vanishes at the stellar surface,   ( = ,  ) = 0, and  matches the core magnetic field with the external multipolar expansion corresponding to a current-free magnetic configuration (see § 2.3) (Armaza et al. 2015).The constants √ 0.6 and √ 0.4, multiply the potentials in order to set 60% and 40% of initial poloidal and toroidal magnetic energy, respectively.These models are shown in Fig. 2. Now, with the aid of equation ( 32), the radial component of the fictitious force can be written as for which equation ( 42) also implies Thus, the condition in equation ( 42) not only guarantees a clear separation of the short and long time scales, but also consistently implies that, at all times, the poloidal component of the fictitious force is much smaller than the non-fictitious poloidal forces, which therefore are close to balancing each other.

EVOLUTION AT CONSTANT TEMPERATURE
The dimensionless code units are summarized in Table 1.These units have been selected to ensure that the system of equations ( 23)-( 28) becomes independent of the magnetic field strength and temperature.
As the time scales    ,   , and   depend on the number densities and other radial functions, hereafter, we use as references for them the expressions and In dimensionless code units, they read as where, as before, we estimated ℓ  /   ≈ 1/4, and ℓ  ≈ 4ℓ  ≈    .We remark that in the dimensionless units used to express  (see Table .1), it must satisfy the condition  ≪ 4 0 / 0 ≈ 0.4 in order to maintain consistency with equation ( 42).We also note that the time   , which in code units becomes proportional to  2  and very small for realistic values of this quantity, does not enter the evolution of the magnetic field at constant temperature, but will become relevant when analyzing the evolution with variable temperature in § 4.3.
As an initial approach to the problem, we used the scheme explained above to evolve the set of equations ( 23)-( 28) at constant temperature in order to study how the NS core reaches chemical equilibrium on the expected time scale.Also, here and hereafter, as initial condition for all the simulations we set the magnetic energy in the NS core to be 60% poloidal and 40% toroidal.

Dependence on artificial friction
In order to separate the time scales in equations ( 47) and (49), we must set  ≪ 1 in dimensionless code units, which is equivalent to equation ( 42) in physical units.In order to make this condition more precise, we note that the long-term evolution must remain independent of the value of .We tested three different values, namely  = 10 −2 , 10 −3 and 10 −4 , which yield the time-scale ratios   /   = 37, 370, and 3700, respectively.The results are summarized in Fig. 3.
In panels (a), (c), and (e), we see that the curves for  = 10 −3 and  = 10 −4 are almost identical from the beginning of the simulation, as expected for small enough , since this causes the non-fictitious poloidal forces to be near hydromagnetic quasi-equilibrium from the very beginning, as discussed in § 2.4.2.More specifically, the good agreement between these curves can be interpreted in the sense that, for small enough , the poloidal force balance (the poloidal component of equation 27) determines the chemical imbalances, but not the poloidal velocity field  Pol , which is determined from the chemical imbalances by the continuity equations ( 25) and ( 26).This good agreement is not seen for  = 10 −2 , suggesting that this value is not small enough for the fictitious force to be ignored in the force balance.
For the toroidal variables, the results are very different.Panel (b) shows that, during the early stages ( ≲    ), the toroidal magnetic energy decreases by the same amount for all values of , but on the different time scales    associated to each value of .On the other hand, the toroidal velocity initially scales ∝  −1 , and hence presents a more prominent evolution for the smaller values of .This happens because, for an arbitrary magnetic field configuration, the toroidal Lorentz force is only balanced by the corresponding component of the fictitious friction force.Thus, the magnetic field has to rearrange itself by moving the fluid in the toroidal direction with velocity  Tor =  Tor  /(   ), until the toroidal hydromagnetic quasi-equilibrium is established on a typical time scale ∼    .
In panels (c), (d), and (e), we see that the velocity components and chemical imbalances strongly decrease towards late times, signaling an approach to hydromagnetic and chemical equilibrium.At this point, these variables become dominated by numerical noise (to a large extent due to the near-cancellation of forces in equation 27), causing discrepancies between the curves for different values of .
Nevertheless, given the otherwise good agreement of the results with  = 10 −3 and 10 −4 , we will hereafter use the value  = 10 −3 , which causes an acceptable separation of the time scales and is numerically cheaper than  = 10 −4 .

Magneto-chemical evolution
In this section, we present and discuss the results of our simulations for the "magneto-chemical" evolution at constant temperature, again using the results of Model 2 to describe some generic aspects of our simulations.The results for the other models (see Table 2 and Fig. 2) are analyzed later in § 3.3.In any case, we have checked that the results for the magneto-chemical evolution are qualitatively the same for all the other models.
The magneto-chemical evolution of Model 2 is presented in Figs. 4  and 5, which show the spatial dependence of different variables at selected times, and in Fig. 6, which shows the time evolution of integrated or rms-averaged quantities.The initial state is clearly not a hydromagnetic quasi-equilibrium, as seen by the fact that the toroidal magnetic field is finite on open poloidal field lines, causing a substantial toroidal Lorentz force, which is balanced only by the corresponding component of the fictitious friction.At the same time, the chemical potential perturbations of the neutrons and charged particles are very different from each other (both in magnitude and in their spatial dependence), implying a substantial departure from chemical equilibrium.
Nevertheless, the initial poloidal component of the artificial friction force,   Pol , is already much weaker than the other poloidal forces, i. e., the simulations start with a very well-established poloidal quasi-equilibrium.We also note the solenoidal appearance of the poloidal velocity field, a consequence of the condition obtained from adding equations ( 25) and ( 26  .However, it can be understood as follows.Since   Pol is from the beginning smaller than the other poloidal forces, the poloidal force balance equation can be approximated as  the chemical potential perturbation of the charged particles; and the chemical potential perturbation of the neutrons.All quantities are given in code units, as defined in Table 1.
.17 0 5.17  .All quantities are given in code units, as defined in Table 1.
On the longer time scale   , Urca reactions become effective, reducing the magnitude of the chemical imbalance Δ (see Fig. 6[a]) by changing the chemical potential perturbations   in order to equate them (Fig. 4[d]), in that way modifying the local poloidal force balance.Therefore, the magnetic field slowly rearranges itself, so that the evolution proceeds through a continuum of consecutive hydromagnetic quasi-equilibria.
As the NS core approaches the final chemical equilibrium state, the chemical potential perturbations become equal, thus   =   and the poloidal artificial friction force  Pol  (orange).The vertical lines correspond, from left to right, to the time scale    (dashed-dotted), and   (dashed), respectively.Here, ⟨ ⟩ denotes the root mean square (rms) average over the volume of the core.All quantities are given in code units, as defined in Table 1.
When the NS core approaches the final magneto-chemical equilibrium state, the velocity decreases substantially in comparison to its initial values, while the chemical potential perturbations become roughly equal.This can be seen in Fig. 4 where, at the final snapshot ( =   ), the velocity field has decreased by two orders of magnitude with respect to the initial state and   ≈   .The strange structures seen in Fig. 4(d), second column, are caused merely by numerical noise due to truncation error in the quasi-cancellation of the forces, suggesting that the final magneto-chemical equilibrium has been reached (see also Fig. 6 [a] and [b]).

Grad-Shafranov equilibria
Figures 5(b,c) and 6(b) show that, within a few times    , a hydromagnetic quasi-equilibrium is reached, in which the fictitious force becomes negligible and all the real forces are close to balancing each other.This implies that the toroidal magnetic force vanishes,

𝒇 Tor
≈ 0, because there is no toroidal pressure gradient or gravitational force available to balance it.Thus, from equation ( 15), ∇ ∥ ∇, so one potential is (at least locally) a function of the other,  = ().One can see this in the first panel from left to right of Fig. 4, where the color scale, indicating the values of , come progressively closer to matching the field lines, which correspond to surfaces of constant .Perhaps even more clearly, Fig. 7 shows that at early stages there is no clear relation between  and , but at  ≳    there is an evident dependence.Once this happens, equation ( 14) can be rewritten as (Here and below, primes denote derivatives with respect to .)At later times ( ≳ 20    ),   also becomes roughly a function of , because Urca reactions have reduced the chemical potential perturbations of the charged particles,   , so the poloidal Lorentz force is mostly balanced by the neutrons,  Pol  ≈ −   =   ∇  , as seen in the last two rows of Fig. 5. Fig. 7 also shows that, at  ∼   , there is a clear relation   ≈   (), and the NS core is very close to chemical equilibrium,   ≈   , so the matter in the core behaves as a barotropic fluid.Therefore, taking into account that, first,  becomes a function of , and later,   =   ≡ (), the poloidal force balance equation, equation ( 50), reduces to a "Grad-Shafranov (GS) equation" (Grad & Rubin 1958;Shafranov 1966): We emphasize that this equation is not satisfied in the previous stages of the evolution, which are also hydromagnetic quasi-equilibrium configurations, but out of chemical equilibrium.In fact, as we shall see in § 4.4, this state is never reached in the actual evolution of a NS core in the strong-coupling regime.

Magnetic energy dissipation
In this section we discuss the magneto-chemical evolution in terms of the magnetic energy dissipated in the core.We start by deriving the expression for the time-derivative of the magnetic energy contained in the NS core for the most general case, i. e., considering both the effects of ambipolar diffusion and non-equilibrium Urca reactions, and also allowing a non-linear dependence of ΔΓ on Δ/(  ) (see e. g.Reisenegger 1995;Fernández & Reisenegger 2005;Gusakov et al. 2005).
Taking the time derivative of the total magnetic energy contained in the core,   = ∫ V  2  3 /(8), where V is the core volume; and using equation ( 8), one obtains Integrating by parts, where S is a surface element outward normal to the surface  defined by the crust-core interface.We identify the first term on the right-hand side of this expression as the mechanical work done by the fluid against the Lorentz force, and the second term (without the minus sign) as the outward Poynting flux integrated over the surface, which we denote hereafter as   .From our equations in § 2.1, one obtains where are the power dissipated by non-equilibrium Urca reactions, ambipolar diffusion, and artificial friction, respectively.When deriving equation ( 57), the terms proportional to / 2 were disregarded, in constancy with Newtonian gravity.Nonetheless, we computed them numerically and checked their insignificance when compared to the other terms.Therefore, the rate of change of the magnetic energy stored in the core reads as Figure 8 shows all the terms that dissipate energy in our reference simulation for the strong-coupling regime, in which   = 0 and ΔΓ = Δ.As might have been expected, initially the fictitious friction dissipation   dominates, since the dynamics is given by quick toroidal fluid motions with velocity  Tor =   Tor /(   ), which unwind the magnetic field around  ∼    .This occurs at the position of the largest peak of (  /) , corresponding to where most of the magnetic energy dissipates.At later times (≫    ), the real physical evolution begins, when non-equilibrium Urca reactions become effective and the chemical power released by them,    , is the main source of magnetic energy dissipation in the second, lower peak.A curious fact of our simulations is the opposite sign of the Poynting flux  P with respect to the other dissipative terms, also seen in Fig. 8, which indicates that external magnetic energy (although a small amount) is being injected into the NS core.One might have anticipated such behavior from the fact that outside the core we are imposing a current-free magnetic field configuration, so the magnetic dissipation occurs only in the core.This is analogous to a wire carrying an electric current: Since most of the magnetic energy is located outside, but the dissipation takes place inside, there is a flux of magnetic energy from outside into the wire.
Our simulations also allow us to determine how much magnetic energy is dissipated in order to establish the final GS equilibrium for different initial conditions.Fig. 9(a,b) shows the evolution of the internal poloidal and toroidal magnetic energies for the different models we simulated.These models are specified in Table 2 and Fig. 2).
The results show that during the early evolution,  ≲    , most of the toroidal energy is dissipated by almost the same amount for all the magnetic field models.This strong dissipation is due to toroidal fluid motions generated by the magnetic field in order to reach the hydromagnetic quasi-equilibrium and thus satisfy the condition  = ().As these toroidal fluid motions proceed, they are opposed by the artificial friction force, so that most of the dissipation is caused initially by the artificial friction,  Tor  ≫  H (Fig. 9[c]), where  Tor  is computed considering only the toroidal component of the velocity (see equation 60).We remark that this pronounced initial evolution of the toroidal energy is a consequence of starting with an initially uncompensated toroidal magnetic field.If the simulations were to begin with a quasi-equilibrium field configuration ( = ()), the toroidal field would not undergo significant changes.This can be corroborated in Fig. 9(d), where we computed the toroidal magnetic flux considering two choices for the surface (S) of integration: the entire meridional cross-section of the NS core, and the piece corresponding to the closed poloidal field lines.
The results show that, as expected, the evolution of the total toroidal flux is particularly significant around  ∼    as the core approaches hydromagnetic quasi-equilibrium.Conversely, the toroidal flux associated with the regions of closed poloidal field lines appears to remain nearly conserved throughout the entire evolution.Therefore, if the initial magnetic field were to be close to a hydromagnetic M agnetic energy dissipation   ) (purple), which is used to check energy conservation, as it should be equal to   /.We note that all these variables are multiplied by time, so they represent energy released per unit ln .The vertical lines correspond, from left to right, to the time scales    (dashed-dotted) and   (dashed), respectively, while the horizontal solid line corresponds to zero.
quasi-equilibrium state, the toroidal field would indeed undergo minimal evolution in order to reach the final GS equilibrium state.All in all, by far the main effect of the artificial friction is to eliminate the toroidal magnetic field in the regions of open poloidal field lines.
On the other hand, during this early evolution, the poloidal energy evolves as it is more sensitive to Urca reactions, which operate on the time scale ∼   .
During the late evolution,    ≲  ≲   , the toroidal magnetic energy evolves very little and settles to a common asymptotic state for all the magnetic field models.For the poloidal energy, however, most of the dissipation occurs in this late evolution, when Urca reactions operate, as can be seen in Fig. 9(a) and, more clearly, for the simulation with  = 10 −4 in Fig. 3(a,b).One can see also that, for all initial configurations, only a small fraction of the poloidal energy is dissipated.This is particularly evident in Model 1, which corresponds to the simplest dipole model.On the opposite extreme, Model 3 exhibits the most pronounced decrease in poloidal energy, likely due to its less symmetric nature.
In terms of the dissipation, at this late stage of the evolution, the artificial friction should be irrelevant, i. e.,  H ≫  Pol  and  H ≫  Tor  , as the system has already reached the hydrostatic equilibrium, and since  Tor  ∝ (∇ ×∇) 2 .Unfortunately, numerical noise affects the last moments of our simulations ≳   , and  Tor  starts to be dominant again (see Fig. 9 [c]), however at this point the core is very close to chemical equilibrium and all dissipation effects become very small.Therefore, we conclude this section by noticing that our results suggest that, for relatively smooth initial magnetic field configura-tions in axial symmetry, reaching the GS equilibrium requires dissipating only a few percent of the poloidal and and most of the toroidal magnetic energies.

MAGNETO-THERMAL EVOLUTION
Up to this point, we have solved equations ( 23)-( 28) in the unphysical situation where the core temperature remains constant, i. e., when  does not change in time.However, in reality, the core temperature evolves, and since  depends strongly on the temperature, it will also evolve.Below, we describe an efficient strategy to include the thermal evolution, using the output of the simulations at constant temperature and incorporating the thermal evolution by reparameterizing the time variable.

Temperature evolution
Before describing the procedure, we introduce the equation needed to evolve the core temperature, which is given by a thermal balance equation in the isothermal approximation (Thorne 1977;Yakovlev & Pethick 2004) where   is the total power released by heating mechanisms,   is the total neutrino luminosity,   is the luminosity associated to the photons emitted from the stellar surface, and  is the total heat capacity of the NS core, which for our stellar model is given by where  9 = /(10 9 K).
In the present case, equation ( 63) takes the form where the only heating mechanisms are the chemical power released by non-equilibrium Urca processes,    ; and the artificial friction dissipation,   (see equation 61), which was included for consistency.We neglected the photon luminosity because this cooling mechanism is only relevant for later stages in the NS life (Yakovlev & Pethick 2004).
The expressions for the neutrino luminosity and the chemical power released by non-equilibrium Urca reactions are (Haensel 1992) where  *  is the equilibrium neutrino emissivity, which is a function of the stellar model and temperature, and  () and  () are dimensionless control functions of  = Δ/(  ).The exact expressions for these functions were found by Reisenegger (1995) for both direct and modified Urca processes.Consistently with our stellar model, we assume only modified Urca process, for which these functions for arbitrary values of  are given by Since in the present problem we will mostly have | | ≪ 1 (the "subthermal" regime) as already assumed in equation ( 21), the net heating rate becomes where  *  is the equilibrium neutrino luminosity, which for our stellar model reads Equation ( 70) shows that non-equilibrium Urca reactions cause enhanced cooling as a consequence of the magnetic feedback (see also Fernández & Reisenegger 2005).However, we will show that no significant magnetic feedback takes place even for very strong magnetic fields.
Incorporating equation ( 65) directly into the existing system of equations ( 23)-( 28) for our simulations is the most natural approach in principle.However, for enhanced efficiency, we employ an alternative strategy outlined in § 4.2.

Time reparametrization
Our previous simulation results, calculated at constant temperature, can be translated into a more realistic evolution with variable temperature, and thus time-dependent (, ), by noting that the system of equations ( 23)-( 28) remains invariant under the following changes of variables where primed variables ( ′ ,  ′ ,  ′ ,  ′ ,  ′ ) correspond to the previous simulations at an arbitrary constant temperature  ′ , while unprimed variables (, , , , ) take into account the thermal evolution and its effect on the parameters.(As seen in equation 22, the functions  and  ′ have the same dependence on density, and thus on , which thus cancels out in their ratio.Therefore, we wrote the latter simply as  ′ /(), which depends only on temperature, and hence on time.)Thus, it is possible to solve the system of equations ( 23)-( 28) for time-independent  ′ () and  ′ , then solve for the evolution of the temperature in order to obtain the time dependence of (, ), and finally reparametrize the time variable in order to obtain the evolution of all variables for the realistic case of variable temperature.Our selfconsistent approach can be understood as follows: (i) We run the simulations that evolve the set of equations ( 23)-( 28) at constant temperature (thus time-independent  ′ and  ′ ) as presented in § 2, calling the time variable  ′ .
(ii) Then, using these results, we solve the equation for the temperature as a function of  ′ , which with the aid of equation ( 70) reads as where  * ′  is the equilibrium luminosity evaluated at the constant reference temperature  ′ .Here, we also used / ′ = (/ ′ ) 6 , and the scaling in equations ( 72)-( 74), so that    =  ′   (/ ′ ) 6 ,   =  ′  (/ ′ ) 6 , and  *  =  * ′  (/ ′ ) 8 where the primed variables are evaluated at the fixed temperature  ′ , not at the true, time-dependent temperature .
(iii) Finally, we obtain the physical time variable that includes the effects of the temperature evolution by integrating equation ( 72).Thus, we can plot any variable of interest as a function of , the real physical time.

Analytic solution for passive cooling
If there is no substantial feedback (heating or enhanced cooling) due to the magnetic field evolution, the thermal evolution equation ( 65) reduces to / = − *  ()/ (), or equivalently  9 / = − 7 9 /(1.32yr),so the temperature evolves as where   is the initial temperature, and is a characteristic cooling time at  =   .
For the initial temperature used in our simulations,   = 2 × 10 10 K,   = 3.4 × 10 −9 yr, and, for any realistically large initial temperature,   ≫    ≈ 5 × 10 8 K, the transition from the strong-coupling to the weak-coupling regime occurs at    ∼ 30 yr.
Applying the time rescaling of equation ( 72) with the temperature given by equation ( 76), the physical time variable  for the evolution with variable temperature becomes an exponential function of the auxiliary time  ′ used for the constant-temperature evolution, (Here and below, all time scales with primes refer to times in our constant-temperature simulations of § 3.) Writing  ′ in terms of the code units defined in Table 1, this reduces to which is plotted in Fig. 10 for two large values of   .Clearly, given by equation ( 78), for two different initial rms magnetic field strengths:   = 2 × 10 16 G (blue) and   = 5 × 10 15 G (red).The dashed vertical line marks the time scale  ′   , which is independent of   .The horizontal lines show, from bottom to top, the time scales    =  ( ′   ) (dashed) mapped from equation ( 79) for the different initial rms magnetic field strengths with their respective colors, and the transition time    from the strongcoupling to the weak-coupling regime (black-solid).The purple-hatched area displays the time in which GS-equilibrium is reached in our simulations at constant temperature, and the green-hatched area shows the time range corresponding to the weak-coupling regime.moderate values of  ′ (≲ 1) and   (≲ 10 15 G) yield huge values of , many orders of magnitude larger than    (or even than the age of the Universe), indicating that very little evolution will occur during the strong-coupling regime,  ≲    .Only for much larger field strengths, a substantial evolution might occur.
As discussed in § § 2.4 and 3.1, the initial stages  ′ ≲  ′   of our constant-temperature simulations are dominated by the fictitious friction force, which allows to set up a hydromagnetic equilibrium that in a real NS core would be reached in a few Alfvén times.Therefore, only for later times,  ′ ≫  ′   , the simulations can be expected to represent the true astrophysical evolution.For  ′   to correspond to a physical time earlier than the transition to the weakcoupling regime, we must have where  ′ = 10 −3 was the dimensionless friction coefficient used in most of our simulations.Thus, our simulations represent the astrophysically relevant regime only for very strong magnetic fields satisfying these conditions.Figure 10 illustrates the correspondence between the time scale  ′   and real time variables    , highlighting the increasing magnitude of the latter values as we diminish the initial rms magnetic field strength.
On the other hand, we have seen in our constant-temperature simulations (Figs. 3,6,7,8,9) that the time needed to reach the final chemical and GS equilibrium typically falls in the range  ′  ≈ 0.1 − 1 in code units, depending on the magnetic field model.Thus, the strength of the initial magnetic field required to reach this equilibrium before however, that these differences are only a factor ∼ 2−3 in the constanttemperature time variable  ′ , whereas the exponential dependence of  on  ′ in equation ( 78) amplifies these time differences in the time variable , making the magneto-thermal evolution for very strong   quite sensitive to the magnetic field geometry.

CONCLUSIONS
In this paper, we have reported the first simulations of the magnetothermal evolution of a NS core in the strong-coupling regime applicable at high temperatures,  >    ≈ 5 × 10 8 K, covering two relevant aspects not considered before, namely the interplay between non-equilibrium Urca reactions and magnetic field dynamics as well as the effect of the thermal evolution on the evolution of the magnetic field.The latter was included by starting from the results of the simulations at constant temperature and implementing an efficient numerical strategy to translate these results to the realistic case of evolving temperature.Like previous work (Castillo et al. 2017(Castillo et al. , 2020)), these simulations are done in axial symmetry, for a two-fluid (neutrons and charged particles) NS core surrounded by a non-conducting medium, and include a fictitious friction force that allows the system to quickly reach a hydromagnetic equilibrium state while ignoring the inertial terms and, in this case, the time derivatives in the particle conservation laws, that would substantially slow down the simulations.The main conclusions can be summarized as follows: (i) In the initial stages, in which the fictitious friction is important, the NS core evolves mostly through toroidal motions that unwind the magnetic field lines, leading to a hydromagnetic equilibrium state in which the toroidal component of the Lorentz force vanishes, confining the toroidal magnetic field component to the region of closed poloidal field lines.In a real NS core, this process is expected to happen very quickly, within a few Alfvén times.
(ii) In the constant-temperature simulations, the magnetochemical evolution generated by the interplay between Urca reactions and the magnetic field dynamics leads the NS core to evolve towards a chemical equilibrium state in which the fluid becomes barotropic and the magnetic field satisfies the non-linear Grad-Shafranov equation.This state is reached in a time which for realistic magnetic fields is much larger than the previous, traditionally underestimated value   provided by equation (37) (Thompson & Duncan 1996;Lander et al. 2021).Moreover, our results indicate that, once the fluid has settled into a hydromagnetic quasi-equilibrium with vanishing toroidal Lorentz force, the magnetic field needs to evolve very little in order to reach the final GS equilibrium state.This may appear surprising at first glance, as one would intuitively expect a significant magnetic evolution accompanied by a strong dissipation during the transition from a wide range of possible equilibria, where two fluid forces proportional to the gradients of two independent scalar functions,   = −  ∇  and   = −  ∇  , balance the poloidal Lorentz force, to the GS equilibrium state, where the latter force is balanced by a force proportional to the gradient of a single scalar function,   +   = −  ∇.Possibly, the results might be different if we considered more complex magnetic field structures or relaxed the axial-symmetry constraint.
(iii) Our results also revealed an intriguing behavior related to the fluid forces: for an arbitrary initial magnetic field, the fluid forces (pressure plus gravity) exerted by the charged particles and the neutrons were found to approximately balance each other,   ≈ −   , and were considerably stronger than the poloidal Lorentz force  Pol  .This observation may seem unexpected since the chemical potential perturbations responsible for the fluid forces originate from  Pol  .However, this phenomenon arises as an indirect consequence of stable stratification.Specifically, it stems from the fact that the radial density profiles,   () and   (), although not identical, exhibit similar smooth decreases over a characteristic length scale ℓ  ∼    , substantially larger than the magnetic length scale ℓ  .
(iv) When the temperature evolution was included, our results showed that even the strongest magnetic fields we considered ( ≳ 10 16 G) cannot generate an appreciable magnetic feedback on the thermal evolution.Consequently, the NS core cools passively, so the effectiveness of Urca reactions gradually declines, the overall magneto-chemical evolution decelerates, and there is no significant magnetic evolution before ambipolar diffusion becomes more important than Urca reactions, thus ending the strong-coupling regime.This was quantified by showing that the physical time  (corresponding to evolution with realistically decreasing temperature) depends exponentially on the time variable  ′ used in the simulations at constant temperature (equation 78).Therefore, for realistic magnetic fields, the time needed to reach the GS equilibrium state would be many orders of magnitude larger than the time to transit to the weak-coupling regime.
Finally, we note some issues that were not included in the present work and might lead to future extensions: (i) Consistently with our stellar model, we neglected the direct Urca process, which might happen in the inner core of massive NSs, speeding up the thermal evolution and thus strongly reducing the time at which the strong-coupling regime ends.The magnetochemical evolution would also be accelerated proportionally in the inner core, but not in the outer layers, where only modified Urca reactions operate.As in the case considered here, the thermal evolution for realistic magnetic field strengths will still be much faster (by a factor ≳ 8/ 2 ) than the magneto-chemical evolution, therefore the chemical and GS equilibrium should not be reached.
(ii) As discussed in the introduction, neutrons and protons are expected to become superfluid and superconducting, respectively, relatively early in the NS life, when the core is at temperatures  ∼ 10 8 − 10 10 K (Migdal 1959;Page et al. 2013).The inclusion of these effects to study the long-term magneto-thermal evolution can be highly complex as the macroscopic manifestation of the interactions between quantized neutron vortices and magnetic flux tubes needs to be taken into account.Qualitatively, one can anticipate that the approach to chemical equilibrium would be slowed down as Cooper pairing drastically reduces the reaction rates (Yakovlev et al. 2001;Yakovlev & Pethick 2004).Additionally, it impacts the cooling rate, initially accelerating it through Cooper pair breaking and formation processes, and eventually delaying it through the suppression of the Urca reactions.Furthermore, disregarding the interactions between vortices and flux tubes, superfluidity and superconductivity suppress inter-particle collisions, potentially resulting in faster ambipolar diffusion.This, in turn, may cause the transition from the strong-coupling regime to the weak-coupling regime to occur earlier.Therefore, including these effects might reinforce the conclusion that no substantial magnetic field evolution will occur during the strongcoupling regime.
(iii) Throughout this work, for simplicity we have used Newtonian gravity.However, the inclusion of general relativity can have significant implications for the overall magneto-thermal evolution.Specifically, when accounting for it in the thermal evolution, it becomes necessary to distinguish between local values of the temperature and the luminosities and their redshifted values as detected by a distant observer.In particular, the two temperatures differ by a factor  −  ∼ 1.3, so the equilibrium neutrino luminosity due to Urca reactions, which is proportional to  8 could change by a non-negligible factor  −8 ≈ 8.2 depending on which temperature is used, with an important quantitative effect on the thermal history (see e. g., Fig. 7 and Fig. 8 in Ofengeim et al. 2017).

Figure 1 .
Figure 1.Schematic representation of the self-consistent method implemented to evolve the magnetic field in the strong-coupling regime.

Figure 2 .
Figure 2. Initial magnetic field configurations used in our simulations given by the potentials listed in Table 2.The lines represent the poloidal magnetic field, labeled by the magnitude of   , and in colors the poloidal current function   .These variables are plotted in code units, see Table.1.

Figure 3 .
Figure 3.Comparison of the evolution of different simulations of the evolution at constant temperature using the same initial condition, Model 2 in Table2, for three different values of  = 10 −2 (red), 10 −3 (blue), and 10 −4 (black), corresponding to the ratios   /   = 37, 370, and 3700, respectively.All quantities are given in code units, as defined in Table1.The vertical lines, from left to right, show the values of the time scales    , with their respective colors in solid for the different  of each simulation, and   (dashed).The panels show the time evolution of the poloidal fraction of the NS core magnetic energy in (a); the toroidal fraction of magnetic energy in (b); the root mean square (rms) averages over the volume of the core, ⟨ ⟩, of the poloidal and toroidal components of the fluid velocity in (c) and (d), respectively; and of the chemical imbalance, Δ in (e).
), ∇ • (  ) = 0, where   ≡   +   is the baryon number density.Figs.5(a) and 6(b) show that initially   ≈ −   , i. e., the fluid forces (pressure plus gravity) on the charged particles and the neutrons roughly balance each other and are substantially stronger (by a factor ∼ 10) than the poloidal Lorentz force  Pol  .This behavior is at first sight quite surprising, because the chemical potential perturbations causing   and   are due to  Pol

Figure 4 .
Figure 4. Magneto-chemical evolution at constant temperature for Model 2 with  = 10 −3 , corresponding to   /   = 370.All panels are meridional cross-sections of the star.Rows (a), (b), (c), (d) correspond to different times,  = 0,    , 10    , and   , respectively.In each row, the four panels display, from left to right: the magnetic field configuration, where lines represent the poloidal magnetic field, labeled by the magnitude of , and colors representd the poloidal current function ; the poloidal component of the velocity field, ;the chemical potential perturbation of the charged particles; and the chemical potential perturbation of the neutrons.All quantities are given in code units, as defined in Table1.

Figure 5 .
Figure 5. Evolution of the force densities the same simulation presented in Figs. 4 and 6.Rows (a), (b), (c), and (d) show snapshots at different times,  = 0,    , 10   , and   , for the different force densities: From left to right, the toroidal Lorentz force  Tor , the poloidal Lorentz force  Pol  , the charged particle force   , the neutron force   , and the poloidal artificial friction force  Pol  .All quantities are given in code units, as defined in Table1.

Figure 6 .
Figure 6.For the same simulation as in Figs. 4 and 5, panel (a) shows the rms values of the chemical potential perturbations    (blue),    (red), and Δ (green); and panel (b) shows the rms values of the different forces: the poloidal Lorentz force  Pol  (green), the charged-particle force   = −  ∇   (red), the neutron force   = −  ∇  (blue), the toroidal Lorentz force  Tor  (purple), which is equal to minus the toroidal artificial friction force −   Tor ,

Figure 7 .
Figure 7. Scatter plot for the same simulation shown in Figs. 4 to 6, with  = 10 −3 and initial conditions given by Model 2, with the variable on the vertical axis written on top and the variable on the horizontal axis written on the bottom of each column: (a)  as function of , (b)   as function of , (c)   as function of , (d)   as function of   .We show the relations at all the grid points at  = 0, 0.1   ,    , 10   , 20   , 40   and   .

Figure 8 .
Figure8.Magnetic energy dissipation for the same simulation as illustrated in the previous figures, whose initial condition is Model 2 and with  = 10 −3 , corresponding to   /   = 370.The plot shows all the dissipation terms (equation 61): the time derivative of the magnetic energy inside the core, (  / ) (blue); the chemical power released per unit time due to nonequilibrium Urca reactions,  H  (red); the artificial friction dissipation,    (green); the Poynting flux,    (orange); and the combination − (  +   +   ) (purple), which is used to check energy conservation, as it should be equal to   /.We note that all these variables are multiplied by time, so they represent energy released per unit ln .The vertical lines correspond, from left to right, to the time scales    (dashed-dotted) and   (dashed), respectively, while the horizontal solid line corresponds to zero.

Figure 9 .
Figure 9. Evolution at constant temperature for the different initial magnetic field configurations listed in Table2(with color coding as shown in panel [a]), always with  = 10 −3 , showing: (a) poloidal and (b) toroidal magnetic energy, both in units of the total initial magnetic energy; (c) quotient between artificial friction dissipation, split into a poloidal and a toroidal parts; the chemical power released due to non-equilibrium Urca reactions,  Pol  / H (solid) and  Tor  / H (dashed); and (d) toroidal magnetic fluxes, where the solid lines are the total fluxes and the dashed-dotted are the fluxes computed only in the regions of closed poloidal field lines.The vertical lines correspond, from left to right, to the time scale    (dashed-dotted), and   (dashed), respectively.

Figure 10 .
Figure10.Physical time variable  (in years) for evolving temperature as a function of auxiliary time  ′ (in code units) for constant temperature, as given by equation (78), for two different initial rms magnetic field strengths:   = 2 × 10 16 G (blue) and   = 5 × 10 15 G (red).The dashed vertical line marks the time scale  ′   , which is independent of   .The horizontal lines show, from bottom to top, the time scales    =  ( ′   ) (dashed) mapped from equation (79) for the different initial rms magnetic field strengths with their respective colors, and the transition time    from the strongcoupling to the weak-coupling regime (black-solid).The purple-hatched area displays the time in which GS-equilibrium is reached in our simulations at constant temperature, and the green-hatched area shows the time range corresponding to the weak-coupling regime.

Table 1 .
Summary of the code units, showing the normalization of the different variables, their notations and physical values obtained from the HHJ EoS, for a NS that has a total radius  = 12.2 km, a core radius   = 11.2 km, and a total mass  = 1.4 M ⊙ .The definitions of   ,   ,   , and  *  are given in § 3.4.Here, the free parameters are the initial rms magnetic field strength,   , and the temperature, , which are needed to recover the physical units from the simulation output.The notation  ,15 means   /10 15 G.