Boson stars are gravitationally bound objects that arise in ultralight dark
matter models and form in the centers of galactic halos or axion miniclusters.
We systematically study the excitations of a boson star, taking into account
the mixing between positive and negative frequencies introduced by gravity. We
show that the spectrum contains zero-energy modes in the monopole and dipole
sectors resulting from spontaneous symmetry breaking by the boson star
background. We analyze the general properties of the eigenmodes and derive
their orthogonality and completeness conditions which have non-standard form
due to the positive-negative frequency mixing. The eigenvalue problem is solved
numerically for the first few energy levels in different multipole sectors and
the results are compared to the solutions of the Schrödinger equation in
fixed boson star gravitational potential. The two solutions differ
significantly for the lowest modes, but get close for higher levels. We further
confirm the normal mode spectrum in 3D wave simulations where we inject
perturbations with different multipoles. As an application of the normal mode
solutions, we compute the matrix element entering the evaporation rate of a
boson star immersed in a hot axion gas. The computation combines the use of
exact wavefunctions for the low-lying bound states and of the Schrödinger
approximation for the high-energy excitations.