Boson Star Normal Modes

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\"odinger 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\"odinger approximation for the high-energy excitations.

Compared with other DM candidates, ULDM shares a similar quality in explaining the DM observations in the large-scale structure, but has novel properties on a small scale that is about its de Broglie wavelength.As a distinct feature of ULDM, coherent soliton solutions, known as axion solitons or boson stars [17], exist due to the interplay between the self-gravity of the axion field and its wave properies.All the axions in a boson star are in the same ground state supported by the gravitational potential.Because the soliton solutions are energetically favorable, they are formed by the kinematic relaxation of axions through gravitational interactions [18] or other mechanisms.The formation of these objects in the central region of axion dark matter halos (also known as the axion mini-clusters [19,20] for the QCD axion) is confirmed by several numerical simulations in the cosmological setting [21][22][23][24], kinetic relaxation region [18,25] or by considering the collision of several seed solitons [26][27][28].Even if the soliton is artificially removed from the halo, evolution will restore it back [29].The size of the formed solitons is about the de Broglie wavelength of the surrounding ULDM.The ULDM with mass m a ∼ 10 −22 ÷ 10 −19 eV, known as fuzzy DM [30], has the wavelength of many parsecs, so is the size of soliton in the center.Tight constraints have been put on fuzzy DM from Ly-α forest data [31][32][33] and dynamics of dwarf galaxies [34,35].
The observational signatures from axion solitons open new opportunities to discover or constrain axion DM.For the fuzzy DM, the presence of solitons in the center of galaxies can modify the galaxy rotation curves [36,37], and heat the stars in the central region via oscillations and random walk [34,38,39].Also, the presence or absence of axion solitons have important implications for lensing search [40].When considering the self-interactions [41] and couplings to the photons [42][43][44], the axion solitons can produce relativistic axions and radio emissions, respectively.The dense axion stars made of inflaton can seed primordial black holes or produce stochastic high-frequency gravitational wave signals [45].
All these interesting signatures are based on the fact that the existence of axion solitons appears as a generic feature of an axion halo.This generic feature is observed through extensive studies on the formation and evolution of axion solitons [18,[21][22][23][24][25], which rely on the interactions between the solitons and their environment, the surrounding axion waves.Refs.[21,26] find an interesting correlation between the soliton mass and the host halo mass in cosmological numerical simulations, and the relationship is confirmed in [22,46].Although the soliton-host halo relation appears intuitive by equating the virial temperature of the soliton and the host halo [47], the underlying mechanism remains unclear.It is not reproduced by other simulations in non-cosmological setting [27,28,48,49], and more recent cosmological simulations [24,48,50] observe a larger scatter of the relation.In the kinetic regime, the soliton-host halo relation is disfavored because a continuous growth of solitons after their formation has been observed [18], while refs.[46,51] argue that the growth slows down when the soliton becomes heavy enough to heat up the inner part of axion halos.
In ref. [25], by combining the analytic approach and numerical simulations we find that, somewhat counter-intuitively, solitons with virial temperature higher than the temperature of the ambient axion gas grow by condensation from the gas, whereas solitons colder than a certain threshold evaporate.Since the temperature of the soliton is proportional to the square of its mass, the former solitons were dubbed in [25] as 'heavy', and the latter ones as 'light'.From the analytic computations, it is evident that the dynamics of solitons arise from the scattering between a soliton and its surrounding axion waves.To explore the dynamics of solitons, we need to understand the eigenstates of axion waves in the soliton background.
The complexity of the axion modes occurs because the modes with positive and negative frequencies with respect to the soliton background are coupled together by gravity through the Schrödinger-Poisson (SP) system of equations.Two approaches to this problem exist in the literature.Refs.[29,49,52] neglect the mixing of modes so that the equation is simplified to a linear Schrödinger equation in fixed gravitational potential of the soliton and the halo.This approximation seems to work well for studies of collective properties of the axion eigenstates.However, it fails in some important details -in particular, it does not yield a dipole zeroenergy mode whose presence is required by the (spontaneously broken) translation invariance.On the other hand, refs.[53][54][55] find numerically several low-lying excited states, for which the mode mixing plays an essential role.
The goal of the paper is to connect the two approaches and to clarify whether and when we can neglect mode mixing.Our analysis follows the lines of refs.[54,55], complementing them in several ways.We lay out a general framework to find all eigenstates and derive their orthogonality and completeness properties.We find numerically the energies and wavefunctions of all low-lying sates with + n ≤ 4, where is the orbital (multipole) quan-tum number and n ≥ 0 is the level number within the given multipole sector. 1 We confirm that the positive-negative mode mixing is important for accurate determination of the lowest energy levels and show how it leads to recovery of the zero-energy dipole mode.On the other hand, we observe that for high-frequency excitations corresponding to + n ≥ 4 the solutions of the SP system are well approximated by the eigenfunctions of the Schrödinger equation in fixed gravitational potential.We further improve the approximation by developing a perturbative expansion for exact eigenvalues around the eigenvalues of the Schrödinger equation.These results are compared with spectroscopic measurements of a perturbed soliton in (3 + 1)-dimensional simulations and are found to be in excellent agreement.
The complete axion wavefunctions and their high-frequency approximation are crucial information to deduce the scattering rate of axion waves on the soliton.We illustrate this point on the example of the light soliton evaporation rate [25] which we evaluate by combining the advantages of the exact and the Schrödinger wavefunctions.
The paper is organized as follows.In section 2, we introduce the SP equations for soliton background and perturbations.We analyze the general properties of the eigenfunctions in section 3. Then we numerically solve the SP equations and Schrödinger equations to obtain the exact and approximate wavefunctions in section 4. Soliton spectroscopy using 3D wave simulation is presented in section 5.In section 6 we explore the interaction between the soliton and axion waves and evaluate the evaporation rate of a light soliton.We conclude in section 7. Appendix contains some details about our numerical procedure to solve the eigenvalue problem.

Soliton and its Perturbations
We consider a non-relativistic complex scalar field ψ(t, x) only with gravitational self-interactions.Its dynamics follows the Schrödinger-Poisson (SP) equations, ) where ∆ denotes the Laplacian.The density of scalar field is ρ(t, x) = m|ψ(t, x)| 2 , and Φ is the gravitational potential given by the density.The system has scaling symmetry which can be used to generate a family of solutions from a single solution.More general scaling transformations given in [25] allow one to eliminate m and 4πG and work with dimensionless quantities.However, we prefer to keep m and 4πG explicitly in the analytic derivations and only switch to the units m = 4πG = 1 in the numerical calculations.We assume spherical symmetry for soliton solutions, then the Ansatz for a soliton with the binding energy E s (E s < 0) takes the form The SP equations reduce to time-independent ordinary differential equations for χ and Φ χ .We impose the regular boundary conditions of χ | r=0 = 0 and Φ χ | r=0 = 0 at the center, apply the scaling symmetry of the equations to normalize χ| r=0 = 1, and look for solutions with the asymptotics χ| r→∞ = Φ χ | r→∞ = 0.The numerical solution with the number of nodes equal to 0, as appropriate for the ground state, provides us with a "standard soliton", whose profile and gravitational potential are shown in fig. 1.
The standard soliton can be rescaled to other soliton solutions through the scaling symmetry (2.2).Restoring the dimensionful variables, a general soliton reads, where k −1 s sets its size.The scaling of coordinates and wavefunctions leads to the scaling of other quantities: the density ρ s ∝ S 4 , the total soliton mass M s = d 3 x ρ s ∝ S, and the binding energy E s ∝ S 2 .In more detail, we have where µ 0 = 25.9 and ε 0 = −0.692are the dimensionless mass and binding energy of the standard soliton.
Let us now treat the soliton as background and consider small perturbations on top of it.It is convenient to factor out the soliton time dependence and write the axion field as (2.6) We obtain the SP equations for the perturbation by keeping only linear terms in δψ(t, x).
Integrating out the perturbation of the gravitational potential yields a non-local equation, Equation (2.7) can be conveniently written in the following matrix form, and (2.9) Note that the soliton background satisfies the equation From the above equations, it is clear that δψ and δψ * are coupled at linear order via gravity.As a consequence, δψ and δψ * will contain both positive and negative frequency components.The situation is similar to the well-known mode mixing in atomic Bose-Einstein condensates which implies a need for a non-trivial Bogolyubov transformation to diagonalize the Hamiltonian [56].
If we neglect V in eq.(2.8), it reduces to the Schrödinger equation i∂ t δψ (Sch) = H 0 δψ (Sch) , or explicitly: This does not have mixing between positive and negative frequency modes.Neglecting V is in general not justified since inside the soliton V and mΦ χ are of the same order.The low-lying bound states whose wavefunctions significantly overlap with the soliton will be affected comparably by these two terms.On the other hand, one expects the Schrödinger approximation to become more accurate for highly excited states with wavefunctions wider than the soliton.Indeed, the Newtonian potential has a long tail, Φ χ = −GM s /r outside the soliton, whereas the operator V is proportional to the soliton profile χ and hence is sharply localized at distances r k −1 s .Therefore, V is expected to have weak effect on the dynamics of high-level modes.Our numerical results below will confirm these expectations.

