Assumptions of axisymmetric and equatorial symmetries as well as rigid rotation are adopted here in order to simplify our investigations.

In a rotating star under a radiative equilibrium, there appears to be meridional flow in special cases. However, we neglect it because the time scale is many orders of magnitude larger than the (magneto)hydrodynamic one (@Tassoul_2000). Also, there is a suggestion that gradual diffusion of the internal magnetic fields drives a meridional flow (@Urpin_Ray_1994). The time scale of this, again, is much larger than the (magneto)hydrodynamic one. Thus it is also neglected.

Under these assumptions, the basic equations are written as follows. The continuity equation is expressed as $$\begin{aligned}
   \nabla \cdot (\rho \Vec{v}) = 0 \ ,\end{aligned}$$ where $\Vec{v}$ is the fluid velocity. The equations of motion in the stationary state are written as: $$\begin{aligned}
   \frac{1}{\rho} \nabla p = -\nabla \phi_g +
   R \Omega^2 \Vec{e}_R+ \frac{1}{\rho} \left( \frac{\Vec{j}}{c} \times
   \Vec{H} \right),
   \label{Eq:eular}\end{aligned}$$ where $\phi_g$, $\Omega$, $\Vec{j}$, $c$ and $\Vec{H}$ are gravitational potential, angular velocity, electric current density, speed of light and magnetic field, respectively. Here we use the cylindrical coordinates $(R, \varphi, z)$ and $\Vec{e}_R$ is the unit vector in the $R$-direction. The gravitational potential satisfies Poisson equation: $$\begin{aligned}
   \Delta 