On Approximation Methods in the Study of Boson Stars

We analyze the accuracy of the variational method in computing physical quantities relevant for gravitationally bound Bose-Einstein condensates. Using a variety of variational ans\"atze found in existing literature, we determine physical quantities and compare them to exact numerical solutions. We conclude that a"linear+exponential"wavefunction proportional to $(1 + \xi)\exp(-\xi)$ (where $\xi$ is a dimensionless radial variable) is the best fit for attractive self-interactions along the stable branch of solutions, while for small particle number $N$ it is also the best fit for repulsive self-interactions. For attractive self-interactions along the unstable branch, a single exponential is the best fit for small $N$, while a sech wavefunction fits better for large $N$. The Gaussian wavefunction ansatz, which is used often in the literature, is exceedingly poor across most of the parameter space, with the exception of repulsive interactions for large $N$. We investigate a"double exponential"ansatz with a free constant parameter, which is computationally efficient and can be optimized to fit the exact solutions in different limits. We show that the double exponential can be tuned to fit the sech ansatz, which is computationally slow. We also show how to generalize the addition of free parameters in order to create more computationally efficient ans\"atze using the double exponential. Determining the best ansatz, according to several comparison parameters, will be important for analytic descriptions of dynamical systems. Finally, we examine the underlying relativistic theory, and critically analyze the Thomas-Fermi approximation often used in the literature.


I. INTRODUCTION
The recent surge in the study of scalar field condensate dark matter (DM) is in part driven by the failure to detect individual Weakly Interacting Massive Particles (WIMPs) at the Large Hadron Collider, or in various direct detection experiments. There is another avenue for dark matter to manifest, as condensates of macroscopic size [1][2][3][4][5][6][7][8][9]. Electrically neutral boson particles, if they are a component of dark matter, can naturally form gravitationally bound Bose-Einstein condensate (BEC) bubbles below a critical temperature, due to quantum statistical effects. These entities are known either as oscillons or boson stars. It is quite interesting and important to understand the physical properties of such condensates to determine if they are a viable alternative to the more popular and widely investigated WIMPs. Thus it is important to develop analytical and numerical methods to study condensate formation, as well as their stability, evolution, and possible decay into component particles.
One widely studied example of scalar field DM forming BECs is the axion. Axion stars were considered first around 30 years ago, originally suggested to form from collapse of overdense miniclusters in the early universe [10,11] (see also more recent simulations [12]). Since then, many properties of axion stars have been studied extensively; these include structural stability [13][14][15][16][17][18][19][20][21][22] (including nonzero angular momentum [23][24][25]), the process of gravitational collapse [26][27][28][29][30][31][32], and their decay through emission of relativistic particles [33][34][35][36][37][38]. There has recently been a significant amount of work regarding relativistic corrections more generally to the classical field description of axion stars [39][40][41][42]. Other authors have investigated possible connections to astrophysical radio sources [43,44]. In some scalar field models, boson stars can be extremely heavy and (if they are stable) could give rise to gravitational wave signatures [45][46][47][48][49][50]. Clearly this is a field booming with new and interesting results. (For boson star reviews, see e.g. [51][52][53].) Besides DM, there are other classes of problems in cosmology where condensate formation is important. Elementary Hermitian boson fields known as inflatons are postulated to drive cosmological inflation, the hypothesized early-time exponential expansion of the universe. In addition to this, a bosonic degree of freedom (termed quintessence) is posited to generate the energy density which causes the observed late time acceleration of the universe. Both inflatons and quintessence can fragment and form BEC bubbles in the same way as dark matter candidates discussed previously [54][55][56]. If sufficiently long lived, these entities can play a crucial role during inflation and also at later stages of cosmological evolution.
Theoretical studies of Bose-Einstein condensation gained prominence in 1990's, after the experimental discovery of atomic BECs in systems of cold atoms [57][58][59]. Atomic condensates are described by the Gross-Pitäevskii (GP) equation, which is a form of the nonlinear Schrödinger equation. The interatomic interactions could be attractive or repulsive, and the atoms could be placed on external potentials. Analytic solutions of the GP equation are difficult to obtain, and various approximation methods had to be employed. Static problems were more amenable to numerical studies, but dynamical questions like expansion, collapse, and decay of condensates required the employment of approxima-tion methods.
The variational method is one approximation technique widely utilized in studies of atomic condensates. Recently, the variational method was also adopted for the study of gravitationally bound condensates of bosons by Chavanis [14]; he used a Gaussian ansatz to approximate the wavefunction and compared the results he obtained thus with numerical solutions of the GP equation [15]. This comparison of ansätze to the numerical solutions is imperative for the study of dynamical problems that are much more difficult to solve numerically. Some such dynamical problems that have been analyzed using the variational method are the collapse [26][27][28][29] and collisions [60,61] of BECs. In subsequent years, numerous authors have presented various ansätze for both static and dynamical problems, either to improve numerical agreement or computational efficiency [19,22,25,37]. 1 To our knowledge, the relative efficacy of one ansatz compared to another has not previously been considered in a rigorous fashion. Some ansätze are simpler computationally, others are more complicated. Some match numerical results for boson star masses but do not reproduce as well the radius. As the literature on boson stars becomes more complex, it becomes increasingly important to have at hand a tractable analytic approach, and to understand the benefits and weaknesses of different choices of approximate wavefunctions. This is the basic goal of the current project.
In this work, we will provide background information for both the time-independent and time-dependent variational methods (Section II). The latter is crucial for understanding dynamical processes, like collapse and decay. We will then analyze various classes of approximate wavefunctions used in the literature, providing comparisons across different ansätze (Section III) as well as comparison to exact numerical results (Section IV). In this discussion, we will also propose a new ansatz with a free parameter, which can be varied to give excellent agreement by whichever measure is relevant to a particular scenario. Finally, we compare the nonrelativistic formulation to the underlying relativistic theory in Section V, and examine where the Thomas-Fermi approximation (often used in the literature) breaks down. We conclude in Section VI.
We will use natural units throughout, where = c = 1.
1 A third approach, in which the exact wavefunction can be computed using an analytical expansion, was developed in [62,63]. This approach has the dual advantage of arbitrary precision and analyticity, though it is still more computationally intensive than the variational approach.