General Properties of Eigenstates
In this section we study the coupled positive and negative frequency modes governed by the linearized SP equations and discuss their general properties, including the orthogonality and completeness relations.

Positive, negative, and zero frequency modes
Equations (2.8) and (2.9) show that the positive and negative frequency modes are coupled.This coupling guides us to look for eigenmodes of the form, Here the subscript a labels different energy levels and a is the energy of the excitation relative to the soliton ground state.The excitations are bound to the soliton as long as their energy is below the soliton binding energy, a < |E s |; otherwise they describe a continuum spectrum of axion waves propagating to infinity.Due to the spherical symmetry of the soliton background, a is a set of orbital, radial, and magnetic quantum number, a = ( , n, m), as we factorize ϕ into the radial and angular components, Here Y m are the spherical harmonics. 2The non-negative integers n ≥ 0 label the discrete energy levels at a given orbital quantum number . 3 For the continuum spectrum, we can choose instead of n the absolute value of the mode momentum at infinity.The SP equations (2.8) are converted to coupled equations of ϕ a and ϕ a , These can be cast in a compact form, and σ 3 = diag(1, −1) being the third Pauli matrix.In this way we have reformulated the problem of finding the spectrum of excitations in the presence of the soliton as an eigenvalue problem in the space of two-component columns ϕ, with ϕ (1) corresponding to positive and ϕ (2) to negative frequencies.Note, however, that the effective Hamiltonian σ 3 H is neither Hermitian, nor positive definite.Instead, its spectrum is symmetric around zero.To see this observe that H commutes with the first Pauli matrix, [σ 1 , H] = 0.Then, given a state ϕ a with positive energy, we can construct its dual with the negative eigenvalue, σ 3 H φa = − a φa .Clearly, the duality interchanges the positive and negative frequency components.Comparing with eq.(3.1), we see that the dual state corresponds to the complex conjugate axion wavefunction.
The Hamiltonian also has zero-energy eigenstates which deserve a separate discussion.These are the Goldstone modes associated to the symmetries of the axion dynamics spontaneously broken by the soliton.It is straightforward to identify these symmetries from the axion action where we have integrated out the gravitational potential Φ.It is invariant under the U (1) phase shifts of the complex scalar field, as well as under spatial translations.Acting with these symmetries on the soliton we obtain other solutions of the field equations with the same energy, The Goldstone modes are the infinitesimal versions of these transformations giving zero-energy states in the monopole and dipole sectors, which in the two-component form read, Recall that the first number in the subscript stands for the angular momentum = 0 or = 1, and the second number 0 stands for the ground state in the corresponding sector.Note that, with a slight abuse of notations, we have denoted by ϕ 10 (x) a triplet of modes differing by the magnetic quantum number m = −1, 0, 1; this triplet can also be written as It is instructive to check that the modes (3.8) indeed satisfy eq.(3.4) with zero r.h.s. by a direct substitution.For the monopole ϕ 00 the terms with V cancel and we arrive at the background equation (2.10).For the dipole, we act with the derivative ∇ on eq.(2.10) and obtain (H 0 + 2V )∇χ = 0, which is precisely the equation given by substituting ϕ 10 into eq.(3.4).This reasoning also tells us that the Schrödinger equation (2.11) has the monopole zero mode, but not the dipole zero mode, because it neglects V and H 0 ∇χ = 0.

