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

Akhilesh Kumar Verma<sup>1,*,†</sup> Rahul Pandit,<sup>1,†</sup> and Marc E. Brachet<sup>2,‡</sup>

¹Centre for Condensed Matter Theory, Department of Physics, Indian Institute of Science, Bangalore 560012, India

²Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université Université de Paris, F-75005 Paris, France

<figure><img src="image_1.png" /></figure>

(Received 19 December 2019; accepted 7 April 2021; published 27 May 2021)

We study self-gravitating bosonic systems, candidates for dark-matter halos, by carrying out a suite of direct numerical simulations designed to investigate the formation of finite-temperature, 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. 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 to a tenuous Bose gas (that is not gravitationally condensed). 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.

DOI: 10.1103/PhysRevResearch.3.L022016

Gravitational effects are important on stellar scales; it might also be possible to mimic such effects in laboratory Bose-Einstein condensates (BEC) [1] and thus emulate gravitationally bound, condensed assemblies of bosons, which are candidates for dark-matter halos [2–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., Refs. [6,7] and references therein), so investigations of other candidates, e.g., axions, boson stars, black holes, and superfluids, have experienced a renaissance [8–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 in 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–18] (see also Supplemental Material [22]). 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 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.

A self-gravitating BEC is described by a complex wave function $\psi(\mathbf{x}, t)$; for weakly interacting bosons, the spatiotemporal evolution of $\psi(\mathbf{x}, t)$ is governed by the GPPE:

$$i\hbar\partial_t\psi = -\frac{\hbar^2}{2m}\nabla^2\psi + [G\Phi + g|\psi|^2]\psi, \\ \nabla^2\Phi = |\psi|^2 - (|\psi|^2), \quad (1)$$

where $m$ is the mass of the bosons, $n = |\psi|^2$ their number density, $\Phi$ is the gravitational potential field, $G = 4\pi G_N m^2$ ($G_N$ is Newton’s constant), and $g = 4\pi a\hbar^2/m$, with $a$ the

<footnote>*akhilesh@iisc.ac.in</footnote>

<footnote>†Also at Jawaharlal Nehru Centre For Advanced Scientific Research, Jakkur, Bangalore, India; rahul@iisc.ac.in</footnote>

<footnote>‡brachet@physique.ens.fr</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI.</footnote>
s-wave scattering length. The subtraction of the mean density $\langle |\psi|^2 \rangle$ can be justified either by taking into account the cosmological expansion [19,20] or by defining a Newtonian cosmological constant [21]. By linearizing Eq. (1) around the constant $|\psi|^2 = n_0$, we obtain the dispersion relation $\omega(k) = \sqrt{-Gn_0/m + k^2gn_0/m + k^4(\hbar/2m)^2}$, which has been studied in detail in Ref. [10]. This dispersion relation displays a low-k Jeans instability for wave numbers $k < k_j = \sqrt{G/g}[(1 + \sqrt{1 + \frac{G\hbar^2}{mg^2n_0}})/2]^{-1/2}$. In the absence of gravity ($G=0$), we identify the speed of sound $c = \sqrt{gn_0/m}$ and the coherence length $\xi = \sqrt{\frac{\hbar^2}{2gn_0m}}$. Units relevant to astrophysics are discussed in Ref. [14].

We solve the GPPE (1) by using the 3D Fourier pseudospectral method, with the 2/3-rule for dealingiasing [17,18,23]: We expand the $2\pi$ periodic wave function as $\psi(x) = \sum_{\mathbf{k} \in \mathbb{Z}^3} \hat{\psi}_{\mathbf{k}} \exp(i\mathbf{k} \cdot \mathbf{x})$ and then we truncate it spectrally by setting $\hat{\psi}_{\mathbf{k}} \equiv 0$ for $|\mathbf{k}| > k_{\max}$, with $k_{\max} = [N/3]$, where $N$ is the resolution and $[\cdot]$ denotes the integer part [15,16]. If we introduce the Galerkin projector $\mathcal{P}_G$ {in Fourier space $\mathcal{P}_G[\hat{\psi}_{\mathbf{k}}] = \theta(k_{\max} - |\mathbf{k}|)\hat{\psi}_{\mathbf{k}}$, with $\theta(\cdot)$ the Heaviside function}, the Fourier-truncated GPPE becomes

