On Profiles of Boson Stars with Self-Interactions

Under the influence of gravity, light scalar fields can form bound compact objects called boson stars. We use the technique of matching asymptotic expansions to obtain the profile for boson stars where the constituent particles have self-interactions. We obtain parametric representations of these profiles as a function of the self-interactions, including the case of very strong self-interactions. We show that our methods agree with solutions obtained by purely numerical methods. Significant distortions are found as compared to the noninteracting case.

It has been noted by many authors that the axions might bind into compact spatial structures (e.g. see the review [16]). These are more generally referred to as boson stars [17][18][19]). Such objects, if they exist, would produce distinctive signatures in axion search experiments, and understanding these signatures requires a description of the profile of the boson stars, which could produce unique time and spatial dependencies of the signal * fkling@uci.edu † arajaram@uci.edu which can distinguish it from backgrounds [20]. It is also important to know how these profiles are distorted by the presence of self-interactions, and also by the gravitational effects of matter. For all these reasons, it is timely to have a precise description of the profiles of boson stars.
In a previous communication [21], we analyzed the profile of such objects in the special case when the bosons had no self-interactions. These objects are solutions to a coupled set of equations called the Gross-Pitaevskii-Poisson equations. In that paper, we showed that one could find approximate solutions to these equations, by using a combination of analytical and numerical methods. We showed that our methods were numerically stable, and that they converged uniformly far away from the core of the star, and furthermore were much less computationally expensive compared to other purely numerical methods [18,[22][23][24][25][26][27][28][29][30][31][32][33][34][35].
In this paper, we extend our previous methods to the case of interacting bosons. This is a more difficult situation, because the self-interactions can be much stronger than the gravitational binding. Nevertheless we show that our methods continue to provide excellent agreement with a fully numerical solution.
Our results show that self-interactions can produce significant distortions of the profile of the star, as compared to the non-interacting star. We parametrize these deformations by finding a parametrization for the star in terms of two matching asymptotic expansions, where the expansions are taken far away from the core and close to the core respectively. The parameters of the solution can be found systematically by matching the expansions in an overlap region. We find these parameters as a function of the interaction strength, which also allows us to find any desired physical quantity (mass, central density etc.) in terms of the coupling. We do this both for weak couplings, where we can perturb around the noninteracting star, as well as for strong couplings, where we can perturb around the Thomas-Fermi limit of the solution. We note that these expansions provide a solution to the axion stars which have an accuracy of at least 10 −3 , and can hence be used in lieu of complicated numerical calculations.
Our results improve on other semi-analytic approaches (e.g. [33,36,37]), which have used variational and other techniques to find approximations to the boson star profile. Our methods are particularly suited to accurate evaluations of the profile away from the center, where the falloff is well described by a Whittaker function.
In the following sections we rederive the Gross-Pitaevskii-Poisson equations satisfied by the boson star. We also describe various limiting cases, including the Thomas-Fermi limit of strong coupling. We then present a series expansion in the two asymptotic regimes, and find the solutions separately for weakly coupled systems and for strongly coupled systems. We show how our solutions apply to the special case of axion stars. We end with conclusions and directions for future work.

II. GROSS-PITAEVSKII-POISSON EQUATIONS
A. The real scalar field In [21] we have derived the structure equations for the ground state of a self-gravitating complex scalar field in the non-relativistic limit. Following the procedure described in [38], we now show that the same set of equations also apply for a real scalar field.
Let us consider the real scalar field φ( r, t) described by the Lagrangian In the presence of gravity, the scalar field can form gravitational bound states, which are called boson stars. A simple solution can be obtained assuming that only the ground state is populated. In this case the field can be expressed in terms of a single real function ψ(r), sometimes called the wave function of the boson star, describing the radial profile of the boson star. We can write where E is the ground state energy and N is number of bosons in the ground state. Note that properly shifting the time coordinate allows us to absorb a possible phase of the wave function and therefore to choose ψ to be real. We have chosen a wave function normalization ψ 2 dV = 1 which allows us to identify ψ 2 with as the probability density.
The equation of motion of the scalar field is the Klein- Assuming that the field couples only weakly to gravity, we can use a Newtonian approximation. This allows us to introduce the Newtonian potential Φ in the metric g µν = diag(1 + 2Φ, −1, −1, −1). We can then rewrite the Klein-Gordon equation as Let us further assume that the ground state is nonrelativistic. In this case we can write E = m + e with binding energy e m. This implies eψ, Φψ, ∇ψ mψ.
Using ∂ 2 t φ = −E 2 φ and E = m + e, we can rewrite the the Klein-Gordon equation in the non-relativistic limit as Inserting the explicit form of the field given in Eq. 2 and rephasing by e iEt , we obtain the Schrödinger-type equation Here we have dropped additional terms with rapidly oscillating phase factor e −inEt where n is a non-zero integer. For the non-relativistic approximation to be consistent, the last term should be sufficiently small, i.e.
N λ 4m 2 m. The Newtonian potential is related to the energy density via the Poisson equation ∇ 2 Φ = 4πGρ. The energy density ρ of the real scalar field is where we used the non-relativistic approximation in the last step. Newton's equation therefore takes the simple form Comparing with the results of [21], we find that both the ground state of a boson star for both a complex scalar field and a real scalar field are described by the same set of equations given in Eq. 5 and Eq. 7. These are often referred to as Gross-Pitaevskii-Poisson equations. Note that the Gross-Pitaevskii-Poisson equations for a real scalar field have also been obtained by the authors of [29], which follow a semi-classical approach considering a quantized scalar field φ.

B. Limits and Validity
The terms on the right-hand side of Eq. 5 represent the contribution to the energy of a scalar particle due to the quantum pressure, gravity and classical pressure respectively. The quantum pressure is a consequence of Heisenberg's uncertainty principle and is always repulsive, preventing the star from gravitational collapse. Gravity on the other hand is always attractive. The classical pressure arises from he contact interaction term and can either be attractive or repulsive, depending on the sign of the self coupling λ. It is illustrative to consider the limits in which one of the three contributions is negligible.
Using scaling relations between the coupling λ, the star's mass M and the star's radius R, we will qualitatively discuss both physical properties of the star and the validity of the solution. Without providing a formal definition of the star's radius R, we know that the probability density inside the star scales like ψ 2 ∼ 1 R 3 . Similarly, a field derivative will scale as ∇ψ = ψ R . For a more rigorous discussion of the mass-radius relations of boson stars see [36].

Non-Interacting Limit λ = 0
In the non-interacting case, λ = 0, the quantum pressure balances gravity. We have obtained a semianalytical solution for this case in [21]. Note that for non-negligible couplings λ = 0, the boson star becomes effectively non-interacting at large radius due to low densities.
We can rewrite Eq. 5 and 7 as Using the scaling behavior discussed above, we see that the radius of a non-self-interacting boson star scales as R ∼ (GM m 2 ) −1 . This is an remarkable result: the star's radius decreases when its mass increases. The binding energy of a scalar particle is The non-relativistic approximation requires e m which

Thomas-Fermi Limit
For strong repulsive self-interactions λ > λ c , the quantum pressure becomes negligible and the classical pressure balances gravity. Here λ c is the critical coupling at which the quantum pressure and classical pressure are equally important: We can use the scaling behavior introduced earlier to solve for the coupling and obtain λ c ∼ Rm 2 M . At the critical coupling, hydrostatic equilibrium requires R ∼ Already a very small coupling is sufficient for the Thomas-Fermi limit to apply.
For λ λ c we can rewrite Eq. 5 and 7 as Using that the gradient scales as ∇ψ ∼ ψ R , we obtain that the radius of a strongly self-interacting boson star scales as R ∼ λ. The radius is independent of the star's mass. This is not surprising, since both gravity and the repulsive self-interaction are proportional to the number of particles and therefore the star's mass. The binding energy of a scalar particle is The non-relativistic approximation requires e m . Larger couplings increase the validity range of the nonrelativistic approximation to higher masses M .
An analytic solution for the Thomas-Fermi limit has been obtained in [23] and is discussed in Sec. IV.

Non-Gravitational Limit
For sufficiently strong attractive self-interactions, λ < λ * , quantum pressure balances the attractive selfinteraction while the effect of gravity becomes negligible. Gravity and the classical pressure become equally important at λ * when Using Φ ∼ GM R and the radial size R ∼ (GM m 2 ) −1 this implies a critical coupling of For a given coupling λ, the non-gravitational limit applies for stars with R < R * where the critical radius with a corresponding critical mass M * = (Gm 2 R * ) −1 = M pl |λ| − 1 2 . For λ < λ * we can rewrite Eq. 5 as Using the scaling behavior, we see that the radial size of the star in the non-gravitational limit is R ∼ M |λ| m 2 . The radius increases linearly with the mass of the star. Since at larger radius R > R * we approach the non-interacting limit in which the mass decreases for increasing radius, we find that there is a maximum possible mass for boson stars M max = M * ∼ M pl |λ| − 1 2 . However, it has be shown in [36], that the solutions for R < R * are unstable with respect to perturbations. Therefore boson stars in the non-gravitational limit cannot be realized in nature, at least for this simple class of interactions. For more complicated interactions, such stars can exist, and fall under the general class of Q-balls [39][40][41][42].

C. Scaling Invariance
Following the terminology of [21] we introduce the dimensionless variables We can then rewrite the Gross-Pitaevskii-Poisson equations given in Eq. 5 and 7 as The wave function normalization condition ψ 2 dV = 1 becomes where M = N m is the mass of the boson star. Let us note that Eq. 16 and 17 are invariant under the scaling where f is a scaling factor. This implies that we can relate different solutions of the Gross-Pitaevskii-Poisson equations corresponding to different boson star masses M and couplings Λ through rescaling. We will make use of this scale invariance and solve for the Gross-Pitaevskii-Poisson equations at a fixed reference scale k.
A particularly useful choice for our discussion is to set −k 2 = V (∞) = e 2m which transforms as k → f k. We can then introduce the scale invariant coordinate z, wave function s, potential v, mass β and coupling γ via Using the scale independent variables, we can write the Gross-Pitaevskii-Poisson equations as The scale choice implies the boundary condition v(∞) = −1. The solution corresponding to a boson star with mass M can be obtained by performing the rescaling given in Eq. 19 with k = GM m 2β . In the following section, we will obtain an approximate analytical form for s, v.

A. Series Expansion
We have seen in the discussion of the non-selfinteracting case, that we can describe the profile of the boson star through an infinite series for the wave function and potential.
Following the same approach as for the non-selfinteracting case [21], we will describe the profile at both small and large radii through a series expansion of the wave function and potential. At small radii, the profile can be described via an (even) polynomial around the center if the boson star z = 0 Eq. 20 then leads to the recursion relations The smoothness of the profile at the origin implies s 1 = v 1 = 0, and therefore, also all odd coefficients s 2n+1 , v 2n+1 vanish. The profile at small radius z can therefore be fully parametrized in terms of the wave function and potential at the origin: s 0 and v 0 . At large radius, we can expand the profile using the series expansion By matching the coefficients, Eq. 20 we obtain the recursion relations n,m p,q=0,0 Let us note the following properties of the solution: i) Normalizability requires s 0 0 = 0. Eq. 24 then implies that all coefficients s 0 m vanish as well. This means that the wave function decays at least exponentially. ii) Eq. 25 then implies that all v 0 m = 0 for m > 1. This means that at large radius, the potential is described by the Here we both used the boundary condition imposed by our scale choice, v 0 0 = −1, and used the notation introduced in Eq. 19, v 0 1 = 2β. All other terms in the expansion of v are at least exponentially suppressed. iii) Setting n = m = 1, Eq. 24 can be written as 2β = v 0 1 = 2(1 − σ). This is a remarkable result: the exponent σ in the series expansion is related to the the total mass of the system σ = 1−β. iv) As derived in [21], the leading order, n = 1, solution for the wave function at large radius is given by the Whittaker function Here we have introduced the normalization parameter α = s 1 0 . v) The large radius solution can be fully parametrized by the expansion parameters α = s 1 0 and β = v 0 1 . The remaining coefficients can then be computed using Eq. 24 and 25. Note however, that the series expansion in Eq. 23 does only converge for m < M , where M is finite. vi) Eq. 24 and 25 further imply that the potential contains only non-vanishing components v n m for even n while the wave function only has non-vanishing components s n m for odd n. For practical purposes, we will truncate the infinite series in Eq. 21 and 23 and only take into account the leading terms with n ≤ N and m ≤ M . Let us define the truncated series expansion at small and large radius z as

