  At the high temperatures present inside recently formed neutron stars
($T\gtrsim 5\times 10^{8}\, \text{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 ($\beta$) equilibrium, so
$\beta$ 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 $\gtrsim 10^{16} \,
\mathrm{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 ($T\lesssim 5\times 10^8 \,
\mathrm{K}$), in which Urca reactions are effectively frozen and ambipolar
diffusion becomes relevant.