$$i\hbar \frac{\partial \psi}{\partial t} = \mathcal{P}_G \left\{ -\frac{\hbar^2}{2m} \nabla^2 \psi + \mathcal{P}_G \left[ (G\nabla^{-2} + g)|\psi|^2 \right] \psi \right\}. \quad (2)$$

Equation (2) conserves exactly the number of particles $N = \int d^3x |\psi|^2$ and the energy $E = E_{kq} + E_{int} + E_G$, where $E_{kq} = \frac{\hbar^2}{2m} \int d^3x |\nabla\psi|^2$, $E_{int} = \frac{g}{2} \int d^3x [\mathcal{P}_G |\psi|^2]^2$, and $E_G = E_{grav} = \frac{G}{2} \int d^3x [\mathcal{P}_G |\psi|^2] \nabla^{-2} [\mathcal{P}_G |\psi|^2]$. If we use the 2/3-rule for dealingiasing, then the momentum $\mathbf{P} = \frac{i\hbar}{2m} \int d^3x (\psi \nabla \bar{\psi} - \bar{\psi} \nabla \psi)$, where the overbar 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–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 by the thermalization of the long-time GPPE dynamics or directly by using the SGLPE,

$$\hbar \frac{\partial \psi}{\partial t} = \mathcal{P}_G \left\{ \frac{\hbar^2}{2m} \nabla^2 \psi + \mu \psi - \mathcal{P}_G \left[ (G\nabla^{-2} + g)|\psi|^2 \right] \psi \right\} + \sqrt{\frac{2\hbar}{\beta}} \mathcal{P}_G[\zeta(\mathbf{x}, t)], \quad (3)$$

where the zero-mean, Gaussian white noise $\zeta(\mathbf{x}, t)$ has the variance $\langle \zeta(\mathbf{x}, t) \zeta^*(\mathbf{x}', t') \rangle = \delta(t-t')\delta(\mathbf{x}-\mathbf{x}')$, with $\beta = \frac{1}{k_B T}$ the inverse temperature; we tune the chemical potential $\mu$ 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–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.

Our pseudospectral DNSs of the GPPE (1) and SGLPE (3) use a cubical computational domain that is $(2\pi)^3$ periodic; we normalize $\psi$ such that $N = (2\pi)^3$ and use units with $\hbar = 1$

<table><caption>TABLE I. The representative runs for which we give plots in Fig. 1.</caption><thead><tr><th>Run Type</th><th>R1 SGLPE (T = 0)</th><th>R2 SGLPE</th><th>R3 GPPE</th><th>R4 GPPE</th><th>R5 GPPE</th></tr></thead><tbody><tr><td>N</td><td>64</td><td>64</td><td>64</td><td>64</td><td>256</td></tr><tr><td>g</td><td>50</td><td>50</td><td>50</td><td>5</td><td>50</td></tr><tr><td>G</td><td>105</td><td>105</td><td>105</td><td>550</td><td>105</td></tr></tbody></table>

and $m=1$. We list below dimensionless variables that we use; here, $G_N$ denotes Newton's constant and $G = 4\pi G_N m^2$:

$$M_a = \frac{\hbar}{\sqrt{G_N ma}}; \quad \frac{M}{M_a} = \frac{M \sqrt{G_g}}{4\pi \hbar^2}; \\ R_a = \sqrt{\frac{a\hbar^2}{m^3 G_N}}; \quad \frac{R}{R_a} = \frac{R}{\sqrt{g/G}}; \\ a_Q = \frac{\hbar^2}{G_N M^2 m}; \quad \frac{a}{a_Q} = \frac{M^2 G_g}{(4\pi \hbar)^2}; \\ R_Q = \frac{\hbar^2}{G_N m^2 M}; \quad \frac{R}{R_Q} = \frac{M G R}{4\pi \hbar^2}.$$