B. Expansion Parameters
As we have seen in the previous section, the series expansion for s and v at small and large radius can be fully parametrized by the four expansion parameters α, β, s 0 , v 0 . In the non-interacting case γ = 0 these were just numbers, while in the general case they will be functions of the coupling γ α(γ), β(γ), s 0 (γ), v 0 (γ).
Following the strategy from [21], we first obtain the expansion parameters from numerical simulations. In a second step, we provide an analytic form for the expansion parameter.
To obtain a numerical solution of the boson star's profile, it is convenient to solve the Gross-Pitaevskii-Poisson equations as given in Eq. 16 using the boundary condition V (0) = 1. The authors of [43] have shown that the solutions of Eq. 16 can then be parametrized by the central value of the wave function, S 0 = S(0) and categorized into three distinct classes: for S 0 > S * 0 the wave function diverges for at large radius towards positive infinity, for S 0 = S * 0 the wave function converges to zero, is positive definite and square integrable, while for S 0 < S * 0 the wave function diverges at large radius towards negative infinity.
Using a Runge-Kutta 4 method with constant step size ∆ x , we perform the numerical integration until the wave function starts to diverge and iteratively optimize the central value of the wave function S 0 to find S * 0 . The precision of the wave function needed for the numerical solution to stay finite until a large value of x, which is needed to fit the large range solution, increases exponentially with the radial coordinate x. The accuracy of the numerical solution is limited by the step size ∆ x . In this study we use a precision of up to 150 significant figures and ∆ x = 10 −3 , providing an accuracy of the solution of order O ∆ 4 x ≈ 10 −12 . The Whittaker solution parametrization given in Eq. 26 always describes the wave function profile at large radius when the density is small enough that the selfinteraction becomes negligible. However, the central mass parameter β and the normalization α depend on the central profile of the star and therefore the coupling parameter γ.
To obtain the expansion parameters α, β as well as the scaling parameter k in Eq. 19, we fit the leading order profiles, the Whittaker solution S(x) = kα 2 β x W β,− 1 2 (2kx) and the Newtonian potential V (x) = −k 2 + 2kβ x , to the numerical solution for V and S at large x. To avoid systematic effects due to the truncation of subleading terms n > 1 of the series expansion in Eq. 23, we restrict the fitting range to x > x * , where the fraction of mass outside radius x * contributed less than 10 −12 to the total mass of the boson star. The expansion parameters at small radius and the coupling are obtained via s 0 = k −2 S 0 , v 0 = k −2 and γ = k 2 Λ.
In Fig. 1 we show the dependence of the expansion parameters α (upper left), β (upper right), s 0 (lower left) and v 0 (lower right) as function the coupling parameter γ as red dots. Note that the horizontal axis switches from a linear scale to a logarithmic scale at γ = 1 as indicated by the dashed black line. This indicated the transition between the weakly coupled regime γ < 1 and the strongly coupled regime γ > 1.
C. Fit for Weak Couplings: −1 < γ < 1 A vanishing self-coupling γ = 0 indicates the noninteracting limit. The corresponding results for the expansion parameters are presented in [21]. If the selfcoupling is weak, |γ| < 1, we can treat the classical pressure as perturbation. In this case we can write the the expansion parameters as a series expansion in the coupling γ around the non-interacting solution. To obtain the coefficients of this expansion, we fit a 6th-degree polynomial to the numerical solutions. We can write the result in an analytic form as The result is shown in Fig. 1 as a green line. The lower panels indicate the accuracy of the analytic form with respect to the numerical solution. We can see that for all four expansion parameters the solution from Eq. 29 reproduces the numerical results with accuracy better of O(10 −4 ) in the range −1 < γ < 1.
As mentioned before, boson stars with attractive selfcoupling and small radius R, or equivalently large negative coupling γ, become unstable with respect to perturbations. We show in Sec. V A, that this happens at γ min = −0.722. The solutions for γ < γ min are unphysical.

