New Solutions for Rotating Boson Stars

It has been shown that scalar fields can form gravitationally bound compact objects called boson stars. In this study, we analyze boson star configurations where the scalar fields contain a small amount of angular momentum and find two new classes of solutions. In the first case all particles are in the same slowly rotating state and in the second case the majority of particles are in the non-rotating ground state and a small number of particles are in an excited rotating state. In both cases, we solve the underlying Gross-Pitaevskii-Poisson equations that describe the profile of these compact objects both numerically as well as analytically through series expansions.

Rotating boson star configurations have also been studied, but all known solutions (that we have found in the literature) have the property that the total angular momentum increases proportionally to the mass of the star (e.g. [23][24][25][26][27]). In these solutions, the ratio of the angular momentum to the number of particles has a minimum value, and hence for a fixed number of particles, these solutions do not include configurations of rotating boson stars with an arbitrarily small angular momentum.
In this paper we remedy this gap, by finding new solutions which carry an arbitrarily small angular momentum for a fixed number of particles. We in fact find two different classes of such solutions.
Our first approach is a generalization of the solutions which exist in the literature, where all the particles are in the same state. However, we impose that the total angular momentum in the bosons is constrained to be fixed at a small value. This produces a state dominated by a spherical component, with a small admixture of a higher harmonic, naturally leading to a star with a small rotation. Our second approach is to take a small number of particles in the star to be in a higher spherical harmonic, while most of the particles are in the non-rotating state. Note that it is clear that such a solution must exist; for instance if a single particle is placed in a = 1 harmonic, there is no lower energy state with this angular momentum. We shall call these two ansätze respectively the one-state and two-state solution. We show that both these approaches successfully yield solutions for a rotating star with a small angular momentum.
his paper is organised as follows: In order to set our notation and to connect to previous work, we first review our previous results for the case of non-rotating boson stars in Sec. II. We then turn to the rotating star: we consider the one-state ansatz in Sec. III and the two-state ansatz in Sec. IV. In each case, we set up the perturbation expansion around the non-rotating star, and solve the equations both numerically and in a series expansion, thereby providing strong numerical evidence that these solutions exist. We conclude in Sec. V.

A. Lagrangian and Structure Equations
Let us consider a real non-interacting scalar field φ(r, t) which is coupled to gravity. This scenario is described by the following Lagrangian The scalar field can form gravitational bound states, or boson stars. In this study, we focus on the case of dilute boson stars, which can be described by the Newtonian non-relativistic limit. For the case of QCD-axions, it has been shown that only dilute axion stars are stable over astronomical time scales [28,29]. In the Newtonian limit, when the field φ couples only weakly to gravity, the metric can be written as g µν = diag(1+2Φ, −1, −1, −1), where Φ is the Newtonian gravitational potential. We are interested in stationary solutions, in which case the gravitational potential is time independent. In this case the Ricci scalar takes the simple form R = −2(∇Φ) 2 . Also in the non-relativistic limit, we can treat the energy as being close to the mass, and we have (∂ t φ) 2 Φ = m 2 φ 2 Φ. The Lagrangian in Eq. (1) then becomes Since the Lagrangian is quadratic in the scalar field, we can quantize the scalar in the usual way. We first find a set of wavefunctions satisfying and quantize by setting the scalar operator equal to φ = a † n φ n + a n φ * n .
The Hamiltonian is then The eigenstates are of the form The gravitational potential interacts with the scalar through the term m 2 φ 2 Φ. This leads to the equation for the potential These field equations, often referred to as Gross-Pitaevskii-Poisson equations, are the structure equations for the boson star.