A. General Formulation
The variational method [64] for finding approximate analytic solutions to eigenvalue problems was utilized as an important calculational tool in the early development of quantum theory. Its success in describing the ground state of Helium [65] played an important role in establishing modern quantum mechanics as a viable theory.
As explained in quantum mechanics texts, the variational method involves the extremization of the expectation value of the Hamiltonian of the system Ψ|H|Ψ with respect to variations of a class of candidate wave functions |Ψ . In the calculation, the norm Ψ|Ψ is held fixed, and in the end one obtains approximate analytic solutions to the time-independent Schrödinger equation. A corresponding method to obtain solutions to the timedependent Schrödinger equation is to extremize the action with respect to a class of time-dependent variational wave functions |Ψ(t) . This time-dependent variational method, first introduced by Dirac in 1934 [66], is not usually described in standard texts of quantum mechanics, but it has found wide use in fields like nuclear physics [67], chemistry [68], and quantum field theory [69].
In condensed matter discussions, the GP formalism is developed starting from an N -particle wave function where each particle is in the state φ( r). Following the derivations in [59], the expectation value of the Hamiltonian in the state χ ( r 1 , ..., r n ) = i φ( r i ) yields the energy functional where V ( r) is an external potential, U 0 represents some contact interaction among the scalars, and φ is normalized to unity. Defining ψ = √ N φ( r) and taking N 1 yields, Extremizing the above equation, subject to the normalization constraint (II.5) where µ, the Lagrange multiplier introduced to maintain the normalization, is defined as the chemical potential.
The constant c = 1 if V ( r) does not depend on ψ * (true for an external trap), but can take other values elsewhere (e.g. for condensates bound by self-gravity). By multiplying eq. (II.5) by ψ * and performing a volume integral, we get Comparing eqs. (II.3) and (II.6), we can conclude the chemical potential is not equal to the energy per particle given in eq. (II.3) if the interaction strength U 0 = 0 or if the trapping potential V ( r) depends on ψ * (i.e. if c = 1).
In the time-dependent case, the GP formalism starts with the action S = dtL, where The Hamiltonian is which is of course identical to the energy functional defined in equation (2.3). Variation of eq. (II.7) with respect to ψ * gives the Gross-Pitäevskii (GP) equation which reduces to the time-independent case of eq. (II.5) if the time dependence of ψ is given by the factor exp(−i µ t). This is also in agreement with the the identification of µ as the chemical potential. The variational method has been applied to study time-independent as well as time-dependent solutions to the GP equation. As we will describe in Section III, in this case one describes the wavefunction by some ansatz ψ( r), so that, for a given physical system, eq. (II.3) can be extremized analytically. Of course, for many physical applications the solutions to the time-independent equations can be found numerically to an arbitrary precision; therefore the variational method offers qualitative understanding of the systems, but falls short of the numerical accuracy.
On the other hand, dynamical problems are more difficult to solve using numerical methods. Time-dependent solutions to the GP equation can be approximated using the variational method formulation of Pethick and Smith [59], which we have described above. A good example of such an application is the collapse of BECs [26][27][28][29]. In this case, the variational parameter σ is taken to be a function of time, σ(t), while the wavefunction is multiplied by some phase. This phase depends on the velocity field of the system, which for spherical collapse is taken to be in the radial direction and proportional to r. It also depends on some parameter (analogous to the Hubble parameter in cosmology) which is related to σ(t) by H(t) =σ(t)/σ(t). It is remarkable that using a simple analytic description of this type, one can assess many of the relevant features of BEC collapse. This illustrates the power of the time-dependent variational formalism.

