to infinity, thus, violating the requirements of finite magnetic energy. Following the notation of @Akgun2016 [@Akgun2018a], we write the axisymmetric magnetic field in terms of its poloidal and toroidal components: $$\begin{aligned}
 {\boldsymbol{B}} = {\boldsymbol{\nabla}} P \times {\boldsymbol{\nabla}}\varphi + T {\boldsymbol{\nabla}}\varphi \ ,\end{aligned}$$ where $\varphi$ is the azimuthal angle in spherical coordinates. Here, $P$ and $T$ are the poloidal and toroidal stream functions. Expressed in the orthonormal spherical basis corresponding to the coordinates $(r,\theta,\varphi)$, the magnetic field can be explicitly computed from the potentials $P$ and $T$ as $$\begin{aligned}
B^r &= \frac{1}{r^2\sin\theta}\partial_\theta P, \label{eq:Br}\\ 
B^\theta &= - \frac{1}{r \sin\theta} \partial_r P, \label{eq:Btheta}\\
B^\varphi &= \frac{T}{r \sin\theta}. \label{eq:Bphi}\end{aligned}$$ For an axially symmetric force-free field, the functions $T$ and $P$ may be expressed in terms of each other and appear as solutions of the force-free GSE: $$\begin{aligned}
\left[\partial_r^2 + \frac{1 - \mu^2}{r^2}\partial_\mu^2\right] P + T \frac{dT}{dP} = 0,\end{aligned}$$ where $\mu=\cos\theta$. $P$ and $T$ are constant on magnetic surfaces or, equivalently, along magnetic field lines. $P$ is related to the magnetic flux passing through the area centered on the axis and delineated by the magnetic surface. Therefore, its value at the poles is zero and increases towards the equator. The function $T$ is related 