(\[eq:toroidal\]) we solve the GSE and obtain a twisted magnetospheric initial model. We would like to note that all equations can be rescaled with $B_{\rm pole}$, hence, the results of our numerical simulations can be normalised to the field strength of interest.

The energy stored in the magnetosphere can be computed as a volume integral $$\begin{aligned}
 {\mathcal{E}}= \frac{1}{8\pi}\int ({\boldsymbol{B}}^2 + {\boldsymbol{E}}^2)\, {\rm d}V .
 \end{aligned}$$ For later reference and in order to normalize the energetic content of our models, we provide the energy stored in the magnetosphere of a pure dipolar magnetic field ($\vec{E}=0$, $B^r=B_{\rm pole}(R_*/r)^3\cos{\theta}$, $B^\theta=(B_{\rm pole}/2)(R_*/r)^3\sin{\theta}$, $B^\varphi=0$): $$\begin{aligned}
\begin{split}
 {\mathcal{E}}_{\rm d} = \frac{1}{12}B_{\rm pole}^2 R_*^3 =8.3\times 10^{46}\,\text{erg}\,\left(\frac{B_{\rm pole}}{10^{15}\,\text{G}}\right)^2\left(\frac{R_*}{10\,\text{km}}\right)^3.
\end{split}
 \label{eq:edipole}
 \end{aligned}$$ Once the surface value of $P$ and the functional relation $T(P)$ are defined, one can solve the GSE iteratively (as it is a non-linear equation), while imposing vacuum boundary conditions at large distances. We use the numerical code described in @Akgun2018a to build our initial data. Using this parametrisation, the boundary condition at the surface of the neutron star for the GSE (values of $P$ and $T$) is fully determined by four parameters $B_{\rm pole}$, $s$, ${P_{\rm c}}$ and $\sigma$. However, the solution of the GSE with 