B. Application to Boson Stars
Dynamics of boson stars are described by the Klein-Gordon equation in the presence of self-gravity. This action, (II.11) in the nonrelativistic limit and with the replacement The gravitational potential originates from the nonrelativistic limit of the Einstein field equations with G = M P −2 Newton's gravitational constant. (See Section V for a more thorough description of the nonrelativistic limit.) In this work, we will consider the case of a quartic self-interaction, where self-coupling λ can be positive (giving rise to a repulsive interaction) or negative (attractive interaction). Variation of eq. (II.12) with respect to ψ * yields a GP-type equation comparable to eq. (II.10). (II.14) Finding analytic solutions to the above non-linear equations (II.13) and (II.14), known as the Gross-Pitäevskii+Poisson (GPP) system, is a challenging task. In order to find an approximate analytic solution one could extremize the action S within a class of variational wave functions as described previously in condensed matter/atomic physics applications.
Assuming a harmonic time dependence for the wavefunction, the chemical potential of eq. (II.6) is given by while the energy of eq. (II.3) is The latter can be extremized by assuming some r dependence for ψ, giving rise to an approximate bound state solution. We present a number of prevalent choices for such an ansatz in the next section. We will focus here on the time-independent case, because it will allow us to analyze the radial dependence of the wavefunction; any relevant time-dependent factors will depend on the particular application one considers (e.g. collapse) and can be added on later.
In order to begin, we will need to impose some conditions on the classes of wavefunction ansätze we consider. Assuming that the wave function ψ(r) is differentiable to all orders at r = 0 and the gravitational potential V g has a Laurent expansion in r, the GPP equations (II.13) and (II.14) imply that the gravitational potential is an analytic function in the variable r 2 at the origin; consequently, the first derivative of the wave function ψ(r) vanishes there. At large r, the interaction term in the GP equation becomes negligibly small, and the GP equation reduces to the linear Schrödinger equation for a particle moving in the gravitational potential V g . The potential at large r takes the GM/r Newtonian form, and the wave functions asymptotically are nearly hydrogen-like, with a exp(−kr) behavior (though see below). We impose these boundary conditions when numerically solving the GPP equations. A successful variational ansatz should also exhibit similar behavior.
Actually, the behavior of the wavefunction at large r is not exactly exponential, but can be calculated in a straightforward way. First, at large r we throw away the terms in the equation of motion that are higher order in the fields than the gravity term; this gives eq. (II.14) with λ = 0 and V g = −GM m/r. This equation has an exact solution in terms of hypergeometric functions, but the leading order solution is proportional to This is in agreement with the results of [62,63].
Compare and contrast this case to that of atomic Bose condensates trapped by external harmonic oscillator potentials. In the latter situation, the GP wave function at large r, when interaction terms are negligible, should approach Gaussian form, making such functions ideal candidates for variational ansätze. This makes clear the fact that the optimal ansatz for the wavefunction should depend on the potential.

A. Generalized Ansatz
In [28,29] we performed an analysis using a general ansatz for the wavefunction of an axion star, under only the (weak) assumptions that the wavefunction is finite at the origin, spherically symmetric, and decreases monotonically with the radial coordinate r. Here, we apply the same method to a general boson star with a 4-point coupling λ.
The starting point is the rescaling of the macroscopic quantities where σ is a variational distance parameter explained in detail in the next section. Then a general ansatz for a boson star can be written in the form where at fixed ξ the function F (ξ) is independent of ρ, n, and λ. Because we can fix F (0) = 1 without loss of generality, we are able to identify w = ψ(0), the central value of the wavefunction. Substituting the ansatz into the normalization condition (II.4) gives the central value as where we introduced the notation Using eq. (II.16) and the general ansatz of eq. (III.2), we obtain for the energy functional and, using eq. (II.15), we obtain for the chemical potential For simplicity, we have defined the dimensionless parameters (III.7) Recall once more that λ > 0 (λ < 0) will correspond to repulsive (attractive) self-interactions. Because the GP equation can be derived from the variation of the total energy, a solution will be a stationary point of (II. 16). Given an ansatz for the wavefunction, we approximate this exact solution by minimizing eq. (III.5) with respect to ρ. This procedure gives the dilute boson star radius wheren is an important scaled particle number given bȳ (III.10) For λ < 0,n determines the mass of the maximum stable configuration, and so we will denote it byn = n c . For λ > 0, there is no maximum mass in the nonrelativistic limit, so we will instead use the notationn = n * ; in this case, self-interactions become increasingly important for n ∼ n * , eventually approaching the region in which the Thomas-Fermi approximation is relevant. We will discuss these points in detail below. Eq. (III.9) is a stable minimum of the energy classically and a metastable solution in quantum theory. There exists another root of the energy at a radius of which we will discuss later. In the case of attractive selfinteractions, ρ u > 0 is an unstable maximum of the energy; for repulsive interactions, ρ u < 0 is an unphysical configuration. Above, we have used the formulation of our previous works [28,29], in which the dimensionful parameters are scaled out; then the numerical constants B 4 , C k , and D 2 are dimensionless and depend only on the shape of the ansatz being employed. Other authors, e.g. [14,70,71], use a different formulation; for ease of comparison we provide the expressions for translating between the two in Appendix A.
This formulation has a number of useful applications. For example, in the limit of no self-interactions (λ → 0) 2 , there is a simple relation between the energy and 2 Of course, in the λ → 0 limit the scaling of eq. (III.1) is not appropriate. In that case one should use σ = ρ/m and N = M P 2 n/m 2 , which gives rise to the analogous equations for E and µ without the interaction term. In the end, one finds the first two terms in eqs. (III.5) and (III.6) do not change, which is all that is needed here. chemical potential: which was derived using other methods in [14,72].
For attractive interactions, we of course recover the standard result that the local minimum and maximum of the energy become equal at n = n c . For n > n c , no stable solutions exist. The value n = n c also corresponds to a critical minimum radius, It is useful to define a parameter which parameterizes the closeness of N to N c . Then substituting ρ d from eq. (III.9) into eqs. (III.5) and (III.6), the energy and chemical potential at the minimum for attractive self-interactions can be written in the form This implies a simple relationship between the energy and the chemical potential that is independent of the choice of ansatz, which is an analogue of the non-interacting result in eq. (III.12), but applied to attractive self-interactions. We can see that when δ → 0 (N → N c ), the ratio goes exactly to 5. Even though the derivation assumed some ansatz for the variational approach, the result does not depend on what form the wavefunction takes, and so it holds even in the exact case. In the case of repulsive self-interactions, there also exists a critical mass, though it is not at n = n * as defined above; it arises due to relativistic effects which we do not consider here [4]. Relativistic effects can be taken into account using the Ruffini-Bonazzola (RB) formalism for analyzing boson stars [2]. One could in principle formulate a variational method which approximated the relativistic equations of motion, in which case these effects would become apparent. We leave such an analysis for future work. Of course, at weak gravity and small binding energy, the RB equations of motion reduce to the GPP system, as we will describe in Section V.

