.\end{aligned}$$ Considering the boundary conditions at infinity and the symmetry of the magnetized stars, the solution of the above equation can be expressed as $$\begin{aligned}
   A_\varphi \sin \varphi
             \equiv  \sum_{n=1}^{\infty} A_{\varphi, n} \sin \varphi
             = \sum_{n=1}^{\infty} b_{n,1} r^{-n-1} Y_{n,1}(\theta,\varphi) \ ,
   \label{Eq:bn}\end{aligned}$$ where $Y_{n,1}(\theta,\varphi)$ is the spherical harmonics of degree $n$ and order $m = 1$. The coefficients $b_{n,1}$ correspond to the magnetic multipoles.

Setting for Numerical Computations
----------------------------------

For numerical computations, the physical quantities are transformed into dimensionless ones using the maximum density $\rho_{\max}$, the maximum pressure $p_{\max}$ and the equatorial radius $r_e$ as follows: $$\begin{aligned}
   \hat{r} \equiv \frac{r}{r_e} = \frac{r}{
   \sqrt{\frac{1}{\alpha}\frac{p_{\mathrm{max}}}{4\pi G
   \rho_{\mathrm{max}}^2}}} \ ,\end{aligned}$$ for polytropic configurations and $$\begin{aligned}
   \hat{r} \equiv \frac{r}{r_e} = \frac{r}{
   \sqrt{\frac{1}{\alpha}\frac{8a}{b}  \frac{1}{4\pi G
   \rho_{\mathrm{max}}^2}}} \ ,\end{aligned}$$ for the Fermi gas models, and $$\begin{aligned}
   \hat{\rho} \equiv \frac{\rho}{\rho_\mathrm{max}} \ .\end{aligned}$$ Here $\alpha$ is introduced so as to make the distance from the center to the equatorial surface of the star to be unity. Arbitrary functions are also transformed into dimensionless ones. Quantities with $\hat{}$ are dimensionless. For example, the dimensionless length is $\hat{r}$ and the dimensionless arbitrary functions are $\hat{\mu}$ and $\hat{\kappa}$, respectively. Dimensionless forms of other quantities are collected in Appendix \[App:dimensionless\].

The computational 