Self-gravitating Bose-Einstein condensates and the Thomas-Fermi approximation

: Self-gravitating Bose-Einstein condensates (BEC) have been proposed in various astrophysical contexts, including Bose-stars and BEC dark matter halos. These systems are described by a combination of the Gross-Pitaevskii and Poisson equations (the GPP system). In the analysis of these hypothetical objects, the Thomas-Fermi (TF) approximation is widely used. This approximation is based on the assumption that in the presence of a large number of particles, the kinetic term in the Gross-Pitaevskii energy functional can be neglected, yet it is well known that this assumption is violated near the condensate surface. We also show that the total energy of the self-gravitating condensate in the TF-approximation is positive. The stability of a self-gravitating system is dependent on the total energy being negative. Therefore, the TF-approximation is ill suited to formulate initial conditions in numerical simulations. As an alternative, we offer an approximate solution of the full GPP system.

Nonrelativistic, dilute BECs are well-described by the Gross-Pitaevskii equation (GPE). In its time-independent form, when kinetic energy can be neglected, the GPE has an approximate solution in the form of the Thomas-Fermi (TF) approximation.
Self-gravitating systems can also be described using the GPE. The long-range gravitational interaction is represented by the external potential. Such systems, referred to in the literature as GPE-Poisson or GPN (Gross-Pitaevskii-Newton) systems, are extensively studied (a few recent examples of interest include [12][13][14][15]).
When applied to self-gravitating systems, the TF-approximation was found to describe an unstable system [16]. This also agrees with this author's experience with numerical simulation code, which failed to yield stable self-gravitating systems when the TF-approximation was used as an initial estimate of the condensate density.
In the present paper, it is shown that this instability is the consequence of the TF-approximation, specifically the well-known issue of the divergence of its kinetic energy [7,8,17,18]. To overcome this issue, we examine the time-independent GPE and attempt to solve it numerically without truncations. We develop a new approximation that has the desirable property that the system has negative total energy and it is stable. This approximation has since been incorporated into our simulation code for self-gravitating GPE-Poisson systems.

