  Axisymmetric magnetohydrodynamics (MHD) can be invoked for describing
astrophysical magnetized flows and formulated to model stellar magnetospheres
including main sequence stars (e.g. the Sun), compact stellar objects [e.g.
magnetic white dwarfs (MWDs), radio pulsars, anomalous X-ray pulsars (AXPs),
magnetars, isolated neutron stars etc.], and planets as a major step forward
towards a full three-dimensional model construction. Using powerful and
reliable numerical solvers based on two distinct finite-difference method (FDM)
and finite-element method (FEM) schemes of algorithm, we examine axisymmetric
steady-state or stationary MHD models in Throumoulopoulos & Tasso (2001),
finding that their separable semi-analytic nonlinear solutions are actually not
unique given their specific selection of several free functionals and chosen
boundary conditions. The multiplicity of nonlinear steady MHD solutions gives
rise to differences in the total energies contained in the magnetic fields and
flow velocity fields as well as in the asymptotic behaviours approaching
infinity, which may in turn explain why numerical solvers tend to converge to a
nonlinear solution with a lower energy than the corresponding separable
semi-analytic one. By properly adjusting model parameters, we invoke
semi-analytic and numerical solutions to describe different kind of scenarios,
including nearly parallel case and the situation in which the misalignment
between the fluid flow and magnetic field is considerable. We propose that
these MHD models are capable of describing the magnetospheres of MWDs as
examples of applications with moderate conditions (including magnetic field)
where the typical values of several important parameters are consistent with
observations. We discuss the challenges of developing numerical extrapolation
MHD codes in view of nonlinear degeneracy.