B. Non-Compact Ansätze
It is well-known by direct solution of eqs. (II.13) and (II.14) that the wavefunction of a boson star does not typically have compact support; it is nonzero for all r ≥ 0 though decreases extremely fast at large r. The standard definition for the "size" of a boson star is R 99 , the radius inside which 0.99 of the mass is contained. A number of non-compact ansätze have appeared in the literature to approximate the exact solution; a few of the most popular ones are [14, 25, 27-29, 37, 61]) N π σ 3 e −r/σ (Exponential (E) [19,22]) [22,37]) On the right we give the long name (e.g. "Gaussian") and short abbreviation (e.g. "G") for each ansatz, and cite a collection of previous works which utilize them in the study of boson stars. We will use the notation that, for example, ψ G (r) is the Gaussian wavefunction, ψ LE (r) is the linear+exponential wavefunction, etc. The first four functions listed here are popular choices in the literature; the last one (the double exponential) is a proposal of ours with some constant parameter a which can be fixed by matching to the exact solutions. We will show in this work that the double exponential can be optimized for a given numerical result. Note the appearance of σ in each of the ansätze. While this parameter has units of distance, it should not be confused with the radius of a boson star. For each ansatz, the parameter σ d of the solution is related by some constant factor to the radius of the boson star R 99 ; we will label this by a real number κ, i.e.
where both σ d and κ depend on the ansatz under consideration. Other reasonable distance scales in a calculation like this one include the expectation values which will be useful for comparing to exact solutions later. Because the translation to a physical length is different for a different ansatz, the parameter σ means something different depending on the ansatz in which it is employed; said a different way, σ itself is unphysical and should not be compared across ansätze.

C. Compact Ansätze
It is sometimes advantageous to employ a compact function to approximate the wavefunction. An example of such a case is the Thomas-Fermi (TF) limit of repulsive interactions, where one neglects the kinetic energy term, and the resulting wavefunction is very close to having an exact radius R. We also found a compact ansatz advantageous when discussing collapse of an axion star through collisions with astrophysical sources [61]. Here, we include a few compact ansätze for completeness: However, we do not include them in our analysis, partly because it is somewhat problematic to compare R 99 for non-compact ansätze to the exact radius for compact ansätze. It should also be noted that one must require a compact ansatz to vanish above r = R. The TF ansatz given above is the exact solution to the GPP equations in the Thomas-Fermi limit [74]. However, given that the solution has an infinite derivative at r = R, we do not recommend its use as an ansatz. To investigate the efficacy of these ansätze further, we have to solve the full equations of motion (II.13) and (II.14) in the stationary limit. We employ the following scaling to dimensionless ("tilde-d") variables: Writing the equations of motion in terms of the rescaled quantities, we solve the resulting system numerically: We employ a shooting method to determine correct boundary conditionsψ → 0 andṼ → constant as r → ∞. In practice, in the numerical routineψ and V converge up to some finite radiusr 0 , which can be as large as the precision of the calculation requires. After solving forψ andṼ , one can calculate the dimensionless macroscopic quantities using n = Before describing the solutions, note that the GPP energy can be written in terms of tilde-d quantities, so that, after dividing by N and m, we get Of course, because µ appears on both the LHS and RHS, this equation can be used as a consistency check on the numerical solutions. We have verified for every numerical solution that the eigenvalue µ satisfies this constraint. The important physical quantities describing a boson star are its mass and radius. In Figure 1, we show the exact numerical relationship between these quantities, along with the result from various ansätze. For attractive interactions, we observe that ψ E tends to be a good fit for small R 99 (on the unstable branch of solutions) but a poor fit for large R 99 (the stable branch). It is also interesting to note that ψ G is used quite often in the literature, but is a poor fit for small R 99 and badly approximates the position of the maximum mass. The double exponential ansatz has a free parameter a, and we find that we can match the value of n c almost exactly by taking the value a = 3π/8.
One interesting feature that is apparent in the right panel of Figure 1 is that, in the case of repulsive interactions, the radius at large M approaches a constant. This is the limit of the Thomas-Fermi approximation, which will be discussed in Section V. When using the variational approach, the radius R T F to which a given ansatz approaches at large M in the repulsive case is equal to R c , the radius at which one obtains the maximum mass in the attractive case. (This is clear by examination of eq. (III.9) as well.) This relationship between these two radii in very different limits is a consequence of the fact that we use the same form for the ansatz in both the repulsive and attractive interaction cases. As a result, it does not hold true for the numerical solutions.
It is interesting that the result for the double exponential ansatz almost exactly approximate the sech and linear+exponential ansätze for the choices a = 2 and a = 3π/8, respectively. We will discuss this in more detail in the following section.