B. The non-rotating Boson Star
For the non-rotating star, we consider an ansatz where we have N particles in the ground state ψ nr ≡ φ 0 , which has an energy eigenvalue e nr ≡ E 0 . The state is then and the corresponding structure equations are given by a Schrödinger type equation for the ground state wavefunction and a Poisson equation for the gravitational potential To solve the structure equations for the boson stars, it is convenient to introduce dimensionless variables. Following Refs. [17,18], we define where M = N m is the star's mass. Using these variables, we can rewrite the Gross-Pitaevskii-Poisson equations as where the derivatives ∇ z are now with respect to the dimensionless coordinate z. In dimensionless variables, we can write the normalization condition of the wavefunction, |ψ| 2 dV = 1, as s 2 0 z 2 dz = 2β and associate 2β with the mass of the star. Note that up to scalings, there is only one ground state solution for non-interacting boson stars.
In [17,18] we have solved the Gross-Pitaevskii-Poisson equations and obtained a semi-analytic solution for the ground state of the boson star. In this approach, the profiles at both small and large radii are separately described through a series expansion of the wavefunction and potential and matched at an intermediate point. At small radii, the profile can be described by an even polynomial around the center of the star (z = 0) At large radii, we take The potential and wavefunction are fully specified by knowing the parameters of the leading expansion The remaining coefficients can be obtained using recursion relations which can be derived from the Gross-Pitaevskii-Poisson equations and have been presented in [17,18]. We can determine the four expansion parameters either through a fit to the numerical solution, or by matching the small and large radius wavefunction and their derivatives at a matching point z * . We have obtained the following solutions [17] s 0 0 = 1.02149303631 ± 1.4 · 10 −10 v 0 0 = 0.93832284019 ± 1.3 · 10 −10 α = 3.4951309897 ± 5.1 · 10 −9 β = 1.7526648513 ± 1.3 · 10 −9 . (16)

C. Slowly Rotating Boson Stars
We now turn to a study of rotating boson stars. In particular, we look for slowly rotating boson stars solution, which can be treated as a perturbation around the non-rotating solution. That is, the non-rotating solution should admit a normalizable perturbation such that the perturbation carries angular momentum. The existence of such a perturbation would indicate that a slowly rotating boson star can be found at least at the linearized level, which is suggestive that the full solution should exist.
To look for these states, we impose a constraint on the total angular momentum of the boson star whereL 2 is the usual total angular momentum operator On astrophysically relevant times scales, the boson star's angular momentum L star is a fixed quantity. We implement this constraint by introducing a Lagrange multiplier µ. The Lagrangian in Eq. (2) then becomes We can repeat the quantization procedure and find resulting equations of motion are the Poisson equation given in Eq. (10) and a modified Schrödinger-type equation In the following, we will present two possible solutions for the slowly rotation boson star, and obtain the corresponding ground state wave-function.

A. The Ansatz
We first look for a solution where all the particles are in the ground state. The state is then This is formally similar to the non-rotating case, but because of the constraints, we must take the ground state in this sector to have non-zero angular momentum. We take the ground state to be a perturbation around the non-rotating spherically symmetric solution ψ nr (r) obtained in Sec. II B. In particular, we choose an ansatz in which the wavefunction and potential perturbation are expanded in spherical harmonics Y m with ≥ 1 and m = 0, as well as −E 0 = e nr + e 1 . The expansion parameter is taken to be parametrically small, which allows us to work in linear order perturbation theory.
The angular momentum constraint in Eq. (17) relates the value of and the star's angular momentum, such that = L star × [N ( + 1) · |ψ 1 | 2 dV ] −1/2 . We then find that the last term in Eq. (19) is of order L 2 star ∼ 2 and can therefore be ignored at linear order in perturbation theory.
We now insert this ansatz into the field equations Eq. (19) and Eq. (10). Collecting terms at zeroth order in , we recover the equations of motion for a non-rotating boson star, whose solution we presented in Sec. II B.
Matching the terms proportional to Y 0 we find the structure equations for the perturbation Finally, collecting the terms proportional to Y 00 implies e 1 = 0, meaning that the rotation does not induce a shift in the binding energy at leading order in perturbation theory. Such a shift first appears at order 2 .
We perform the change of variables in Eq. (11) and further define The resulting structure equations for the dimensionless field and potential perturbations s 1 and v 1 then read