A. Thomas Fermi Limit
In the previous section we have discussed the weakly coupled scalar field. Let us now consider the case of a large repulsive self-coupling γ 1. In this case the quantum pressure becomes negligible and the classical pressure balances gravity. This scenario is known as Thomas-Fermi limit and has been examined in [23,44] We can write the Gross-Pitaevskii-Poisson equations as The (normalizable) solution for the profile is given by where sinc(x) = sin(x)/x. The wave function becomes zero at Z = π √ γ, implying that the boson star is compact and has a radius Z. At z > Z the wave function remains zero and the potential is described by the Newtonian potential v(z) = −1 + 2β z . Using that exterior and interior solution for the potential of the star have to match at surface, v(Z) = 0, we can solve for the mass parameter β = 1 2 Z = π 2 √ γ. We can further use the normalization condition (see Eq. 17) to obtain the central density s 0 = γ − 1 2 . Using Eq. 30, this implies that the central value for the potential is v 0 = 1.
So far, we have ignored the effects of quantum pressure to the boson star's profile. However, close to the star's radius Z, the density drops until eventually the self-interaction becomes negligible. The outer part of the star can therefore be described by the non-interacting solution, which at leading order can be written in terms of the Whittaker function. In the following, we estimate the remaining expansion parameter α by matching the Thomas-Fermi solution for the interior of the star with the Whittaker solution for the exterior part.
We define the matching point z * as the radius where the quantum and classical pressure terms become equal, ∇ 2 s(z * ) = γs(z * ) 3 . Using the Thomas-Fermi profile from Eq. 31, we can write this condition as We can write the matching point as z * = (π − δ) √ γ and expand the matching condition in δ. Keeping only the linear terms and using s 0 = γ − 1 2 , we can solve for δ and obtain δ = (4γ/π) − 1 3 . The wave function at z * has a value s(z * ) = s 0 sinc(π − δ) Following [45] (see Sec. 8, Eq. 18) we can write the Whittaker function W β,− 1 2 (2z) for z = 2β = π √ γ as Using the form of the exterior profile as given in Eq. 26 and matching it to the Thomas-Fermi solution at z * in Eq. 34 allows us to extract the expansion parameter α. We find We have already seen that in the Thomas-Fermi limit the remaining expansion parameters are given by We have therefore obtained a simple analytic form for all expansion parameters in the limit of large self-couplings. The results from Eq. 36 and 37 are shown in Fig. 1 as dashed cyan lines. Note that for increasing self-coupling γ, the expansion parameter α exponentially decreases leading to a sharp drop of the wave function profile at z = Z. The bosons become more and more confined in the inner part while the large radius tails of the wave function vanish.

