  Axion stars are hypothetical objects formed of axions, obtained as localized
and coherently oscillating solutions to their classical equation of motion.
Depending on the value of the field amplitude at the core $|\theta_0| \equiv
|\theta(r=0)|$, the equilibrium of the system arises from the balance of the
kinetic pressure and either self-gravity or axion self-interactions. Starting
from a general relativistic framework, we obtain the set of equations
describing the configuration of the axion star, which we solve as a function of
$|\theta_0|$. For small $|\theta_0| \lesssim 1$, we reproduce results
previously obtained in the literature, and we provide arguments for the
stability of such configurations in terms of first principles. We compare
qualitative analytical results with a numerical calculation. For large
amplitudes $|\theta_0| \gtrsim 1$, the axion field probes the full non-harmonic
QCD chiral potential and the axion star enters the {\it dense} branch. Our
numerical solutions show that in this latter regime the axions are
relativistic, and that one should not use a single frequency approximation, as
previously applied in the literature. We employ a multi-harmonic expansion to
solve the relativistic equation for the axion field in the star, and
demonstrate that higher modes cannot be neglected in the dense regime. We
interpret the solutions in the dense regime as pseudo-breathers, and show that
the life-time of such configurations is much smaller than any cosmological time
scale.