Orthogonality and completeness relations
Since the operator in the eigenvalue problem (3.4) is not Hermitian, it leads to unusual orthogonality and completeness relations.By using eq.(3.4) and the identity a Hϕ b ) † = 0, we find the first orthogonality condition Assuming a = b we can write the condition in terms of positive and negative frequency components, Note the minus sign in front of the negative-frequency contribution.Likewise, we use eq.(3.4) and another identity d 3 x (ϕ T b σ 1 Hϕ a ) − (ϕ T a σ 1 Hϕ b ) T = 0 to obtain the second orthogonality condition, In terms of components, it reads Note that it is also satisfied for the zero modes (3.8), despite the fact that for them a = b = 0.
The first orthogonality condition (3.9) implies that we can set the inner product and normalization of positive-energy wavefunctions by using σ 3 , This raises an issue of the completeness of the eigenmode basis.Indeed, consider e.g. a column ϕ = (χ(r), 0) T .If it were possible to decompose it as a linear combination of eigenmodes, its product with ϕ 00 would be zero.However, the explicit evaluation gives ϕ, ϕ 00 = d 3 x χ 2 (r) > 0 . (3.17) Thus the expansion (3.16) must be supplemented by a basis state having non-zero inner product with ϕ 00 .We will refer to this vector as conjugate to ϕ 00 and will denote it by φ00 .
Repeating the argument for ϕ = (∇χ(r), 0) T we conclude that we also need to include a vector φ10 conjugate to ϕ 10 .
To find φ00 , φ10 we again use the symmetries of the problem.Recall that the SP system is invariant under the scaling (2.2) which is broken by the soliton background.The scaling transformation changes the binding energy of the soliton and hence does not give rise to a Goldstone mode.Instead, it provides us with φ00 .To see this, apply an infinitesimal transformation S = 1 + α, α 1, to the soliton.This generates a perturbation, Note that it linearly grows with time.Since both the rescaled and the original solitons satisfy the SP system, the perturbation (3.18) obeys the dynamical equation (2.8).Substituting and using Hϕ 00 = 0 we obtain, From (σ 3 H) 2 φ00 = σ 3 H ϕ 00 = 0 and (σ 3 H) 2 ϕ a = 2 a ϕ a , we infer that φ00 is orthogonal to all eigenmodes with strictly positive energies and their duals.Clearly, its inner product with itself also vanishes.On the other hand, φ00 has non-zero inner product with ϕ 00 , and thus is the sought-after conjugate state to ϕ 00 .Note that eq.(3.19) implies that in the subspace spanned by ϕ 00 and φ00 the effective Hamiltonian σ 3 H has the form of a Jordan cell ( 0 1 0 0 ).To summarize, the orthogonality relations for the sates in the (00) sector are One more symmetry spontaneously broken by the soliton is the Galilean invariance, under which the axion field changes as Applying it to the soliton, we find a linearly growing perturbation, Substituting into eq.(2.8) we arrive at We see that in the (10) sector the effective Hamiltonian σ 3 H is a also Jordan cell.Reasoning as before, we conclude that the state φ10 has zero norm and is orthogonality to all eigenstates of σ 3 H, except ϕ 10 for which we have, where we have recalled that ϕ 10 and φ10 are actually triplets of states.As a summary, the orthogonal conditions for the (10) sector are We have now constructed the complete basis in the space of two-component columns.It consists of the eigenmodes of σ 3 H and the two conjugate states φ00 and φ10 .An arbitrary column can be represented as and the coefficients in this expansion are obtained by taking the inner product of ϕ with the basis states: This implies the following decomposition of the δ-function: Returning to the physical axion perturbations δψ, we observe from eq. (2.8) that they correspond to self-dual columns.This imposes restrictions on the coefficients in eq.(3.27): αa must be complex conjugate of α a , whereas α 00 , α10 must be purely imaginary and α00 , α 10 real.Restoring the time dependence, we obtain the most general axion perturbation in the soliton background, with Re α 00 = Re α10i = 0, Im α00 = Im α 10i = 0. Finally, let us note that for the eigenmodes of the Schrödinger equation (2.11) the orthogonality and completeness relations are standard, so that any solution of eq.(2.11) is expanded as The energies here are the eigenvalues satisfying and are in general different from a .