B. Fit for Strong Couplings: γ > 1
In the above discussion of the Thomas-Fermi limit, we have neglected the effect of the quantum pressure on the structure of the star. For large but finite self-couplings γ, we can consider the quantum pressure as a perturbation to the Thomas-Fermi limit. We include this perturbation as a correction factor to the Thomas-Fermi solution, which we can express asseries expansion in the (inverse) coupling parameter. The coefficients of this expansion are obtained from a fit to the numerical solution. We can write the result in an analytic form as with the leading order Thomas-Fermi solutions Note that the correction factor for the expansion parameter α contains a non-unity constant term. This should not be surprising, considering that we obtained α in the Thomas-Fermi limit by matching the Thomas-Fermi solution and the Whittaker-solution at a matching point z * , even though both solutions describe the profile only poorly at this point. The results are shown in Fig. 1 as a blue line. We can see that for all four expansion parameters the solution from Eq. 38 reproduces the numerical results with accuracy of O(10 −5 ) for self-couplings γ > 1.
Using the analytic expression for the expansion parameters, we can now obtain the profile of the boson star for a given value of the self-coupling γ. For very large selfcouplings, γ 10, we can directly use the Thomas-Fermi and Whittaker solution to describe the central and outer profile of the star. This is shown in Fig. 2 for γ = 10 (red), 100 (blue) and 1000 (green). The upper panel shows the wave function profile s(z) normalized by its central value s 0 . The numerical solution is shown in gray, the Thomas-Fermi solution for the central profile given by Eq. 31 as dashed lines and the Whittaker solution given by Eq. 26 as solid lines. The lower panel shows the accuracy of the approximate solutions with respect to the numerical solutions.
We can see that for large couplings, a combination of the Thomas-Fermi and Whittaker solutions describe profile with an accuracy of O(10 −6 ) in the center and outside the star, corresponding to the accuracy of the expansion parameters. However, even for large couplings, there remains a surface region around z ≈ Z that is poorly described by both the Thomas-Fermi and Whittaker solutions. Even though the wave function, and therefore the mass density, is small in the surface region, the fraction of the star's mass contained in it can still be sizeable.
To obtain a better description, in particular for the surface region, we can use the series expansion derived in Sec. III A. This is shown in Fig. 3 for an intermediate

