  We study the nature of phase transitions between dilute and dense axion stars
interpreted as self-gravitating Bose-Einstein condensates. We develop a
Newtonian model based on the Gross-Pitaevskii-Poisson equations for a complex
scalar field with a self-interaction potential $V(|\psi|^2)$ involving an
attractive $|\psi|^4$ term and a repulsive $|\psi|^6$ term. Using a Gaussian
ansatz for the wave function, we analytically obtain the mass-radius relation
of dilute and dense axion stars for arbitrary values of the self-interaction
parameter $\lambda\le 0$. We show the existence of a critical point
$|\lambda|_c\sim (m/M_P)^2$ above which a first order phase transition takes
place. We qualitatively estimate general relativistic corrections on the
mass-radius relation of axion stars. For weak self-interactions
$|\lambda|<|\lambda|_c$, a system of self-gravitating axions forms a stable
dilute axion star below a general relativistic maximum mass $M_{\rm
max,GR}^{\rm dilute}\sim M_P^2/m$ and collapses into a black hole above that
mass. For strong self-interactions $|\lambda|>|\lambda|_c$, a system of
self-gravitating axions forms a stable dilute axion star below a Newtonian
maximum mass $M_{\rm max,N}^{\rm dilute}=5.073 M_P/\sqrt{|\lambda|}$, collapses
into a dense axion star above that mass, and collapses into a black hole above
a general relativistic maximum mass $M_{\rm max,GR}^{\rm dense}\sim
\sqrt{|\lambda|}M_P^3/m^2$. Dense axion stars explode below a Newtonian minimum
mass $M_{\rm min,N}^{\rm dense}\sim m/\sqrt{|\lambda|}$ and form dilute axion
stars of large size or disperse away. We determine the phase diagram of
self-gravitating axions and show the existence of a triple point
$(|\lambda|_*,M_*/(M_P^2/m))$ separating dilute axion stars, dense axion stars,
and black holes. We make numerical applications for QCD axions and ultralight
axions.