Discussion
A self-gravitating Bose-Einstein condensate is described by a combination of the Gross-Pitaevskii equation and Poisson's equation for gravity. In units such that the BEC particle mass is m = 1 and alsoh = 1, the time-independent Gross-Pitaevskii equation can be written in an attractively simple form: where Ψ is the BEC wavefunction, V is the gravitational potential, c is the BEC coupling coefficient and µ is the chemical potential, the presence of which guarantees the conservation of energy. We normalize the wavefunction such that the number of particles is V |Ψ| 2 = N. The energy functional from which with the GPE (1) can be derived using the variational principle (cf. [18,19]; note the additional factor of 1/2 in front of V, required to avoid double counting the gravitational potential energy between two regions of the condensate as V is itself a function of |Ψ| 2 ) is given by In the Thomas-Fermi (TF) approximation [6][7][8], kinetic energy is neglected. Therefore, the time-independent GPE takes on the following simplified form: If V is not dependent on Ψ, this equation can be solved directly for |Ψ| 2 : Moreover, if we require the wavefunction to vanish at infinity, we must have at infinity. If V is dependent on Ψ, the situation becomes somewhat more complicated. In particular, in the GPP system, the relationship between V and Ψ is given by Poisson's equation: where G is the gravitational constant. Solving the GPE (1) for V in the TF-approximation, and substituting this solution back into Poisson's equation (6), we get If µ = const., we are left with which is an homogeneous Helmholtz-type equation for |Ψ| 2 , spherically symmetric solutions of which are where k 2 = 4πG/c, while C 1 and C 2 are integration constants. To avoid solutions that are singular at the origin r = 0, we must set C 2 = 0. On the other hand, sin kr/r (and thus, |Ψ| 2 ) vanishes at r = π/k. Therefore, we set r 0 = π/k as the radius of the condensate. This determines C 1 since we require that This integral can be readily evaluated: hence C 1 = N/4r 2 0 . Therefore, the TF-approximation for the GPP is given by the Lane-Emden type solution for 0 ≤ r ≤ r 0 = √ πc/4G. Given |Ψ| 2 , we can solve the GPE (1) for V: again for 0 ≤ r ≤ r 0 . At r 0 and beyond, the condensate vanishes, and the gravitational potential becomes that of a point mass M (where M = Nm is the total mass of the condensate), i.e., V = −GN/r (r 0 ≤ r). At the boundary, these two forms must agree. This can be achieved by setting This clarifies the role of the chemical potential in the case of the GPP system in the TF-approximation: its presence ensures that the gravitational potential takes on the standard form outside the condensate and vanishes at infinity. The energy density of the time-independent GPE in the Thomas-Fermi limit is given by or, after substituting the solution for V from the GPE (1), To find the total energy, we integrate over the condensate volume: or, after restoring units, The positive sign of the total energy implies that the solution for a self-gravitating BEC using the TF-approximation is inherently unstable. This result is based on the assumption that the kinetic energy can be neglected. Now that we have an explicit solution for |Ψ| 2 , this assumption can be verified by direct substitution into the energy functional (2). When we do so we find that, using the solution given by Eq. (13), the condensate kinetic energy, is divergent for any r 0 > 0. The implication of this divergence is that the assumption behind the TF-approximation, namely that the kinetic term in the GPE (1) can be neglected, is maximally violated.
The divergence of the kinetic energy density (represented by the integrand in Eq. (21)) in the Thomas-Fermi approximation, near the condensate surface, and the logarithmic divergence of the total kinetic energy are well known 1 [7,8,17,18]. However, in the case of self-gravitating systems, as we have seen, the problem gets worse: the total energy of a self-gravitating condensate is positive, hence the condensate is gravitationally unstable. This has especially important implications for numerical simulations of self-gravitating Bose-Einstein condensates that use this approximation to model the initial state (see, e.g., [16,20]). Not only is the magnitude of the total kinetic energy ill-defined (which makes it dependent on nonphysical simulation parameters, such as the numerical integration step size or even small rounding errors), the positive total energy is especially troublesome, as the stability of a self-gravitating system is dependent on E < 0.
This finding agrees with the author's experience using numerical simulation code [20] that was designed to model the (time-dependent) GPP system, using the TF-approximation to model the initial condensate density. It was perplexing that versions of the code ran differently in different programming environments (e.g., single vs. double precision, FORTRAN vs. C), processors (Intel x86 vs. GPGPU) and operating systems (Linux vs. Windows). This is clearly not permissible: the results, apart from accuracy and rounding issues, should not be dependent on such factors. Indeed, the present study arose as a result of systematically analyzing the failure of these algorithms to produce consistent results.
Numerically stable simulations require an initial state that is not dependent on nonphysical parameters and has a well-defined total energy that is consistent with stability conditions. This is clearly not the case when the TF-approximation is used for a self-gravitating condensate. Therefore, 1 The author wishes to thank the anonymous referees and the Academic Editor of Galaxies for stressing this point.
we now aim to find an approximate solution of the GPP system in the spherically symmetric case without resorting to the TF-approximation (see also [7] for another approximation and [8] for a numeric solution in terms of hydrodynamics using the Madelung-representation). Assuming spherical coordinates and a spherically symmetric condensate, the GPE (1) becomes an equation of a single independent variable r. Let us denote ∂/∂r = ∂ r for brevity, while noting the form of the Laplacian in spherical coordinates, ∇ 2 = ∂ 2 r + (2/r)∂ r . We can then write the GPE (1) as whereas Poisson's equation (6) becomes: Let us now consider writing the wavefunction as in which case and the GPE (1) can be rewritten, after dividing through by e iφ , as Since |Ψ|, φ, V, c and µ are all real, the real and imaginary parts of this equation can be separated: Equation (29) can be integrated: where C is an integration constant with the dimensions of r|Ψ| 2 . It seems that C = 0 is not only a valid choice but the only choice that does not result in φ becoming singular at the origin (assuming the wavefunction does not vanish at the origin). Therefore, φ = const. Under these circumstances, the GPE (1) will read which is readily solvable for V algebraically: This result can be substituted back into Poisson's equation, yielding a fourth-order ordinary differential equation in |Ψ|: This equation can be solved numerically. Given that it is a fourth-order homogeneous differential equation in |Ψ|, it has a very large solution space, parameterized by boundary or initial conditions, such as the values of |Ψ| and its first three derivatives at some value of r. A hint for a suitable solution comes from numerical simulation [20], where we find that apparently stable nonrotating spherically symmetric solutions converge on |Ψ| 2 ∝ [sin(r/r 0 )/(r/r0)] 2α , with 3 α 4. This approximate solution has many desirable properties. It is smooth in the interval 0 ≥ r/r 0 ≥ π. The corresponding kinetic energy (21) is finite in the same interval. Moreover, it is possible to compute the condensate mass, which is given by and this, too, is finite and well-behaved. Therefore, we find that the following initial approximation for the magnitude of the BEC-Poisson wavefunction: agrees well with a numerical solution (except for very small values of r), as shown in Fig. 1. Furthermore, this choice yields a corresponding solution of Eq. (6) for the gravitational potential that is finite, negative, and vanishes at infinity, as expected. The stability of these solutions is confirmed by numerical simulation of a condensate using Eq. (35) as an initial approximation. An example result is shown in Fig. 2.

Conclusions
We have demonstrated that when the Thomas-Fermi approximation is used to describe self-gravitating Bose-Einstein condensates in astrophysical contexts, the resulting systems have divergent (positive) total energy and are unstable. However, this behavior is a specific consequence that arises from the use of the TF-approximation; GPE-Poisson systems are not inherently unstable. By investigating the untruncated Gross-Pitaevskii equation, we found that the total energy is, in fact negative. Furthermore, using an approximate numerical solution as a guide, we developed a simple approximation formula that can be used to provide an initial density estimate for BECs in numerical . Density cross-section of a stable simulated 1 M ⊙ , r ≃ 50 km Bose-star (or stellar core) after approximately 300,000 numerical iterations that corresponds to 3 seconds [21]. (For comparison, the period of a circular orbit at r = 50 km is approximately 0.006 s.) Axes are in km, density is in units of 10 30 kg/km 3 ≃ 0.5 M ⊙ /km 3 . NB: The spatial grid used in this simulation is low resolution (60 × 60 × 60) and image smoothing was used to improve the presentation quality.
simulations. These simulations can be applied to describe, e.g., BEC dark matter stars or stellar cores; these results and analysis will be reported elsewhere [21] as they become available.