Energy Eigenstates
In this section we present the numerical results for the energy eigenstates in the soliton background and compare the levels and wavefunctions of the full SP system with those of the Schrödinger equation in the fixed gravitational potential.We work in dimensionless units m = 4πG = 1 and choose the standard soliton discussed in section 2 as the background.
Solutions for any other soliton can be obtained by the rescaling (2.2).

Numerical results for the Schrödinger-Poisson and Schrödinger equations
We solve numerically the eigenvalue problem (3.3) in various multipole sectors labeled by the orbital number .To resolve the non-locality of the operator V (see eq. (2.9)), we integrate back in the perturbation of the gravitational potential which obeys the equation ∆δΦ = χ (δψ + δψ * ).Decomposing it into spherical harmonics similarly to eq. (3.2), we arrive at the system of 3 coupled second-order ODEs: Here χ 0 , Φ 0 and ε 0 are the profile, gravitational potential and the binding energy of the standard soliton.
For the Schrödinger equation (3.33) we obtain a single ODE, We solve the system (4.1) and eq.(4.2) numerically subject to the regularity conditions at the center r = 0 and vanishing conditions at infinity.Under these conditions all wavefunctions can be chosen real.The details of our numerical procedure are described in appendix A. The numerical results for the wavefunctions and the energy levels are summarized in fig. 2 and table 1.We obtain the first few soliton eigenmodes for = 0, 1, 2, 3, 4 and compare them to the Schrödinger solutions.In the last row of the table we also show the results of perturbative improvement of the Schrödinger approximation, to be discussed in the next subsection.Evidently, the positive frequency component of the exact eigenmodes gets closer to the Schrödinger wavefunction for higher levels and large .At the same time, the negative frequency component becomes negligible.This justifies the replacement of the exact SP problem by the Schrödinger approximation with fixed gravitational potential at high energy levels to simplify some computations.Note also the degeneracies between the high energy levels with different and n, but the same sum + n.They stem from the fact that the wavefunctions of highly excited states are much wider than the soliton (see fig. 2) and effectively feel only the Coulomb-like tail of the soliton gravitational potential Φ 0 (r) −µ 0 /(4πr), r 1.For the low energy levels and small , the two solutions of the SP system and the Schrödinger equation are essentially different.For instance, the numerical solutions of the energy levels in table 1 confirm that there are monopole and dipole zero energy modes for the exact solutions, but the Schrödinger equation has only the monopole zero-energy mode.This is consistent with the fact that the Schrödinger equation explicitly breaks the translational symmetry.We have also checked that the exact zero modes of the SP system have the expected form (3.8).
Our results for the eigenfrequencies in the monopole sector agree with the previous studies [53,54] within O(10 −3 ).    1.The energy levels in the background of the standard soliton obtained using the full Schrödinger-Poisson (SP) system and their estimates from the Schrödinger (Sch.)approximation.
The last row shows the Schrödinger eigenvalues improved using perturbation theory described in section 4.2.

Perturbation theory
We have seen above that the eigenvalues of the full SP system and the Schrödinger approximation approach each other for high energy levels.Since the single Schrödinger equation is significantly easier to solve than the SP system, it motivates us to look for the solution of the latter perturbatively, starting from the Schrödinger wavefunction at the zeroth order.This is similar to the perturbation theory in quantum mechanics with H 0 playing the role of the unperturbed Hamiltonian and the terms V in eq.(3.3) treated as the perturbation.
Since the eigenfunctions of H 0 form a complete basis, we can expand the positive and negative frequency components of the exact eigenstate as We insert this expansion into eq.( 3.3), multiply by ϕ (Sch) c * (x), integrate over x and use the orthonormality of the Schrödinger basis.This gives the conditions for the coefficients, with the matrix elements where j (z) are the spherical Bessel functions and we have used the identity By subtracting eq.(4.4b) from eq. (4.4a), we find a simple relation, 5 Inserting the above relation back to eqs.(4.4) we find, 5 This relation does not apply to the monopole ground state, for which 00 = (Sch) 00 = 0.In this case, however, we do not need perturbation theory since the corresponding mode is known exactly, ϕ Until this point, all the formulas have been exact.Let us now assume that the matrix elements V (Sch) cb are small.In the leading order we set them to zero in eqs.(4.9) and obtain that the only solution is ) ac = 0 .(4.10) As expected, at this order the negative frequency component vanishes, while the positive frequency one coincides with the Schrödinger wavefunction.The next step is to take into account the corrections linear in V .Using (4.10) in the terms in eqs.(4.9) that are already proportional to V , we find at the next-to-leading order, The results of the numerical evaluation of the energy levels using eq.(4.11a) are presented in the last row of table 1.They are in excellent agreement with the exact eigenvalues for almost all the modes, except the lowest levels ( = 0, n = 1) and ( = 1, n = 0).Using eqs.(4.11b), (4.11c) to find the wavefunctions is more challenging since the expansions (4.3) involve summation over infinitely many Schrödinger modes.We have found that the sums converge rather slowly and we do not attempt to evaluate them here.