We list the parameters for different runs in Table I. We use the following random-phase initial condition in Fourier space: $|\hat{\psi}_{\mathbf{k}}| = 1$ for $k = 0$, $|\hat{\psi}_{\mathbf{k}}| = 10^{-2}/k$ for $1 \le k \le 10$ and $|\hat{\psi}_{\mathbf{k}}| = 0$ otherwise. Note that more evolved initial conditions have been studied in Refs. [24–27] in a cosmological context (with Hubble expansion and $g=0$). These cosmological initial condition (power spectra with a cutoff at small scales) are defined in Section 3 of Ref. [27]. Our simple random-phase initial conditions are sufficient for triggering the Jeans instability and then generating gravitational collapse.

In columns 1–3 of Fig. 1, we show 10-level contour plots of $|\psi(\mathbf{x}, t)|^2$ to illustrate, at representative times, the spatial organization of $|\psi(\mathbf{x}, t)|^2$ that we obtain via the SGLPE ($T=0$) (top row, run R1), the SGLPE (second row, run R2), and three GPPE runs (third row, run R3; fourth row, run R4; fifth row, run R5); Supplemental Videos S1–S5 (Supplemental Material [22]) show, respectively, the complete spatiotemporal evolution of $|\psi(\mathbf{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 $R/L = \frac{1}{L} \sqrt{\frac{f_V \rho(r) r^2 dr}{f_V \rho(r) d\tau}}$ (blue curves) and the scaled gravitational energy $E_{grav}/E_a$ (red curves) versus the scaled time $t/(\xi/v)$, where $E_a = 2^5 \pi^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 [22]) 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
<figure><img src="image_2.png" /><figcaption>FIG. 1. Columns 1–3 show 10-level contour plots of $|\psi(\mathbf{x}, t)|^2$ at representative times: SGLPE ($T=0$) (top row, run $R1$), SGLPE (second row, run $R2$), GPPE (third row, run $R3$), GPPE (fourth row, run $R4$), and $256^3$ GPPE (fifth row, run $R5$; Supplemental Videos S1–S5 in the Supplemental Material [22] show, respectively, the complete spatiotemporal evolution for these cases). Column 4: Plots of the scaled radius of gyration $R/L = \frac{1}{L}\sqrt{\frac{\int_V \rho(r)^2 dr}{\int_V \rho(r)d\mathbf{r}}}$ (blue) and the scaled gravitational energy $E_{grav}/E_a$ (red) versus the scaled time $t/(\xi/v)$ for the different runs, where $E_a = 2^5\pi^4(G/g^3)^{1/2}$.</figcaption></figure>

approximated as $E(R) = \frac{\hbar^2 N}{2mR^2} + \frac{gN^2}{2R^3} - \frac{GN^2}{2R}$ and the equilibrium radius follows from $dE/dR|_{R=R_0} = 0$, whence $R_0 = \frac{R_Q}{4\pi}(1 + \sqrt{1 + 48\pi^2(\frac{R_a}{R_Q})^2})$, where $R_Q = \frac{\hbar^2}{G_N m^3 N}$ and $R_a = \sqrt{\frac{a\hbar^2}{G_N m^3}}$. The details of our extensive DNSs for the GPPE and the SGLPE are given in the Supplemental Material [22] (see, especially, Table I there). Figure 2(a) shows plots 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 [22]). 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 \simeq 3.25 \times 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 summarized 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,
<figure><img src="image_3.png" /><figcaption>FIG. 2. The plots of (a) the scaled radius of gyration ($R/R_Q$) versus scaled scattering length ($a/a_Q$) for the GPPE Runs A1-A30, (b) the scaled mass ($M/M_a$) versus the scaled radius of gyration $R/R_a$ for GPPE Runs B1-B30 on a log-log scale (the dashed line indicates a scaling exponent of 3/4), (c) the scaled energy components ($E_{kq}/E_a$, $E_{int}/E_a$, $E_G/E_a$), and the total energy ($E/E_a$) versus the scaled temperature $k_B T/E_a$.</figcaption></figure>

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 $\psi(\mathbf{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 noncollapsed state. In the right-side panels of Fig. 3, we show 10-level contour plots of $|\psi(\mathbf{x})|^2$ and the associated spectra $|\psi(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 \times 10^{-5}$ and $k_B T/E_a = 3.62 \times 10^{-5}$ in the top panels (a) and (b), respectively, and $k_B T/E_a = 2.3 \times 10^{-5}$ and $k_B T/E_a = 3.16 \times 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, noncondensed assembly. The panel of figures at the very bottom of Fig. 3 shows 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).

Does 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

<figure><img src="image_4.png" /><figcaption>FIG. 3. Top left panel: plots of the dimensionless radius ($R/L$) versus the dimensionless temperature ($k_B T/E_a$), for heating (red) and cooling (green) runs showing a hysteresis loop. We show 10-level contour plots of $|\psi(\mathbf{x})|^2$ and the associated spectra $|\psi(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 \times 10^{-5}$ and $k_B T/E_a = 3.62 \times 10^{-5}$ in (a) and (b) of the top panels, respectively, and $k_B T/E_a = 2.3 \times 10^{-5}$ and $k_B T/E_a = 3.16 \times 10^{-5}$ in (c) and (d) of the 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 $|E_G|$ at $T=0$ to make the temperature dimensionless.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 4. Ten-level contour plots of the $|\psi(\mathbf{x}, t)|^2$ (a) at $t = 0.018$, (b) at $t = 0.025$, and (c) at $t = 0.03$, for the initial condition for $\psi(\mathbf{x}, t)$ given in Eq. (1) of the Supplemental Material [22], showing the rotating binary system (see Supplemental Video S6 in the Supplemental Material).</figcaption></figure>

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 Supplemental Video S6 in the Supplemental Material [22] (this also gives the initial data).

In astrophysics, temperature effects are expected to be extremely small as cosmological dark matter, away from collapsed structures, is very cold. However, typical dark matter, in today's universe, is found principally as collapsed-dark-matter haloes, whose temperature is determined just by the properties of the halo [28]. The maximal mass of stable BEC stars is insensitive to the temperature of the matter. However, the mass of the boson star, with a relatively large radius, can depend significantly on the temperature of the boson matter [29].

We have presented a first-principles computational study of the possibility of taking into account thermal effects in superfluid gravitational systems by using spectral truncation. The results have been given in terms of dimensionless scaled temperature

Should a concrete astrophysical situation be studied using a quantitative spectrally truncated model, the physical temperature scale could be estimated as done for $G = 0$ in Ref. [30], which shows that the classical truncated dynamics is a good approximation to the physical dynamics only for phonon modes with a high occupation number; these modes are in equipartition at that given temperature. The basic idea to match the physical temperature is to count these modes and adjust the truncation wave number accordingly. We must generalize the study of Ref. [30] to the $G > 0$ case to adjust

the physical temperature scale in a quantitative study. This will be addressed in future work.

Our study goes well beyond earlier studies [31] that employ the Thomas-Fermi approximation and include finite-temperature effects at the level of a noninteracting Bose gas ($g = 0$). We have shown that the truncated GPPE can be used effectively to study the gravitational 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 noncondensed 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.

We note that a GP-Poisson system has been used, in 2D and at $T = 0$, in the context of ultracold plasmas, which have electrostatic interactions, instead of gravity; here, the Poisson equation is essentially the same as in our paper; however, in electrostatics, like charges repel, whereas masses have an attractive gravitational interaction. We refer the reader to Refs. [32,33] for a discussion of this 2D system at $T = 0$.

We thank DST, CSIR (India), and the Indo French Center for Applied Mathematics (IFCAM) for their support.

[1] D. O'Dell, S. Giovanazzi, G. Kurizki, and V. M. Akulin, Phys. Rev. Lett. 84, 5687 (2000).

[2] R. Ruffini and S. Bonazzola, Phys. Rev. 187, 1767 (1969).

[3] G. Ingrosso, D. Grasso, and R. Ruffini, Astron. Astrophys. 248, 481 (1991).

[4] P. Jetzer, Phys. Rep. 220, 163 (1992).

[5] W. Hu, R. Barkana, and A. Gruzinov, Phys. Rev. Lett. 85, 1158 (2000).

[6] H. Abdallah, A. Abramowski, F. Aharonian, F. Ait Benkhali, E. O. Angüner, M. Arakawa, M. Arrieta, P. Aubert, M. Backes, A. Balzer et al., Phys. Rev. Lett. 120, 201101 (2018).

[7] C. Kachulis, K. Abe, C. Bronner, Y. Hayato, M. Ikeda, K. Iyogi, J. Kameda, Y. Kato, Y. Kishimoto, Ll. Marti et al., Phys. Rev. Lett. 120, 221301 (2018).

[8] https://physics.aps.org/articles/v11/48?utm_campaign=weekly&utm_medium=email&utm_source=emailalert.

[9] M. Lawson, A. J. Millar, M. Pancaldi, E. Vitagliano, and F. Wilczek, Phys. Rev. Lett. 123, 141802 (2019).

[10] P.-H. Chavanis, Phys. Rev. D 84, 043531 (2011).

[11] P.-H. Chavanis and L. Delfini, Phys. Rev. D 84, 043532 (2011).

[12] P.-H. Chavanis, Phys. Rev. D 94, 083007 (2016).

[13] P.-H. Chavanis, Phys. Rev. D 98, 023009 (2018).
[14] P.-H. Chavanis, *Quantum Aspects of Black Holes*, edited by. X. Calmet (Springer, Cham, 2015), Chap. 6, pp. 151–194.

[15] N. P. Proukakis and B. Jackson, *J. Phys. B* **41**, 203002 (2008).

[16] N. G. Berloff, M. Brachet, and N. P. Proukakis, *Proc. Natl. Acad. Sci. USA* **111**, 4675 (2014).

[17] G. Krstulovic and M. Brachet, *Phys. Rev. E* **83**, 066311 (2011).

[18] V. Shukla, M. Brachet, and R. Pandit, *New J. Phys.* **15**, 113025 (2013).

[19] P. J. E. Peebles, *The Large-Scale Structure of the Universe* (Princeton University Press, Princeton, NJ, 1980).

[20] M. Falco, S. H. Hansen, R. Wojtak, and G. A. Mamon, *Mon. Not. R. Astron. Soc.* **431**, L6 (2013).

[21] M. K.-H. Kiessling, *Adv Appl Math* **31**, 132 (2003).

[22] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevResearch.3.L022016 for details of the following: (a) the gravitational collapse that we obtain for our model (see the Main Paper) at $T = 0$ from the stochastic Ginzburg-Landau-Poisson equation (SGLPE) [Sec. I]; (b) videos from our DNSs of the Gross-Pitaevskii-Poisson equation (GPPE), SGLPE at $T = 0$ and at finite $T$ [Sec. II]; (c) the initial conditions that we use to obtain a binary condensate that comprises a pair of collapsed objects rotating around their center of mass [Sec. III]; and (d) parameters of our simulations [Sec. IV] for the GPPE and the SGLPE at $T = 0$.

[23] D. Gottlieb and S. A. Orszag, *Numerical Analysis of Spectral Methods* (SIAM, Philadelphia, 1977).

[24] H. Y. Schive, T. Chiueh, and T. Broadhurst, *Nat. Phys.* **10**, 496 (2014).

[25] J. Veltmaat, J. C. Niemeyer, and B. Schwabe, *Phys. Rev. D* **98**, 043509 (2018).

[26] P. Mocz, A. Fialkov, M. Vogelsberger, F. Becerra, M. A. Amin, S. Bose, M. Boylan-Kolchin, P.-H. Chavanis, L. Hernquist, L. Lancaster et al., *Phys. Rev. Lett.* **123**, 141301 (2019).

[27] P. Mocz, A. Fialkov, M. Vogelsberger, F. Becerra, X. Shen, V. H. Robles, M. A. Amin, J. Zavala, M. Boylan-Kolchin, S. Bose et al., *Mon. Not. R. Astron. Soc.* **494**, 2027 (2020).

[28] C. Armendariz-Picon and J. T. Neelakanta, *J. Cosmol. Astropart. Phys.* **03** (2014) 049.

[29] S. Latifah, A. Sulaksono, and T. Mart, *Phys. Rev. D* **90**, 127501 (2014).

[30] G. Krstulovic and M. Brachet, *Phys. Rev. Lett.* **106**, 115303 (2011); **107**, 099602 (2011).

[31] T. Harko and E. J. M. Madarassy, *J. Cosmol. Astropart. Phys.* **01** (2012) 020.

[32] H. Sakaguchi and B. A. Malomed, *Phys. Rev. Res* **2**, 033188 (2020).

[33] J. Qin and G. Dong, and B. A. Malomed, *Phys. Rev. A* **94**, 053611 (2016).