B. Fitting the Double Exponential
Every variational calculation, regardless of the ansatz, can be improved by the addition of an extra parameter. The double exponential ansatz proposed in this paper has an extra variational parameter that can be tuned to fit the numerical solution and is computationally efficient. The extra variational parameter allows one to choose how to fit the ansatz to the numerical solution. For example, one can choose to fit the ansatz and numerical wavefunctions, or to fit various expectation values. In the comparisons above, we found that the choice a = 3π/8 reproduced the exact value of the maximum mass. This choice also resulted in an ansatz, denoted ψ DE 3π/8 , that was in very good agreement with ψ LE . We also considered the variational parameter a = 2, which resulted in an ansatz ψ DE2 that was in good agreement with ψ S . By varying the parameter a, we found we could optimize agreement to the numerical solutions on both branches of attractive solutions, as well as for repulsive interactions.
As noted previously, ψ DEa can be tuned to fit ψ LE and ψ S . We can explain this by looking at the double exponential in different limits of a. First of all, it is clear that in the limit a 1, we have ψ DEa → ψ E , the ordinary exponential ansatz. On the other hand, in the limit a → 1, where f (a) is the original prefactor in eq. (III.17). Because the prefactor in any ansatz is finite and equal to ψ(0) by definition, we find that ψ DEa → ψ LE in the limit a → 1. In this sense, the double exponential ansatz interpolates between ψ E and ψ LE as a varies between 1 and 1. Slightly more perplexing is the sech ansatz, which seems to be reproduced approximately when a = 2. The reason for that is that the second derivatives at r = 0 coincide for those ansätze. The behavior near the origin is the most important for all integrals over the wave functions, because the largest contributions to all integrals come from that region. This explains this coincidence and further illustrates the versatility of ψ DEa .
The double exponential can be easily further generalized, as it is straightforward to add additional exponential functions with further fitting parameters a 1 , a 2 , a 3 , .... Because computations on exponential functions are relatively fast, these additional functions do not increase the computation time very much (also see next section). Of course, one could generalize any of the other ansätze by introducing a fitting parameter, but the increase in computation time is more problematic for these other ansätze.

C. Methods for Quantitative Comparisons
In this section, we report on the direct comparison of the exact numerical solutions of the previous section with    Table I, for each ansatz. The top row are the absolute times in seconds on a standard laptop; the lower row is the time normalized to the Gaussian case.
the ansätze we presented in Section III. One criterion on which one can compare different ansätze is computational efficiency; in our analysis we have seen, for example, that the sech ansatz is significantly more difficult to employ than a function with exponential dependence. This is an important consideration because the purpose of using an ansatz in the first place is to simplify the calculation. The sech ansatz, for example, is orders of magnitude slower than the others (and, as we will show below, it does not pay off sufficiently in numerical accuracy). The double exponential is relatively fast, regardless of the value of a. Another criterion which is potentially more important is proximity to the exact solution. To make such comparisons simple to understand, we will rewrite the ansätze slightly so that the radial variable is the same as the numerical one, given in eq. (IV.1); that is, we write where A labels the ansatz being considered, as before. Theñ can be compared directly to the numerical resultψ as a function ofr. A simple comparison of the shapes of the wavefunctions can be performed as follows: for a given N , compute r , r 2 , and R 99 for the ansatz and compare it to the exact value computed numerically. These three radial variables should be a reasonable proxy for the shape of the wavefunction. More precisely, we will compute (IV.12) A potentially more robust method of comparison is to evaluate the numerical solution and a given ansatz at the same value of N , and see how large the wavefunction deviations are over the whole range of r. To make this quantitative, we compute the difference integral (IV.13) We will use both ∆ ψ and ∆ r to compare how well a given ansatz matches the exact result. In the following sections, we will show the result of these analyses for attractive interactions on both the stable and unstable branches, as well as for repulsive interactions.

D. Attractive Interactions: Comparison for Stable Branch
We begin with an analysis of the stable branch of solutions for attractive self-interactions. This branch is relevant for so-called dilute boson stars, including axion stars. In the top row of Figure 2, we show the exact and approximate wavefunctions, with the vertical axis on a log scale. The top row correpsonds to attractive interactions along the stable branch. The most commonly used ansätze are the exponential and Gaussian functions, which turn out to be the worst fits to the numerical solutions, and are particularly bad in the large-r tail. The other ansätze do reasonably well both near the core and in the tail.
Near the maximum mass, both the Gaussian and the exponential ansätze approximate the exact solution exceptionally badly. This is in part clear from Figure 1: because these functions do not well-approximate the exact value of the maximum mass, the radius will also be very different, resulting in larger overall deviations (we will quantify these statements in a more precise analysis below). The other ansätze provide a much better fit in this region of parameter space. At larger, the double exponential does slightly better than the linear+exponential function.
For a more quantitative comparison, we turn to the bottom row of Figure 2. At every n, the smallest ∆ ψ deviations are found using the LE ansatz. In the limit a → 1, the same is found for the double exponential function, as discussed in an earlier section. In fact, we have scanned over values of the parameter a and found that the best agreement to numerical results (i.e. smallest ∆ ψ and ∆ r ) is obtained in this limit. In this sense, the linear+exponential function provides the best fit to the data, and is also computationally efficient. It is also the most similar to the leading-order wavefunction of [62,63], whose wavefunction is explicitly calculated to have the correct asymptotic behavior as r → ∞.