Soliton Spectroscopy with 3D Simulations
In this section we observe how the energy levels of the soliton bound states manifest themselves in the spectrum of soliton perturbations in dynamical wave simulations on a 3-dimensional grid.We inject = 0, 1, 2 perturbations separately in the simulations and perform the frequency analysis of the density field at several points in space.Its power spectrum contains resonant lines corresponding to the soliton eigenmodes with the three orbital angular momenta .The resonant frequencies excellently match the prediction given by the theoretical and numerical studies before.We follow the same approach and conventions as ref.[25] to perform the wave simulations.A soliton wavefunction is introduced in a cubic lattice of size N = 256 with periodic boundary conditions.The axion field is evolved by solving the dimensionless SP equations (we set 4πG = m = 1), (5.1b) We further fix the lattice spacing to unity, dx = 1, by the remaining scaling symmetry of the Schrödinger-Poisson equations.
We use the soliton background shown in fig. 1 as an input for the lattice initial condition, rescaling it by a factor S = 0.428 in order to achieve a better spatial and temporal resolution.The parameters of the soliton are fixed accordingly: the peak density ρ c = 0.0337, the radius at half peak density r 1/2 = 3.03, the total mass M s = 11.1, and the binding energy E s = −0.127.The rescaling also applies to all the energy levels of the soliton perturbations, some of which are listed in table 1: An energy level is rescaled to S 2 n in the simulation setup.
To analyze the soliton spectrum, we follow ref.[55] and introduce a perturbation δψ with a fixed orbital angular momentum of the form The perturbation represents a wavepacket of width σ p in the radial direction and centered at r p .Its angular dependence is a pure spherical harmonic Y m .We choose m = 0, so that the perturbation is axially symmetric (no φ-dependence).In the simulations, we take r p = 20, σ p = 1 and A 0 = 0.01, 0.005, 0.0025 for = 0, 1, 2, which perturbs the soliton mass by no more than ∼ 0.5%.We evolve a soliton and its perturbation on the lattice using the leapfrog integration algorithm [57,58], (5. 3) The time evolution with the operator e i ∆ dt/4 is computed in the momentum space, while the phase e −i Φ(t+dt/2,x) is evaluated in real space.We choose dt = 1/(2π) and the simulation time t max = 80, 000 dt.In order to perform the spectroscopic analysis of the system, we Fourier transform the density ρ(t, x) = |ψ(t, x)| 2 into its frequency space ρ(ω, x).The resolution in frequency space is determined by the total simulation time, dω = π/t max ∼ 4 • 10 −5 .The power spectra for = 0, 1, 2 perturbations P (ω, x) = |ρ(ω, x)| 2 are presented in figs.3,4,5, in which the four panels show the power spectra measured at x = (0, 0, 0), (5, 0, 0), (0, 5, 0) and (0, 0, 5) in Cartesian coordinates.The visible resonant peaks are consistent with energy eigenvalues predicted in the previous study.The summary and analysis of the three perturbations are given in order: = 0 perturbation: In fig.3, the vertical dotted red lines correspond to the predicted energy levels of = 0 modes with n = 1 to 4 (see table 1), rescaled to the simulation units.In each panel (r = 0, x = 5, y = 5 and z = 5), there are four visible peaks, though the fourth one is hardly seen.The locations of the peaks match well with the vertical red lines.The fundamental n = 1 mode is prominent in the soliton center (r = 0), but gets weak at the distance r = 5 because its wavefunction quickly decreases outside the soliton.Also, it is less excited than the n = 2 mode since the injected perturbation is localized far away from the center.As a consequence, the most prominent mode at x = 5, y = 5 and z = 5 is the second mode n = 2.Note that the amplitudes of the n = 2 peaks are similar in all directions which is consistent with the monopole nature of the perturbation.

= 1 perturbation:
In fig.4, the vertical dotted blue lines correspond to the predicted energy levels of = 1 modes and n = 1 to 3. We see a prominent peak corresponding to the ( , n) = (1, 1) mode in the z = 5 panel, together with much smaller peaks corresponding to ( , n) = (1, 2) and (1, 3) modes.The peaks are essentially absent in all other panels -notice the much smaller vertical scale in these panels.This is consistent with the structure of the injected perturbation which is proportional to Y 1,0 (θ) and vanishes in the z = 0 plane.
At r = 0 and x = 5, y = 5 we still observe some weak resonant lines corresponding to = 0 and = 2 modes.This may be due to the mode mixing caused by the discretization and boundary effects from the lattice.
Ref. [55] reports a resonant peak in the dipole sector at very low frequency, about 20 times lower than the frequency of the ( , n) = (1, 1) mode.We do not see such peak in our simulations.We believe that its appearance is due to an explicit breaking of the translation invariance by the code used in [55] which discretizes the equations in cylindrical coordinates and thus introduces a preferred origin at the level of the lattice.As a results, the zero dipole mode associated to this symmetry acquires a small non-zero frequency.Our code, on the other hand, preserves exact lattice translation invariance thanks to the use of Cartesian coordinates and periodic boundary conditions, so the zero frequency of the translational mode remains protected.It is worth stating, however, that in a physical situation of a soliton embedded in a dark matter halo, the translation invariance is expected to be broken by the gravitational potential of the halo, so that the translational mode will have non-zero frequency.The precise value of the frequency will depend on the halo density and its profile.= 2 perturbation: In fig. 5, the vertical dotted green lines correspond to the predicted energy levels of = 2 modes and n = 0 to 2. There are three prominent peaks in the density power spectrum measured away from the soliton center at the points on the x, y and z axes.They are due to the quadrupole perturbation proportional to Y 2,0 (θ).Since Y 2,0 (0) = −2Y 2,0 (π/2), the power spectrum in the z direction is four times larger than the power spectrum in the x and y directions, as observed in these panels.There is no resonance corresponding to any of the quadrupolar modes at the origin since their wavefunctions vanish at r = 0. Due to the mode mixing induced by the lattice, the fundamental monopole line is present in the first panel (r = 0) but with a much smaller power spectrum.
It is truly satisfying to see that the study of = 0, 1, 2 perturbations in the wave simulation confirms the energy eigenvalues deduced from solving the Schrödinger-Poisson equations.