A. Using the Solution
In the above discussion, we have solved the Gross-Pitaevskii-Poisson equations in Eq. 20, expressed in terms of the dimensionless variable s and v. We have seen that the solution of these equations is fully parametrized by the dimensionless coupling parameter γ and then obtained an analytic form for profile of boson stars. To use the solution to describe a boson star with mass M made of boson with mass m and self-coupling λ, we therefore need to obtain the corresponding value of the dimensionless coupling parameter γ.
Combining our scale choice k = GM m 2β with Eq. 15 and 19, we can see that Note that the value of the self-coupling parameter γ is independent of the boson mass m and only depends on the star's mass M and the quartic self-coupling λ. We can obtain the coupling parameter γ by solving for γ as a function of ξ where we have introduced the short-hand notation ξ = λGM 2 . Alternatively, we can follow our previous approach used to obtain the expansion coefficients and obtain an expression for γ by finding a suitable fit to the numerical solution. This is shown in Fig. 4. The red dots correspond to the numerical solution using the same data set as Fig. 1.
In the Thomas-Fermi limit, the expansion parameter β is given by β = π 2 √ γ. We can solve Eq. 41 and obtain as shown in Fig. 4 as dashed cyan curve. To obtain a more accurate result, we follow our approach used to obtain the expansion parameters and fit the numerical solution with a suitable series expansion. For large values of ξ > 100 we expand around the Thomas-Fermi limit and obtain with ξ min = −51.523602. Note that domain of ξ is restricted to ξ > ξ min , and therefore, the coupling parameter γ is bound from below γ > γ min = −0.722. Physical boson stars therefore need to fulfill the condition λGM 2 > ξ min . This implies that the maximal mass of a boson star with attractive self-interaction, λ < 0, is given by There exists a second branch of the solution corresponding γ < γ min . However, the corresponding boson star would have a higher total energy than the solution for γ > γ min . Such configuration is unstable with respect to perturbations and therefore unphysical [36]. The accuracy of the expressions given in Eq. 44 and Eq. 43 with respect to the numerical solution is shown in the lower panel of Fig. 4. We can see that the analytic expressions reproduce the numerical results with accuracy better of O(10 −4 ) for all physical values of ξ.
Knowing the value of γ for a given boson star, we can obtain the profile for physical wave function ψ and grav- itational potential Φ by rescaling the dimensionless solution s and v according to Eq. 15 and 19. We can write Furthermore, we can simply read off the binding energy and the central density

