Not yet reviewed by Pith; the record is open.

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.