E. Attractive Interactions: Comparison for Unstable Branch
We move now to an analysis of the unstable branch which is relevant for collapsing boson stars, an important topic both for ordinary QCD axion stars [28,29], as well as possible galaxy core collapse induced by as-trophysical interactions [38]. In the top row of Figure 3, we show the exact and approximate wavefunctions, with the vertical axis on a log scale for attractive interactions along the unstable branch. It can be seen that for larger, ψ G is a poor fit for both small and large masses, while ψ S and ψ DE2 tend to be better fits than ψ LE and ψ DE 3π/8 . It is also interesting to note that ψ E does better than all ansätze for some range of n, though we do not have a satisfactory explanation for this behavior at this time. We do observe in the numerical solutions that the wavefunctions become increasingly "squeezed" towards r = 0 on the unstable branch, and while the wavefunctions remain cored, they more closely approach an exponential behavior closer to the origin.
In the bottom row of Figure 3, we show the ∆ ψ and ∆ r deviations for the unstable branch. For a range of masses, ψ E does better than all ansätze, while ψ S and ψ DE2 do better than ψ LE and ψ DE 3π/8 for all masses except near the maximum mass. We conclude that over the most relevant parameter space, an exponential function is the ideal choice for analysis on the unstable branch, though close to n c it is preferable to use the double exponential with a = 2. We do not recommend use of the sech function, only because of the computational inefficiency pointed out in Section IV C.

F. Repulsive Interactions: Comparison
Finally, we move to an analysis of repulsive interactions which is relevant for axion-like particles; for recent model-building of dark matter scalar with repulsive interactions, see [75]. For such self-interactions, a maximum mass arises due to relativistic effects, to which our nonrelativistic analysis is not sensitive. This maximum mass [4] M rep max = 0.22 is much greater than the masses resulting from our numerical solutions, and so the comparisons made here are all for physical masses. For repulsive interactions, as the mass of the condensate increases, it approaches the Thomas-Fermi limit at which its kinetic energy becomes negligible and its radius becomes independent of mass [73]. From eqs. (III.9) and (III.10), one can see that as the number of particles n increases, the solution for the equilibrium radius does, indeed, approach a constant. This constant radius happens to be the critical radius for attractive self-interactions (eq. III.13), given that the ansatz considered is held constant. It is also interesting to note that N = N * can be seen as a scale at which the condensate begins to approach the Thomas-Fermi limit.
In the top row of Figure 4, we show the exact and approximate wavefunctions, with the vertical axis on a log scale. One can see that for n n * all ansätze become increasingly worse fits in the tail. As is the case for attractive self-interactions, ψ LE and ψ DE 3π/8 do better than ψ S and ψ DE2 . It is interesting to note that for n/n * ∼ 40, the wavefunction of the numerical solution falls off more rapidly at larger than ψ G . This is in agreement with the results of Böhmer and Harko [74] who found an exact solution to the equations of motion for the Thomas-Fermi limit to be of the form sin(ξ)/ξ. This exact solution has been included as a possible compact ansätze to be used for axion-like condensates with repulsive interactions and large numbers of particles (see Section III C). For N N c , the condensate approaches the noninteracting limit; this is why the top-left panel of Figure 4 is extremely similar to the top-left panel of Figure 2.
In the bottom row of Figure 4, we show the ∆ ψ and ∆ r deviations for repulsive interactions from which one can see that ψ G does better than all other ansätze at large n. It is also interesting to note that the deviations have a considerable change in behavior near n * (most notably for ψ G ), which is presumably due to the condensate approaching the Thomas-Fermi limit. For repulsive interactions at large n, the Gaussian function is thus the appropriate choice.

V. RELATIVISTIC FORMULATION
In the preceding sections, we have analyzed various approximations to the GPP formalism. However, even the exact GPP equations are themselves an approximation to the underlying relativistic theory. The theory at high energies is defined by the Klein-Gordon equation for a scalar field, which is coupled to a curved spacetime metric described by the Einstein equations. In this section, we describe precisely the limit of the Einstein+Klein-Gordon (EKG) equations in which one recovers the GPP equations. We will then apply the relativistic formulation to the case of repulsive interactions, to analyze the limitations of the TF approximation.