B. Series Expansion
We have seen in Sec. II B that we can describe the profile of the non-rotating boson star through an infinite series for the wavefunction and potential. We will follow the same approach to obtain a solution for Eq. (24).
At small radii, the profiles for s 1 and v 1 can be described via a polynomial around the center of the boson star z = 0, By matching the coefficients in Eq. (24) we obtain the recursion relations Requiring the left hand side of Eq. (24) to be defined at z = 0 implies that the perturbation vanishes at the origin and hence s 1 0 = v 1 0 = 0. The profile at small radii can therefore be fully parameterized in terms of the derivative of the wavefunction and potential at the origin ∂ z s 1 = s 1 1 and ∂ z v 1 = v 1 1 . At large radii, we will once again use the series expansion ansatz and obtain the recursion relations and Setting m = 0, Eq. (30) can be written as s 1 1,0 Γ = 0, which either implies s 1 1,0 = 0 or Γ = 0. Although both possibilities will lead to a solution, we will mainly focus on the Γ = 0 solution. For m = 1, we find that 1 − σ = β, where the σ originates from the s far 1 expansion. This is the same relation we found for the non-rotating boson star in Ref. [17], justifying our ansatz to use the same σ for both the s far 0 expansion in Eq. (14) and s far 1 expansion in Eq. (27). Finally, setting m = M +1 we can use Eq. (30) to obtain the recursion relation This means, that all coefficients can be determined recursively from s 1 1,0 and v 1 0, . More generally, we can use Eq. (28) and Eq. (29) to recursively calculate all coefficients s 1 n,m and v 1 n,m in the expansion of s far 1 and v far 1 . We can now determine the expansion parameters by matching the near and far field wavefunction and potential and their derivatives at a matching point z * . We have performed such a matching using the near field solution in Eq.
where the perturbation is normalized such that v 1 1 = 1. To estimate the uncertainty associated with the matching procedure, we performed multiple matchings for 3 ≤ z * ≤ 3.5.

D. Numerical analysis
In Sec. III B we have shown that the wavefunction and potential profile of the rotating boson star can be described by a series expansion, which is characterized by the expansion parameters given in Eq. (32). In the following, we will compare this result to the numerical solution of Eq. (24), focusing on the case = 1.
As we have seen before, near z = 0 the solution takes the form s 1 ∼ s 1 1 z + . . . and v 1 ∼ v 1 1 z + . . . . To obtain a numerical solution, it is convenient to normalize the field s 1 and the potential v 1 such that v 1 1 = 1, so that the solution is only parameterized by s 1 1 . Using a Runge-Kutta 4 method, we then perform the numerical integration of Eq. (24). For most values of s 1 1 , the wavefunction profile will diverge to positive or negative infinity at large radii z 1. Using a shooting point method analogous to those used by the authors of Ref. [31] and [17], we adjust s 1 1 such that the wavefunction converges and becomes square integrable.
The numerical solution for Γ = 0 is shown in the left panel of Fig. 1 We can see that already such few terms in the series expansion are sufficient to describe the wavefunction well.
The dotted black curve shows the Whittaker function solution of Eq. (36), which is already well described by the first few terms of the far field expansion. The right panel of Fig. 1 shows the numerical solution for both Γ = 0 and Γ = 0, alongside with the non-rotating ground-state solution s 0 discussed in Sec. II B. In particular, we found that solutions exist for Γ = −0.667, −0.812, −0.887. These values are consistent with the relation Γ = β 2 /κ 2 − 1 found in Sec. III C for κ = 3, 4, 5. Notably, κ−2 also characterizes at how many radii the wavefunction vanishes identically, s 1 = 0.

IV. ROTATING AXION STARS: TWO-STATE SOLUTIONS
A. The Ansatz We now consider a second approach to find rotating boson star solutions. In this ansatz, we look for a state where N particles are in the ground state a † 0 |0 , and k particles are in the excited state a † 1 |0 . The state is then This leads to the Poisson-type equation for the potential which should be solved along with the Schrödinger-type equations We will assume k N , and perturb in the small parameter = k/N . For this reason, we dropped the L 2 star term in Eq. (41) which only contributes at subleading order in . Now, to zeroth order in , Φ will just be equal to the potential for the non-rotating star Φ nr , φ 0 is equal to the wavefunction for the non-rotating star ψ nr , and E 0 = e nr . As before, we will consider a single Y 0 mode i.e.
We again perform the change of variables in Eq. (11) and further define and obtain the structure equation The angular momentum of the boson star is equal to L star = [k· ( +1)· |ψ 1 | 2 dV ] 1/2 . Note that unlike for the one-state case, in this case the rotation does induce a shift in the binding energy at leading order in perturbation theory.

