  Taking a comprehensive view, including a full range of boundary conditions,
we reexamine QCD axion star solutions based on the relativistic Klein-Gordon
equation (using the Ruffini-Bonazzola approach) and its non-relativistic limit,
the Gross-Pitaevskii equation. A single free parameter, conveniently chosen as
the central value of the wavefunction of the axion star, or alternatively the
chemical potential with range $-m<\mu< 0$ (where $m$ is the axion mass),
uniquely determines a spherically-symmetric ground state solution, the axion
condensate. We clarify how the interplay of various terms of the Klein-Gordon
equation determines the properties of solutions in three separate regions: the
structurally stable (corresponding to a local energy minimum) dilute and dense
regions, and the intermediate, structurally unstable transition region. From
the Klein-Gordon equation, one can derive alternative equations of motion
including the Gross-Pitaevskii and Sine-Gordon equations, which have been used
previously to describe axion stars in the dense region. In this work, we
clarify precisely how and why such methods break down as the binding energy
increases, emphasizing the necessity of using the full relativistic
Klein-Gordon approach. Finally, we point out that, even after including
perturbative axion number violating corrections, solutions to the equations of
motion, which assume approximate conservation of axion number, break down
completely in the regime with strong binding energy, where the magnitude of the
chemical potential approaches the axion mass.