A. Equivalence between GPP and EKG Equations
A relativistic method of describing boson stars was pioneered by Ruffini and Bonazzola (RB) [2]. The crucial idea is to take the field A describing the boson star to be linear in ground-state creation and annihilation operators a ( †) 0 , A(r, t) = R(r) e −i µ0 t a 0 + e i µ0 t a † 0 , (V.1) and use these operators to build N -particle states defined by Here, µ 0 is the energy eigenvalue of the relativistic equations; it is related to µ in the GPP formulation by which effectively defines the particle number N .
The gravitational metric describing the curvature of spacetime in the presence of the boson star is which is spherically symmetric, as is the ground state wavefunction R(r) in eq. (V.1). Then the rr and tt components of the Einstein equations (taken as expectation values), along with the Klein Gordon equation, constitute a complete system of equations for R(r), A(r), and B(r): In the original RB paper, the self-interaction potential was trivial V (A) = m 2 A 2 /2, and only the N |KG[A]|N − 1 = 0 expecation value of the KG equation is nontrivial. In that case the ansatz of eq. (V.1) is an exact solution. For the φ 4 potential we have been considering, we have Evaluating the expectation values as in [2], we find 3 One can simplify the expression by defining [18] X(y) = 2 |λ| N m R(r), y = m r, (V. 8) in which case the EKG equations simplify to Note that primes in eq. (V.9) now indicate derivatives with respect to y (rather than r).
In [18], our original investigation of axion stars, we identified two small parameters and expanded the EKG equations in both. The analogue of these small parame- 3 The third equation of (V.7) is the expectation value N |KG[A]|N − 1 = 0, but in the self-interacting theory higherorder expectation values N |KG[A]|N − (2k + 1) = 0 with k > 0 will not be satisfied, and so eq. (V.1) is not an exact solution.
However, as argued in [40], taking the leading order is a good approximation for all but the most strongly bound boson star configurations. In the context of the axion potential, we have presented a method of going beyond the RB ansatz to calculate relativistic corrections perturbatively in [40]. For other work on relativistic corrections to scalar field theory, see [39,41,42].
(V.10) The requirement δ 1 corresponds to the weak gravity limit, as the δ → 0 limit recovers the Einstein equations for the vacuum. It is valid only when |λ| 8 π m 2 /M P 2 . As an example, axions have λ = −m 2 /f 2 and so this condition is equivalent to 8 π f 2 M P 2 , a condition that is easily satisfied in nearly all applications. In particular, for QCD axions, δ ≈ 10 −14 , and even theories of axionic "fuzzy dark matter", which have much larger f , still δ 1 is satisfied easily. The other parameter ∆ is small precisely when the axion star is weakly bound, i.e. when the eigenenergy µ 0 is of the same order as the particle mass m.
In applications with weak gravity, we can expand the metric components as We also rescale the axion wavefunction X(y) using with x = y ∆ (which is to say, the wavefunction and coordinate scale with ∆ as their scaling dimension). The resulting equations of motion, to leading order in δ and ∆, take the form The constant κ ≡ δ/∆ 2 controls the effective coupling to gravity, and is finite even though δ, ∆ 1. We have integrated this set of equations previously in [18] to find the spectrum of weakly bound axion stars.
Finally, we can integrate the first two equations to eliminate a(x) and obtain a Poisson-like equation for b(x): This implies that b(x) is proportional to the Newtonian gravitational potential V g . The solution of this equation is Using this expression in the third equation in (V.13), we arrive at a self-contained expression for the rescaled wavefunction Y (x): (V.16) Because this equation is true only to leading order in ∆, we refer to this equation as the infrared, or low energy, limit of the Klein-Gordon equation for the axion.
We may manipulate the GPP equations (II.13) and (II.14) into a very similar form in a straightforward way. First, rewriting the gravitational potential V g in integral form, we find a single integro-differential equation of motion (V.17) In the nonrelativistic limit, the wavefunction has a time dependence which is approximately harmonic, such that iψ = µ ψ = −(m − µ 0 )ψ ≈ −m ∆ 2 ψ/2. Further, we rescale to the dimensionless wavefunction and dimensionless coordinate x = m ∆ r as before. Then eq. (V.17) takes the form Finally, we recognize the prefactor on the gravitational term as Thus, eq. (V.19) is exactly equivalent to eq. (V.16) with the identification Thus the GP formalism is equivalent to the RB formalism used in the infrared limit. We summarize the identifications between the two paradigms in Table III. Because these formalisms are precisely equivalent, one can work with whichever is more convenient for the application at hand. For example, in analyzing axion star decay, it is better to use the RB approach, as the transition matrix elements are more directly calculable [33]. It is also more straightforward to generalize if higher-order relativistic corrections are needed [40]. For the analysis of collapsing axion stars, we instead found it more convenient to use the GP formalism [28].

