Expansion in Higher Harmonics of Boson Stars using a Generalized
  Ruffini-Bonazzola Approach, Part 1: Bound States
  The method pioneered by Ruffini and Bonazzola (RB) to describe boson stars
involves an expansion of the boson field which is linear in creation and
annihilation operators. In the nonrelativistic limit, the equation of motion of
RB is equivalent to the nonlinear Schr\"odinger equation. Further, the RB
expansion constitutes an exact solution to a non-interacting field theory, and
has been used as a reasonable ansatz for an interacting one. In this work, we
show how one can go beyond the RB ansatz towards an exact solution of the
interacting operator Klein-Gordon equation, which can be solved iteratively to
ever higher precision. Our Generalized Ruffini-Bonazzola approach takes into
account contributions from nontrivial harmonic dependence of the wavefunction,
using a sum of terms with energy $k\,E_0$, where $k\geq1$ and $E_0$ is the
chemical potential of a single bound axion. The method critically depends on an
expansion in a parameter $\Delta \equiv \sqrt{1 - E_0{}^2/m^2} < 1$, where $m$
is the mass of the boson. In the case of the axion potential, we calculate
corrections which are relevant for axion stars in the transition or dense
branches of solutions. We find with high precision the local minimum of the
mass, $M_{min}\approx 463\,f^2/m$, at $\Delta\approx0.27$, where $f$ is the
axion decay constant. This point marks the crossover from the transition branch
to the dense branch of solutions, and a corresponding crossover from structural
instability to stability.