An Application: Light Soliton Evaporation
The soliton bound states play an important role in the dynamics of a soliton immersed in the gas of axion waves [25].Scattering of a wave on the soliton can lead to excitation/deexcitation of such states which reduces/increases the number of particles in the soliton itself.In this section we use the bound state wavefunctions found above to calculate the numerical coefficient in the evaporation rate of light (cold) solitons discussed in [25].It will be clear from this example that the use of the exact wavefunctions is essential to obtain a finite result: The Schrödinger approximation would lead to infrared divergences.The example will also illustrate the efficiency of using the Schrödinger approximation for the high-level and continuum modes to simplify the derivations.

Review of light soliton and axion gas scattering
We start by briefly reviewing the derivation of the light soliton evaporation rate presented in [25].We consider a process when an axion from the gas (g) scatters on the soliton (s) and kicks out another axion into the gas.The soliton evaporation rate is given by the difference of the rates of this process and its inverse, Importantly, we include into the soliton only the particles in the ground state ( = 0, n = 0), whereas all excited levels, including the discrete ones are ascribed to the gas.This division is justified because transferring a particle from the ground state to an excited state, even if the latter is bound, depletes the coherent condensate in the soliton.The soliton-gas scattering is described by the non-linear term in the axion action 7 which contains a cubic contribution in the gas field, ) where the soliton wavefunction has the form (2.4).The gas field here should be understood as the second-quantized operator of perturbations in the soliton background.It is expanded in the eigenmodes as δψ g (t, x) = a>0 a a ϕ (1)  a (x) e −i at + a † a ϕ (2)   a * (x) e i at , ( where a a and a † a are annihilation and creation operators for an axion in the ath state.Note that we do not include in the expansion the zero modes or their conjugates since they correspond to coherent scaling of the soliton wavefunction or motion of the soliton as a whole and thus do not represent genuine soliton excitations.Inserting eq. ( 6.3) into the action we obtain, With this action one can compute the rates of the processes (6.1) using the standard quantummechanical perturbation theory.Note that the term in the second line leads to an s-channel exchange diagram.This term is generally subdominant compared to the contribution from the first line.For example, it is small if the particle a belongs to the continuum since then its negative-frequency wavefunction ϕ a is negligible.On the other hand, in the light soliton case considered here, it is suppressed by the hard momentum transfer, as will be clear shortly.
We assume that the occupation numbers of the unbounded gas modes are given by the Maxwell distribution, where k is the momentum of the mode far away from the soliton.The case of light/cold soliton corresponds to the setup when the characteristic gas momentum is much larger than the inverse soliton size, In this case the dominant contribution into the rates comes from the first term in eq. ( 6.4) and its hermitian conjugate describing the processes c + s ↔ a + b where the particles b and c are 'hard', i.e. have momenta of order k g , but the momentum transfer is 'soft', i.e. the difference between their momenta, is only of order k s .This 'collinear' scattering is enhanced compared to other channels due to the presence of inverse Laplacian in the vertex which makes the amplitude inversely proportional to the square of the momentum transfer.For hard particles of the continuous spectrum we can neglect the negative-frequency component of the wavefunctions and set ϕ = 0. Furthermore, their positive-frequency components can be chosen to be simple plane waves, ϕ c = e ikcx .On the other hand, we cannot drop the negative-frequency part for the particle a, which has momentum of order k s and can be in an unbound or bound state.
With these considerations, and including the Bose enhancement factors for initial and final states we get the rate of the light soliton mass change, where f a,b,c are the occupation numbers in the respective states and A as,bc is the matrix element following from the action (6.4) stripped off the energy conservation δ-function, A as,bc = (4πGm 2 ) Here the vertex form factors are Due to the collinear scattering kinematics we have ) which leads to the cancellation of the first two terms in brackets in eq. ( 6.7) at the leading order in ν.In this way the occupation number of the soft state a drops out of the rate.The remaining rate happens to be negative implying that the mass of the soliton decreases, i.e. it evaporates.
To proceed, it is convenient to use the characteristic momentum of particles in the soliton k s as a unit to rescale other quantities: the coordinate, momentum, energy, and wavefunction By converting all the variables to the rescaled ones, the scattering rate takes the form where and the rescaled form factor is All quantities in eqs.(6.12), (6.13) are dimensionless and refer to the standard soliton background χ 0 (ξ).In particular, µ 0 is the standard soliton mass introduced in eq.(2.5).Taking the integral over the average momentum (κ b + κ c )/2 we obtain at the leading order in ν, (6.14) We now turn to the evaluation of the coefficient C ls .