B. The Thomas-Fermi Approximation
We will use the RB formulation to critically analyze the Thomas-Fermi (TF) approximation, which is commonly used in studies of repulsively interacting boson stars. At leading order, the equations of motion are where we have chosen to work with the wavefunction X(y) rather than Y (x) so that the appearance of the small parameters δ, ∆ is manifest. In the TF approximation, one neglects the kinetic term ∇ 2 y X(y) compared to the other terms. The consistency of that assumption must then be checked after the solution for b(y) and X(y) have been found. Suppose the scale for r is R, the radius of the boson star. Then the kinetic term is If the TF approximation is valid, then the Klein-Gordon equation simplifies to which can be directly solved for X(y) to obtain Substituting this back into the Poisson equation gives The solution of this equation which is regular at the origin is where p 0 = 2 √ δ and c is a dimensionless constant which will be determined below. Finally, we can calculate the wavefunction in the TF limit by substituting eq. (V.27) into eq. (V.25): We can now observe directly that the constant p 0 = π/m R T F defines the radius R T F of the condensate in the TF limit, given by  and matches the standard result [74], previously derived by considering a polytropic equation of state. Note that the radius of the condensate in the TF limit does not depend on the particle number N ; this behavior of the cutoff is characteristic of the TF approximation. The coefficient c can be determined by normalization of the wavefunction (c.f. Table III) (V.30) which implies c = N λ/(4π 2 ). Thus, we can write down the final expression for the rescaled wavefunction Finally, we consider the breakdown of the TF approximation. To check the self-consistency of omitting the kinetic term from the equation of motion, (V.21), we must substitute the solutions of eqs. (V.28) and (V.27) back into (V. 21), to see whether the kinetic term is much smaller than the rest of the terms. Substituting those solutions into the kinetic term gives a complicated expression, which simplifies considerably if we evaluate it just at the boundaries of the physical range for r, at r = 0 and r = R T F = π / (m p 0 ). In both limits we obtain, up to numerical constants, Similarly, evaluating the interaction term at the the and of the range of r we obtain The validity of the Thomas-Fermi approximation requires that the ratio of eq. (V.32) and (V.33) be small, i.e.
This implies and comparing to eq. (III.1) we see that this condition simply implies n 1. A further constraint on the Thomas-Fermi solution is the relationship between the Schwarzschild radius, R SCH and R. Unless R > R SCH the boson star collapses to a black hole. This constraint leads to the inequality In other words, there is an upper limit for N , Note that this differs from the exact result of Colpi et al. [4] only by a small numerical factor. Of course, if N N max is not satisfied, then our calculations, based on the Newtonian approximation to gravity, would not be acceptable. Therefore we obtain upper and lower (V. 34) bounds for N in the Thomas-Fermi approximation: providing a wide range for the applicability of the approximation, provided The TF approximation has been generalized to describe rotating [76] or charged [77] boson stars, as well as boson stars comprised of N scalar fields [78]. Note that the behavior of the TF wavefunction near R T F is not physical at finite N . One could in principle calculate corrections to the TF equations of motion which take this into account, which should give rise to the standard exponentially falling wavefunction at large r. This would have consequences in certain applications where the tail behavior is important, including calculations of the classical decay rate. An analysis of this type is beyond the scope of the present paper.

VI. CONCLUSIONS
In this analysis we have analyzed a number of approximate methods for describing boson stars. We focused on the Gross-Pitäevskii+Poisson (GPP) equations, which have a broad range of applicability for weakly-bound, nonrelativistic boson stars. Using a timeindependent variational formalism, we compared various ansätze which describe gravitationally bound BECs with self-interactions. These ansätze allow the GPP system to be solved analytically, alleviating the obligation to cumbersome numerical solutions. Moreover, numerical solutions exist primarily for stationary BEC configurations, and are much more difficult to use in dynamical applications; ansätze are powerful tools for solving dynamic BEC problems such as collapse, collisions, and expansion.
We have treated the numerical solution to the stationary GPP system as a benchmark for comparing different ansätze, prioritizing factors such as computational ease, fit of the wavefunction profile, and value of maximum mass (for attractive interactions). We found that a linear+exponential wavefunction is the best fit for attractive self-interactions along the stable branch, as well as for repulsive self-interactions at small N . For attractive self-interactions along the unstable branch a single exponential is the best fit for small N while a sech wavefunction fits better for large N , though the latter is computationally inefficient. A Gaussian wavefunction, which is used often in the literature, is exceedingly poor across most of the parameter space, with the exception of repulsive interactions for large N .
We found that our proposed double-exponential ansatz is much more tunable compared to other ansätze in the literature. By choosing various values of the parameter a, one may optimize a given fit parameter over others, or near-perfectly replicate more computationally complex ansätze. In particular, on the stable branch for attractive interactions, the limit a → 1 gives rise to the lin-ear+exponential ansatz, which we found to be in closest agreement to the exact case. On the unstable branch, one can take the a 1 limit to obtain the exponential ansatz, appropriate at large central densities, while the choice a = 2 very nearly reproduces the sech ansatz, which is quite inefficient computationally. We also showed how to generalize the addition of free parameters in order to create more computationally efficient ansätze using the double exponential.
There remain many unstudied applications of the ansätze we have compiled here, and analytic solutions to the time-dependent GPP system will undoubtedly have myriad uses. As mentioned in the introduction, there are many open questions regarding oscillons formed by inflation and quintessence fields, and axion BECs remain the subject of active study. Along with a rigorous, multicriterion comparison of ansätze from the relevant literature, we have presented arguments for the necessity of time-dependent ansätze, and have introduced a tunable, computationally simple ansatz.
Finally, we also highlighted the relevant differences between the relativistic Einstein+Klein-Gordon (EKG) equations, as derived using the Ruffini-Bonazzola formalism, and the nonrelativistic GPP system. In particular, we showed explicitly the relevant expansion parameters which characterize the nonrelativistic limit. We used this formulation to critically analyze the Thomas-Fermi approximation, which is a large-N limit for repulsive selfinteractions in which the kinetic energy is taken to be small.
An alternative formulation, used for example in [14,70,71], is as follows: The energy per particle is where by comparison with (III.5), we identify (A. 2) The energy per particle of (A.1) is minimized at In the special case of λ < 0 (attractive selfinteractions), we find the critical number and radius beyond which no bound state solutions exist, We can also write the bound state solution in terms of N c and σ c ,  where δ ≡ 1 − N 2 /N c 2 . In these variables the energy per particle can then be written as (A.6) The corresponding value for the chemical potential is (A.7) Evaluated at the critical values (i.e. N = N c and δ = 0), the energy per particle and chemical potential are, which will be true independently of the choice of ansatz.

Appendix B: Full Solutions for Double Exponential Ansatz
Tabulated in Table IV are the full solutions for the parameters in Table I for the double exponential ansatz for arbitrary values of a. Note thatn and ρ c can be calculated from these coefficients by using eqs. (III.10) and (III.13).