ABSTRACT At the high temperatures 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 these states, the particle species deviate from chemical (β) equilibrium, which tends to be restored by β decays (Urca reactions), 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. Here, we present a numerical scheme that decouples the magnetic and thermal evolution, enabling to efficiently perform, for the first time, long-term magnetothermal 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 state is not attained. Therefore, any substantial evolution of the core magnetic field must occur later, in the ‘weak-coupling’ regime ($T\lesssim 5\times 10^8 \, \mathrm{K}$), when Urca reactions are frozen and ambipolar diffusion becomes relevant.