B. Series Expansion
As before, we will parameterize the wavefunction via an infinite series expansion. At small radii, the profile for s 1 can be described via a polynomial around the center of the boson star z = 0, By matching the coefficients in Eq. (44) we obtain the recursion relation As in the one-state case, requiring the left hand side of Eq. (44) to be defined at z = 0 implies that the perturbation vanishes at the origin and hence s 1 0 = 0. The profile at small radii can therefore be fully parameterized in terms of the derivative of the wavefunction at the origin ∂ z s 1 = s 1 1 . At large radii, we will use the series expansion Note that the form of this ansatz is slightly different than for the non-rotating boson star in Eq. (14) and the onestate solution in Eq. (27). As we will see later, two-state solutions only exist for Γ = 0, and the resulting far field solution would approximately follow the Whittaker function W κ,µ (2z √ 1 + Γ). In order to match the asymptotic behaviour of this Whittaker function solution, the additional √ 1 + Γ factor as well as a new parameter σ have been included in the series expansion ansatz.
The coefficients of the expansion are related by the recursion relation Here we have used the approximate ground-state potential v f ar 0 ≈ −1 + 2β/z, such that the Cauchy product v 0 s 1 is well defined.
As before, requiring the wavefunction to be normalizable implies that all coefficients s 1 0,m vanish. The first non-vanishing terms appear for n = 1, in which case we can simplify Eq. non-rotating boson star expansion in Eq. (14). Setting m = M + 1 we find Following the same procedure as in Sec. III B, we determine the expansion coefficient s The uncertainty was estimated by performing multiple matchings for 3.2 ≤ z * ≤ 3.5.

C. Leading Order Analytic Far-Field Solution
Similar to the ground-state and the one-state solution, we can obtain an approximate analytical solution for the far field at leading order n = 1. Using v 0 ≈ −1 + 2βz −1 and performing a change of variables to w = 2zs 1 and y = 2z(1 + Γ) 1/2 , we can rewrite Eq. (44) in the familiar Whittaker equation form where κ = β · (1 + Γ) − where µ 2 = 1 4 + ( +1). The function M κ,µ (y) diverges at large y unless Γ = β 2 /κ 2 −1 with κ being a natural number ≥ 2. We will see in the next section that the solutions and corresponding values of Γ do not fulfill this condition. Normalizability of the wavefunction then requires c = 0. The solution must therefore be solely described by W κ,µ (y) which allows us to write Expanding this function leads to the ansatz in Eq. (47) and matching the coefficients of the leading terms allows us to identify c = s 1 1,0 [2 −σ √ 1 + Γ] −1 with σ = 1 − κ.

D. Numerical results
As for the one state solution, we also obtain a numerical solution of Eq. (44), focusing on the case = 1. The equation is linear in s 1 which allows is to choose s 1 1 = 1 without loss of generality. We then use a Runge-Kutta 4 method to perform the numerical integration of Eq. (44) and apply a shooting point method to find the values of Γ for which the wavefunction converges at large radii.

V. CONCLUSIONS
Light scalar fields can form gravitationally bound compact objects, called boson stars. In the Newtonian limit, the profiles of boson stars are described by the Gross-Pitaevskii-Poisson equations.
In previous works, we presented a semi-analytic solution to these equations describing the profile of boson stars formed by scalar fields [17,18]. The solution was based on a series expansion which is parametrized by four expansion parameters that were obtained from numerical simulation at high accuracy. In this paper we have extended our methods to find new solutions which allow for slowly rotating boson stars; specifically, we have found solutions for boson stars where the ratio of the angular momentum to the number of particles can be made arbitrarily small. We considered two possibilities; in one case, all the particles are in the same state and in the second case the majority of the particles are in the zero angular momentum ground state and a small number of particles are in an excited state containing angular momentum. In each case, we obtained accurate numerical and semi-analytic profiles (about 1% precision), thereby establishing the existence of these slowly rotating boson stars.
The results and methods presented in this paper allow for systematic studies of the properties of boson stars in an analytic way without further relying on numerical simulations. There are several directions for further research; in particular, it would be interesting to extend these solutions to interacting scalars and to relativistic stars. It would also be interesting to see how the profiles are modified in the presence of other astrophysical objects like planets. We hope to return to these questions in future work.