B. Axion Stars
In the following, we illustrate the use of the obtained solution on one particularly well-motivated scenario: axion stars. The axion is a real pseudo-scalar field, which was initially introduced to solve the strong-CP problem.
Furthermore it also provides a natural dark matter candidate if its mass in the range between m = 10 −5 eV and 10 −2 eV. The axion potential can heuristically be described by the instanton potential where a denotes the axion field, m π = 135 MeV is the pion mass, f π = 92 MeV is the pion decay constant and f is the axion decay constant. After expanding the potential and matching it to the form of Eq. 1, we find that the axion mass and decay constant are related by mf = m π f π . The quartic coupling is given by where the negative sign indicates that the self-interaction is attractive. Note that there exist higher terms in the expansion. However, those terms will not contribute in the non-relativistic approximation and have been shown to be insignificant for dilute axions stars even in the relativistic limit [22].
In the previous section, we have seen that a boson star with attractive self-interaction has an upper mass. Using the axion parameters this implies We have seen in Sec. II B that the non-relativistic approximation is valid for axion stars with mass For axions in the dark matter axion window of masses between m = 10 −5 eV and 10 −2 eV, the upper bound on the axion star's mass is well below this limit, and therefore, the axion star is well described by the nonrelativistic approximation.
In Fig. 5 we show the binding energy e given in Eq. 47 as a function of the axion star's mass M for a axion mass m = 10 −5 eV. The green line shows the result for a non-interacting axion field, λ = 0. In this case the binding energy is given by e = −0.1627 · G 2 M 2 m 3 , which increases quadratically with the star's mass. The red lines indicates the axion field including the attractive self-interaction, λ = − m 2 2f 2 = −3.2 · 10 −53 . For small axion star masses M 10 −12 M , it approaches the noninteracting limit. In this case the classical pressure due to the self-interaction is negligible since the axion density is small even in the center of the star. For masses M > 10 −12 M the self-interaction term becomes important and the binding energy increases relative to the non-interacting case. The upper bound on the axion mass star mass, M max = 2.74 · 10 −12 M , is indicated by the gray dashed line. The authors of [46] have shown that FIG. 6. The upper panel shows the energy density profile ρ for an axion star with mass M = 2.5 · 10 −12 M and axion mass m = 10 −5 eV. We show the results for an axion-star without self-interaction (green), with attractive self-interaction (red) and a wrong-sign axion-star with a repulsive self-interaction (blue). We use the truncated series expansion of the wave function at small radius with N = 20 (dashed) and with at large radius with N = 3, M = 2 (solid) as given in Eq. 27. The lower panels shows the corresponding differential mass distributions dM/dr. axions with mass m = 10 −5 eV will form axion miniclusters with masses around M ∼ 10 −12 M and therefore in the regime in which self-interactions are important.
For comparison, we also show the results for a wrongsign axion star, in which the self-coupling has the same magnitude as the axion field but is repulsive, λ = + m 2 2f 2 = +3.2 · 10 −53 , as a blue line. At high masses, the binding energy approaches the Thomas-Fermi limit, shown as a dashed cyan line. In this case the binding energy is given by e = − 4 √ λπ G 3 2 M m 3 . In the upper panel of Fig. 6 we show the density profile, ρ(r) = M ψ 2 (r). For illustration, we choose an axion mass m = 10 −5 eV and axion star mass M = 2.5 · 10 12 M which is slightly below the maximal mass M max . Following the previous discussion, we consider an axion-star without self-interaction (green), including the attractive self-interaction (red) and a wrong-sign axionstar with a repulsive self-interaction (blue). To calculate the density profile, we use the truncated series expansion as given in Eq. 27 with N = 20 at small radius, as indicated by the dashed lines, and N = 3, M = 2 at large radius as indicated by the solid lines. We can see that the two series expansions match well at intermediate values of the radius. Note that the attractive self-interaction of the axion field leads to a significant deformation of the density profile with respect to the a non-interaction axion star of same mass. This indicates the importance of including the axion self-interaction when considering axions stars close to their maximal mass M ∼ M max .
The lower panel of Fig. 6 shows the differential mass distribution dM/dr = 4πr 2 ρ for the three considered stars. We can see that the mass distribution peaks at intermediate radii and approximately half of the star's mass is described by both the small and large radius expansion of the wave function as given in Eq. 27. This once again shows the importance of an accurate descriptions of the tails of the wave function profile in order to correctly describe the axion star.

VI. CONCLUSION
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 a previous study [21], we presented a semi-analytic solution to these equations describing the profile of boson stars formed by a noninteracting scalar field. The solution is based on a series expansion which is parametrized by four expansion parameters that have been obtained from numerical simulation at high accuracy.
In this study, we have generalized our semi-analytic approach to boson stars where the constituent particles have self-interactions. In this case, the expansion param-eters are functions of the quartic self-coupling. Based on results from numerical simulations, we found a corresponding analytic expression for all expansion parameters.
This allows to simply obtain profiles of boson stars in an analytic form for arbitrary self-couplings at high precision directly from the series expansion. In particular, no further time consuming and computational expansive numerical integration is needed.
We have also applied our methods to axion stars, and shown how the the mass and density profiles can be obtained for both weak and strong interactions. The profiles are significantly modified from the case of noninteracting bosons. The methods developed in this paper allows for systematic studies of the properties of boson stars in an analytic way without further relying on numerical simulations.
Finally, we note that there are several possible generalizations of these results. In particular, we can extend our results to rotating boson stars. It would also be interesting to see how the profiles are modified in the presence of other astrophysical objects like planets. We leave these and other questions to future work.