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(vertical bar psi vertical bar(2)) involving an attractive vertical bar psi vertical bar(4) term and a repulsive vertical bar psi vertical bar(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 <= 0. We show the existence of a critical point vertical bar lambda vertical bar(c) similar to (m/M-p)(2), where m is the axion mass and M-p is the Planck mass, 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 vertical bar lambda vertical bar < vertical bar lambda vertical bar(c), a system of self-gravitating axions forms a stable dilute axion star below a general relativistic maximum mass M-max,GR(dilhute) similar to M-P(2)/m and collapses into a black hole above that mass. For strong self-interactions vertical bar lambda vertical bar < vertical bar lambda vertical bar(c), a system of self-gravitating axions forms a stable dilute axion star below a Newtonian maximum mass M-max,N(dilhute) = 5.073M(p)/root vertical bar lambda vertical bar [Phys. Rev. D 84, 043531 (2011)], collapses into a dense axion star above that mass, and collapses into a black hole above a general relativistic maximum mass M-max,GR(dence) similar to root vertical bar lambda vertical bar M-P(3)/m(2). Dense axion stars explode below a Newtonian minimum mass M-min,N(dense) = 98.9m/ root vertical bar lambda vertical bar 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 (vertical bar lambda vertical bar*, 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. Our approximate analytical results are in good agreement with the exact numerical results of Braaten et al. [Phys. Rev. Lett. 117, 121801 (2016)] for Newtonian dense axion stars. They are also qualitatively similar to those obtained by Helfer et al. [J. Cos mol. Astropart. Phys. 03 (2017) 055] for general relativistic axion stars, but they differ quantitatively for weak self-interactions presumably due to the use of a different self-interaction potential V(vertical bar psi vertical bar(2) ).We point out analogies between the evolution of self-gravitating axions (bosons) at zero temperature evolving from dilute axion stars to dense axion stars and black holes and the evolution of compact degenerate (fermion) stars at zero temperature evolving from white dwarfs to neutron stars and black holes. We also discuss some analogies between the phase transitions of Newtonian axion stars at zero temperature and the phase transitions of Newtonian self gravitating fermions at nonzero temperature. Finally, we suggest that a dense axionic nucleus may form at the center of dark matter halos through the collapse of a dilute axionic core (soliton) passing above the maximum mass M-max,N(dilute). It would have a mass 1. 11 x 10(9)(f/m)M-circle dot, a radius 0.949/(mf(1/3)) pc, a density 2.10 x 10(-8) (m(2)f(2)) g/m(3) a pulsation period 8.24/(mf(1/3)) yr, and an energy 5.59 x 10(62)(f/m) erg, where the axion mass In is measured in units of 10(-22) eV/c(2) and the axion decay constant f is measured in units of 10(15) GeV. This dense axionic nucleus could be the remnant of a bosenova associated with the emission of a characteristic radiation
