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 |θ₀| ≡ |θ(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 |θ₀|. For small |θ₀| [< over ~] 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 |θ₀| [> over ~] 1, the axion field probes the full non-harmonic QCD chiral potential and the axion star enters the 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.