Scattering rate coefficient
To find C ls , we need to sum over all the possible states of particle a.This task may appear prohibitively complicated since it requires knowledge of all eigenfunctions in the soliton background entering into the form factor (6.13).We overcome this difficulty by the following strategy.We use the exact eigenfunctions for the modes of the few lowest levels, but the higher modes are replaced by the Schrödinger approximation.This allows us to evaluate the infinite sum over states by using the sum rules for the Schrödinger wavefunctions.Note that it is impossible to directly apply the completeness relations (3.29) for the exact wavefunctions since they do not provide a sum rule for the combination ϕ (1) a entering into the matrix element (6.13).On the other hand, one cannot use the Schrödinger approximation for all modes, even as a rough estimate: We are going to see that this approach would produce the divergence of the integral (6.14) at κ → 0 and thus would lead to a qualitatively wrong answer.
Upon decomposing the states into multipoles and using the spherical symmetry of the background, we can integrate over angles which leaves us with the expression, 8 where j is the spherical Bessel function of the order .The sum here is taken over all excitations with strictly positive energy n > 0. It includes discrete modes and the states from the continuum.The monopole and dipole sectors contain also zero modes which are not included in the sum and require special treatment.We evaluate the contributions into C ls separately for = 0, = 1 and ≥ 2.
Monopole = 0. We calculate the = 0 contribution as follows: First, we replace the exact wavefunctions by the Schrödinger wavefunctions and apply the completeness relation of the latter to obtain C (Sch) ls ( = 0).Second, because for the low-lying modes the two wavefunctions differ significantly, we consider the corrections to C (6.16) Note that we do not include the n = 0 state since its energy is degenerate with the soliton.The same procedure applies to calculating the other multipole contributions.
In more detail, before replacing the wavefunctions by the Schrödinger ones, we use the orthogonality with respect to the inner product (3.13) between the ground state ϕ 00 (see eq. (3.8)) and the excited states which implies dr r 2 ϕ (1) 0n (r) + ϕ 0n (r) χ 0 (r) j 0 (κr) − 1 .(6.18) 8 In deriving this formula we have used the decomposition of the plane wave into spherical harmonics, where κ, ξ are unit vectors in the direction of κ and ξ.
This modification makes this matrix element manifestly of order O(κ 2 ) at small κ and will allow us to interchange the order of integration in C ls ( = 0) without encountering infrared divergences.
We are now ready to make the replacement ϕ and use the sum rules following from eq. (3.31), where on the r.h.s.we have explicitly subtracted the contribution of the ground state wavefunction ϕ (Sch) 00 and used that it is proportional to the soliton profile χ 0 (r).Inserting this completeness relation into eq.( 6.15) we obtain, It is convenient to interchange the order of integrations over κ and r, so that the integrals over κ involving only the Bessel functions can be performed analytically.Then the remaining integrals over the radial coordinate are taken numerically and give Having the numerical solutions of the exact and Schrödinger wavefunctions given in section 4 and shown in fig.2, we evaluate the corrections ∆C s ( = 0, n) for the first few modes.For a single mode, the correction takes the form 0n (r) + ϕ Note that by far the largest correction comes from the fundamental mode n = 1, whereas the corrections from higher modes quickly decrease.Adding these corrections to eq. ( 6.21) we will obtain the monopole contribution C ls ( = 0).Before presenting its numerical value it is natural to ask about its uncertainty coming from the fact that we are still approximating the high-level modes by the Schrödinger wavefunctions.To estimate the uncertainty, we compute the relative correction for the highest mode, which happens to be about 44%.We expect the relative correction to decrease for higher levels since the exact and Schrödinger wavefunctions approach each other for higher frequencies. 9herefore, the uncertainty is estimated as, where C (Sch) ls ( = 1) excludes n = 0 mode, like for the monopole, since the zero mode has degenerate energy with the soliton.However, a naive replacement of the exact wavefunctions by the Schrödinger wavefunctions leads to a spurious divergence.Indeed, in the absence of any orthogonality relations between the dipole Schrödinger modes and χ 0 (r), the matrix element dr r 2 ϕ (Sch) 1n (r)χ 0 (r)j 1 (κr) behaves as O(κ) at small κ for any n.The integral over κ in eq.(6.15) then logarithmically diverges at the lower end.
The divergence is spurious and is removed by modifying the matrix element before substituting the Schrödinger modes.We observe that the exact modes with n > 0 are orthogonal to the zero-mode conjugate vector φ10 (see eq. (3.24)), implying 1n (r) rχ 0 (r) = 0 , n > 0 .(6.28) Thus, we can subtract the leading term in the Bessel function, without changing the value of the matrix element.This leads us to modify the form of the dipole contribution as follows: 1n (r ) χ 0 (r ) j 1 (κr ) .(6.30) Note that we have made the replacement (6.29)only in one matrix element, but kept j 1 (κr) in the other.This is sufficient to avoid divergence at κ → 0, and does not introduce a new spurious divergence at κ → ∞, which would occur if we made the replacement (6.29) in both factors.
Starting from the expression (6.30) the rest of the analysis proceeds in the same way as in the monopole case.Substituting the Schrödinger wavefunctions and using their completeness relation we obtain C The corrections from the exact low-lying modes are ∆C ls ( = 1, n) = 0.119 , 0.019 , 0.006 , 0.002 for n = 1, 2, 3, 4 .(6.32) The relative correction for the highest mode with n = 4 is 20%.Putting everything together we arrive at C ls ( = 1) = 0.67 ± 0.07 .(6.33)This is about a factor of two smaller that the monopole contribution (6.26), with a comparable uncertainty.
Higher multipoles ≥ 2. For the other multipoles, since there are no small or largeκ divergences, we sum the contributions together with the completeness relation without modifications, and include the corrections from the exact wavefunction, The exact wavefunction corrections are: ∆C ls ( = 2, n) = 0.181 , 0.005 , −0.002 for n = 0, 1, 2 , ( ∆C ls ( = 3, n) = 0.003 , 0.001 for n = 0, 1 , ( ∆C ls ( = 4, n = 0) ∼ 10 −5 .(6.38c) We observe that the corrections quickly decrease with the increase of and n.Still, it happens that the relative correction remains sizeable -of order 10%.Using it to estimate the uncertainty, we obtain C ls ( ≥ 2) = 2.43 ± 0.14 .(6.39)This is larger than the joint contribution of the monopole and dipole.We conclude that excitation of modes with higher multipoles is as important as that of monopole and dipole modes.Though contributions of individual states beyond quadrupole are tiny, the cumulative effect is significant due to a large phase space consisting mostly of unbounded waves.
We now have all the multipole contributions in hand.Summing eqs.(6.26), (6.33), (6.39) we obtain our final result C ls = 4.3 ± 0.2 , (6.40) where we have added the uncertainties in quadratures.We believe this treatment of uncertainty is reasonable since we observed that the difference between the exact and Schrödinger contributions can have different signs for different modes.The light soliton evaporation rate found in this section is consistent with the result of 3D wave simulations of soliton immersed in axion gas [25].

Conclusions
We systematically studied the normal modes of boson stars.The negative and positive frequency perturbations are coupled together through gravity.It causes significant differences between the exact normal mode wavefunctions and the solutions of the Schrödinger equation which neglects the mixing between the negative and positive frequencies.One of the differences shows up in the zero modes.The exact spectrum contains monopole and dipole zero modes, corresponding respectively to the U (1) and translation symmetry breaking by the soliton background, whereas only a monopole zero mode is present in the Schrödinger spectrum.
We derived the general properties of the soliton perturbations.The coupling between positive and negative frequencies leads to non-standard orthogonality and completeness relations.We laid out a numerical approach to systematically find exact wavefunction solutions by solving the Schrödinger-Poisson equations.As far as we know, this is the first time that the complete first few wavefunctions are presented.For higher-level modes, the solutions are getting closer and closer to the Schrödinger wavefunctions.This observation turns out very useful for evaluating the scattering rates of soliton and axion gas.Further, we performed the 3D wave simulation of soliton and its perturbations.The power spectrum of the perturbations has clear resonant peaks at the frequencies of the few first normal modes, in excellent agreement with the numerical wavefunction solutions.
Our results are useful for understanding the evolution of axion solitons in various environments.As a concrete example, we studied the evaporation rate of a light soliton immersed in a gas of axion waves.We demonstrated that the properties of the exact normal mode wavefunctions are essential to obtain an infrared finite result.We also exploited the fact that the Schrödinger wavefunctions become a good approximation for high-level modes to evaluate an infinite sum over final states in the expression for the rate.The final results are consistent with the 3D wave simulations [25].It will be interesting to generalize this study to other physically relevant situations, such as the growth of the heavy soliton by kinetic relaxation or dynamics of fuzzy dark matter solitons in galaxies.The approach developed in this work can also be applied to boson stars in more complicated models, such as scalar theories with non-linear potentials [59][60][61], or massive vector fields [62][63][64].

Figure 2 .
Figure 2. The numerical solutions of the linearized Schrödinger-Poisson (SP) system and the Schrödinger equation for = 0, 1, 2, 3, 4. Solid magenta lines represent the positive frequency component of the SP eigenstate, orange dashed lines -its negative frequency component, and black dashdotted lines the perturbation of the gravitational potential.The turqoise lines show the Schrödinger wavefunctions.For each multipole, the first few energy eigenstates are displayed.For visualization purpose we normalize ϕ(1) 10 to have unit derivative at r = 0, which differs by a constant factor from eq. (3.8).

=
−χ, where the coefficient β00 is due to different normalization of the Schrödinger and exact modes. β

Figure 3 .
Figure3.The power spectrum of = 0 perturbation (black lines) at four locations: r = 0, x = 5, y = 5, and z = 5.The vertical dotted red lines correspond to the predicted binding energy of = 0 modes and n = 1 to 4; the vertical dotted blue lines correspond to the binding energy of = 1 modes and n = 1 to 3; the vertically dotted green lines correspond to the predicted energy levels of = 2 modes and n = 0 to 2. The grey-shaded region represents the start of the continuum spectrum.The same color coding applies to the power spectrum analysis of = 1 and = 2 shown in figs.4 and 5.

Figure 4 .
Figure 4.The power spectrum of = 1 perturbation (black lines) at four locations: r = 0, x = 5, y = 5 and z = 5.The same color coding as fig. 3 is used.Notice the vast difference in vertical scale in different panels.

Figure 5 .
Figure 5.The power spectrum of = 2 perturbation (black lines) at four locations: r = 0, x = 5, y = 5 and z = 5.The same color coding as fig. 3 is used.Note the difference in vertical scale between the first and the remaining panels.

. 25 )Dipole = 1 .
In this way we obtain the monopole contribution with the uncertainty, C ls ( = 0) = 1.21 ± 0.08 .(6.26) Similarly, the dipole contribution to C ls has the Schrödinger wavefunction summation and the exact wavefunction corrections,