field structure in the interior depends completely on how currents are internally distributed.

In the astrophysical scenario we are considering, the magnetar reaches the initial state in which we start our numerical simulation after a slow magnetothermal evolution that proceeds in a long timescale compared to the dynamical timescales (cf. section \[sec:timecales\]) of the magnetosphere ($\sim 1$ms) or the crust ($\sim 10$ms). On such long timescales, any current close to the surface of the NS is expected to be dissipated by Ohmic diffusion. Therefore, we consider that initially all fields are continuous across the surface. We build our interior solution by extrapolating the exterior magnetic field towards the stellar interior across a number of grid cells as needed by the reconstruction algorithm used for the magnetospheric evolution in our simulations. Since the neutron star is basically a perfect conductor, the initial charge density and electric field in the interior (and the magnetosphere) are set to zero.

The surface values of $B^r$ and $B^\varphi$ are coincident for degenerate models (e.g. within the series C1, C2 and C3 in Figure \[fig:initial\_model\]) because $P$ and $T$ at the surface are identical. However, since $P$ and $T$ may have a different radial dependence outside 