domain is defined as $0\leq\theta\leq\frac{\pi}{2}$ in the angular direction and $0\leq \hat{r}\leq 2$ in the radial direction. Since the equation of magnetohydrostationary equilibrium is defined only inside the star and the source terms of the elliptic equations for the gravitational potential and the magnetic flux function vanish outside the star, our computational domain covers a region of the space that is sufficient for obtaining equilibria. In order to resolve the region near the axis sufficiently, we use a special coordinate in actual numerical computations. Total mesh numbers in $r$-direction and in $\theta$-direction are 1025 and 1025, respectively. We describe details of the computational grid points in Appendix \[App:grid\].

Numerical method
----------------

We use the scheme of [@Tomimura_Eriguchi_2005]. This scheme is based on the Hachisu Self-Consistent Field (HSCF) scheme (@Hachisu_1986), which is the method for obtaining equilibrium configurations of rotating stars. We define the ratio of the equatorial radius to the polar radius as the axis ratio $q$. This quantity $q$ characterizes how distorted the stars are due to non-spherical forces. The stronger the non-spherical force becomes, the more distorted the stellar shape is. The non-spherical force can be the centrifugal force, the magnetic force or both of them. We fix 