The formation of compact objects at finite temperatures in a dark-matter-candidate self-gravitating bosonic system

We study self-gravitating bosonic systems, candidates for dark-matter halos, by carrying out a suite of direct numerical simulations (DNSs) designed to investigate the formation of finitetemperature, compact objects in the three-dimensional (3D) Fourier-truncated Gross-Pitaevskii-Poisson equation (GPPE). This truncation allows us to explore the collapse and fluctuations of compact objects, which form at both zero temperature and finite temperature. We show that the statistically steady state of the GPPE, in the large-time limit and for the system sizes we study, can also be obtained efficiently by tuning the temperature in an auxiliary stochastic Ginzburg-Landau-Poisson equation (SGLPE). We show that, over a wide range of model parameters, this system undergoes a thermally driven first-order transition from a collapsed, compact, Bose-Einstein condensate (BEC) to a tenuous Bose gas without condensation. By a suitable choice of initial conditions in the GPPE, we also obtain a binary condensate that comprises a pair of collapsed objects rotating around their center of mass.

Gravitational effects are important on stellar scales; it might also be possible to mimic such effects in laboratory Bose-Einstein condensates [1] and thus emulate gravitationally bound, condensed, assemblies of bosons, which are candidates for dark-matter halos [2][3][4][5]. Although many experiments have been carried out to establish the identity of dark matter, there is still no unambiguous dark-matter candidate. For many years the front runners have been weakly interacting massive particles (WIMPs); but their existence has not been established convincingly (see, e.g., [6,7] and references therein), so investigations of other candidates, e.g., axions, boson stars, black holes, and superfluids, have experienced a renaissance [8][9][10][11][12][13][14].
The Gross-Pitaevskii-Poisson equation (GPPE), for a self-gravitating assembly of weakly interacting bosons, is the natural theoretical model for dark-matter-candidate bosonic assemblies, both in laboratory and astrophysical settings. We address finite-temperature (T > 0) effects in such condensation, a problem of central importance in this challenging field. Several ways have been suggested for including T > 0 effects in the Gross-Pitaevskii (GP) model [15] without gravity; one important way uses the Fourier-truncated GP model, in which this truncation generates a classical-field model [15][16][17][18]. We generalize these studies by using the Fourier-truncated GPPE to study T > 0 effects during the gravitational collapse of a system of self-gravitating bosons. In addition, we define an algorithm that directly reconstructs the thermalized state of such a system by using an auxiliary stochastic Ginzburg-Landau-Poisson equation (SGLPE).
We obtain several interesting results by pseudospectral direct numerical simulations (DNSs) of the truncated GPPE and the SGLPE in three dimensions (3D). We follow the spatiotemporal evolution of different initial conditions for the density of bosons. If we start from a very-nearly uniform density in the truncated GPPE, this system undergoes gravitational collapse and thermalizes to a Bose-Einstein condensate (BEC), at low T . Our SGLPE study yields a hitherto unanticipated, thermally driven, first-order phase transition from a collapsed, compact, BEC to a tenuous Bose gas without condensation, for a wide variety of parameters in the GPPE. Finally, by a suitable choice of initial conditions in the GPPE, we obtain a binary condensate that comprises a pair of collapsed objects rotating around their center of mass.
We solve the GPPE (1) by using the 3D Fourier pseudospectral method, with the 2/3-rule for dealiasing [17,18,23]: We expand the 2π periodic wave function as ψ(x) = k∈Z 3ψk exp(ik · x) and then we truncate it spectrally, by settingψ k ≡ 0 for |k| > k max , with where N is the resolution and [·] denotes the integer part [15,16]. If we introduce the Galerkin projector P G [in Fourier space P G [ψ k ] = θ(k max −|k|)ψ k , with θ(·) the Heaviside function], the Fourier-truncated GPPE becomes Equation (2) conserves exactly the number of particles If we use the 2/3-rule for dealiasing, then the momentum where the over-bar denotes complex conjugation, is also conserved [17].
This spectral truncation generates a classical-field model that allows us to study finite-T effects in the GPPE (Refs. [16][17][18] for the GP case). We show that this spectrally truncated GPPE can describe dynamical effects and, at the same time, yield thermalized states, which we can obtain both by the thermalization of the long-time GPPE dynamics or directly, by using the SGLPE where the zero-mean, Gaussian white noise ζ(x, t) has the variance ζ(x, t)ζ * (x , t ) = δ(t − t )δ(x − x ), with β = 1 k B T the inverse temperature; we tune the chemical potential µ at each time step to conserve the total number of particles N . The finite-T SGLPE dynamics does not describe any physical evolution; but it converges more rapidly, than does the GPPE dynamics, toward a thermalized state (Refs. [16][17][18] for the GP case). Clearly, the SGLPE leads to a state with a given temperature; but the GPPE yields a state with a given energy.
videos V1-V5 (Supplemental Material [19]) show, respectively, the complete spatiotemporal evolution of |ψ(x, t)| 2 for these five runs. In column (4) of Fig. 1 we give, for these runs, plots of the scaled radius of gyration V ρ(r)dr (blue curves) and the scaled gravitational energy E grav /E a (red curves) versus the scaled time t/(ξ/v), where E a = 2 5 π 4 (G/g 3 ) 1/2 . If we tune T in the SGLPE (3), it yields a statistically steady state whose properties (like R/L and E grav /E a ) are close to their counterparts in the thermalized state of the GPPE (e.g., by comparing rows (2) and (3) in column(4) of Fig. 1 we see that the final state of the SGLPE has nearly the same energy and radius as the corresponding statistically steady GPPE state). Furthermore, convergence to this thermalized state is more rapid in the (canonical) SGLPE than in the (microcanonical) GPPE. To validate our DNSs, we have checked explicitly (Supplemental Material [19]) that, at zero temperature, our results agree with those of the T = 0 study of Refs. [10,11], which yield spherically symmetric ground states ( row (1) of Fig. 1) with radius R and N bosons. Their ground-state energy can be approximated as E(R) = 2 N 2mR 2 + gN 2 2R 3 − GN 2 2R and the equilibrium radius follows from dE/dR| R=R0 = 0, GNm 3 . The details of our extensive DNSs for the GPPE and the SGLPE are given in the Supplemental Material [19] (see, especially, Table I there). Figure 2(a) shows plots of of R/R Q versus a/a Q for three different values of G. In Fig. 2(b) we give plots of M/M a versus R/R a , from our GPPE DNSs for three different values of g; the trends in these plots are markedly different than those at T = 0 (see Fig.(1) in the Supplemental Material [19]). In Fig. 2(c) we plot the scaled energies E kq /E a , E int /E a , E G /E a , and their total E/E a versus the scaled temperature k B T /E a from our SGLPE DNSs; the jumps in these curves at k B T /E a 3.25 × 10 −5 suggest a first-order transition from a collapsed BEC to a tenuous, non-condensed assembly.
To confirm such a transition, we carry out an SGLPE hysteresis study, whose results are summarised in the plots in the panels of Fig. 3. The left panel of this figure shows R/L versus k B T /E a ; we obtain the red and green curves by, respectively, increasing and decreasing k B T /E a (often referred to as heating and cooling runs in statistical mechanics). In these SGLPE runs, we use the final steady-state configuration for ψ(x), from the previous temperature, as the initial condition at the next temperature; clearly, there is significant hysteresis at the first-order transition from the -collapsed BEC to the non-collapsed state. In the right side panels of Fig. 3, we show ten-level contour plots of |ψ(x)| 2 and the associated spectra |ψ(k)| 2 to illustrate, at representative points on heating and cooling curves in the hysteresis plot, the real-space density distribution and the k−space density spectra (k B T /E a = 2.7 × 10 −5 and k B T /E a = 3.62 × 10 −5 in the top panels A and B, respectively, and k B T /E a = 2.3 × 10 −5 and k B T /E a = 3.16 × 10 −5 in the bottom panels C and D, respectively). From these density distributions and spectra, we conclude that our system undergoes a first-order transition from a collapsed BEC to a tenuous, non-condensed assembly. The panel of figures at the very bottom of Fig. 3 show that this first-order transition occurs at g = 0 too. We expect that this also occurs when g < 0, which is the appropriate parameter range for axion stars [13,14]; the g < 0 case requires a quintic nonlinearity in Eq. 1 for stability; finite-temperature effects can be studied for this case by using the methods that we have described above (as we will show in future work).
Do our Fourier-truncated GPPE yield only single collapsed objects? No. We now show, for the illustrative parameter values g = 20, G = 1000, and 128 3 collocation points, that this GPPE can also yield long-lived states with temporal oscillations. In particular, by using an initial condition with two rotating spherical compact objects, the truncated GPPE dynamics yields a rotating binary system, which we depict at representative times in Fig. 4 and in the Video V6 in the Supplementary Material [19] (this also gives the initial data).
Our study goes well beyond earlier studies [24,25] that employ the Thomas-Fermi approximation and include finite-temperature effects at the level of a noninteracting Bose gas (with g = 0). We have shown that the truncated GPPE can be used effectively to study the graviational collapse of an assembly of weakly interacting bosons at finite temperature. Our study of this collapse shows that there is a clear, thermally driven first-order transition from a non-condensed bosonic gas, at high T , to a condensed BEC phase at low T ; this transition should be contrasted with the continuous BEC transition in the weakly interacting Bose gas, which is described by the GPE, in the absence of gravitation. Furthermore, we have shown that gravitationally bound, rotating binary objects can be obtained in our GPPE simulations. Therefore, our work opens up the possibility of carrying out detailed finite-temperature studies of self-gravitating bosonic systems, which are potentially relevant for studies of dark-matter candidates like boson stars and axions. FIG. 2: The plots of (a) the scaled radius of gyration (R/RQ) versus scaled scattering length (a/aQ) for the GPPE Runs A1-A30, (b) the scaled mass (M/Ma) versus the scaled radius of gyration R/Ra for GPPE Runs B1-B30 on a log-log scale, (c) the scaled energy components (E kq /Ea, Eint/Ea, EG/Ea), and the total energy (E/Ea) versus the scaled temperature kBT /Ea.

FIG. 3:
Top left panel: plots of the dimensionless radius (R/L) versus the dimensionless temperature (kBT /Ea), for heating (red) and cooling (green) runs showing a hysteresis loop. We show ten-level contour plots of |ψ(x)| 2 and the associated spectra |ψ(k)| 2 to illustrate, at representative points on heating and cooling curves in the hysteresis plot, the real-space density distribution and the k−space density spectra (kBT /Ea = 2.7 × 10 −5 and kBT /Ea = 3.62 × 10 −5 in panels A and B of top panels, respectively, and kBT /Ea = 2.3 × 10 −5 and kBT /Ea = 3.16 × 10 −5 in panels C and D of bottom panels, respectively). The analogs of these plots, for the case g = 0, are shown in the panels at the very bottom. In the bottom left panel we use |EG| at T = 0 to make the temperature dimensionless.