# Kinetic relaxation and nucleation of Bose stars in self-interacting wave
dark matter

Mudit Jain<EMAIL_ADDRESS>Department of Physics and Astronomy, Rice
University, Houston, Texas 77005, U.S.A.  Theoretical Particle Physics and
Cosmology, King’s College London, Strand, London, WC2R 2LS, United Kingdom
Wisha Wanichwecharungruang<EMAIL_ADDRESS>Jonathan Thomas
<EMAIL_ADDRESS>Department of Physics and Astronomy, Rice University,
Houston, Texas 77005, U.S.A.

###### Abstract

We revisit kinetic relaxation and soliton/Boson star nucleation in fuzzy
scalar dark matter featuring short-ranged self-interactions $\mathcal{H}_{\rm
int}=-\lambda|\psi|^{4}/2m^{2}$, alongside gravitational self-interactions. We
map out the full curve of nucleation timescale for both repulsive
($\lambda<0$) and attractive ($\lambda>0$) short-ranged self-interaction
strength, and in doing so reveal two new points. Firstly, besides the two
usual terms, $\propto G^{2}$ and $\propto\lambda^{2}$, in the total relaxation
rate $\Gamma_{\rm relax}$, there is an additional cross term $\propto
G\lambda$ arising due to interference between gravitational and short-ranged
self-interaction scattering amplitudes. This yields a critical repulsive
interaction strength $\lambda_{\rm cr}\simeq-2\pi Gm^{2}/v_{0}^{2}$, at which
the relaxation rate is smallest and serves as the transition point between
typical net attractive self-interaction ($\lambda\gtrsim\lambda_{\rm cr}$),
and net repulsive self-interaction ($-\lambda\gtrsim-\lambda_{\rm cr}$).
Secondly, while in the net attractive regime, nucleation time scale is similar
to inverse relaxation time scale $\tau_{\rm nuc}\sim\Gamma^{-1}_{\rm relax}$,
in the net repulsive regime nucleation occurs at a delayed time $\tau_{\rm
nuc}\sim(\lambda/\lambda_{\rm cr})\Gamma^{-1}_{\rm relax}$. We confirm our
analytical understanding by performing 3D field simulations with varying
average mass density $\bar{\rho}$, box size $L$ and grid size $N$.

††preprint: APS/123-QED

###### Contents

  1. I Introduction
  2. II Model
  3. III Wave kinetics and relaxation
  4. IV Nucleation and behavior of solitons
    1. IV.1 Nucleation
    2. IV.2 Eventual behavior
      1. IV.2.1 Attractive short-ranged interactions $\lambda>0$
      2. IV.2.2 Repulsive short-ranged interactions $\lambda<0$
  5. V Field simulations
    1. V.0.1 Net attractive interactions ($\lambda\gtrsim\lambda_{\rm cr}$).
    2. V.0.2 Net repulsive interactions ($-\lambda\gtrsim-\lambda_{\rm cr}$)
  6. VI Summary and Discussion
    1. VI.1 Comparison with earlier work
    2. VI.2 Implications
  7. A Statistical convergence
  8. B Peculiar appearance of high frequency modes for the repulsive case

##  I Introduction

Understanding the nature of dark matter (DM) is one of the main quests of
modern cosmology. It could be multi-faceted in the sense that there are many
degrees of freedom in the whole dark sector, for instance the String theory
Axiverse [1, 2, 3] or other confined sector(s) (e.g. see [4, 5], and also
[6]). Or it could be that there is a dominant degree of freedom, such as the
QCD axion [7, 8, 9, 10, 11, 12], that comprises all (or most) of the dark
matter. Furthermore, while the DM appears to interact only gravitationally
with the Standard Model degrees of freedom (or very weakly if it does
otherwise), it can still have appreciable non-gravitational self-interactions
(nGSI) besides the usual gravitational self-interactions (GSI). Such is the
case even for the above mentioned examples.

For bosonic particles (of any integer spin) and high enough occupation
numbers, which is indeed the case for particle masses below a few eV,
classical description of the associated field suffices and the dynamics is
described by a non-linear Schrödinger equation in the non-relativistic regime.
The non-linear Schrödinger equation entails novel wave dynamics owing to the
de-Broglie scale becoming manifestly important. As a few examples, suppression
of structure on small scales [13, 14], turbulence [15], superradiance [16],
vortices [17], bound states called solitons/Bose stars [18, 19, 20, 21, 17,
22, 23, 24, 25, 26, 27, 28], interference patterns [29, 28, 30], field
correlation scales depending upon the nature of self-interaction [31], etc.
For comprehensive recent reviews in the case of scalar DM, see [32, 33].

Of particular interest to us in this paper, is the phenomenon of kinetic
relaxation and associated nucleation of Bose stars within a bath of DM waves
[34, 35, 36, 37, 38, 39, 40, 41, 42]. The term “kinetic” implies two key
aspects: (a) The self-interactions in the field are small, allowing wave modes
to freely evolve (at leading order) with the non-relativistic dispersion
relation $\omega_{k}=k^{2}/2m$. This enables a kinetic treatment of the mode
occupation number function; (b) the size of the ‘box’ ($\sim$ the size of a DM
halo for practical purposes), is much larger than the typical fluctuation
scale $\ell_{\rm dB}\sim\pi/\bar{k}\sim\pi/(mv_{0})$ in the bath of DM waves.
The process of kinetic relaxation is attributed to these self-interactions of
the DM field which although small, over large time scales $\tau_{\rm
relax}\gg\omega_{\bar{k}}^{-1}$ drive the occupation number function to
develop increasing support towards smaller wavenumbers ${\bm{k}}\rightarrow
0$. See [34] and [41] for a relevant discussion for the cases of point-like
quartic self-interactions and gravitational self-interactions respectively.
Once enough particles condense into lower momentum states, their collective
net attractive self-interaction becomes strong enough to counter balance their
wave pressure resulting in the nucleation of a Bose star.

In this paper, we focus on investigating kinetic relaxation and subsequent
Bose star nucleation for a single scalar Schrödinger field with both GSI and
point-like quartic nGSI. Employing wave-kinetic Boltzmann analysis and 3D
simulations, we demonstrate the presence of a previously overlooked cross term
$\propto G\lambda$ in the rate of relaxation $\Gamma_{\rm relax}$. (Here $G$
denotes Newton’s constant and $\lambda$ represents the point-like nGSI
strength). It arises due to interference between the gravitational and point-
like self-interaction scattering amplitudes. The presence of this cross term
gives rise to a critical nGSI (repulsive) strength $\lambda_{\rm
cr}\simeq-(2\pi G)m^{2}/v_{0}^{2}$, at which the rate of relaxation reaches
its minimum value (corresponding to maximum nucleation time). This critical
value also serves as the transition point from typical net (contributions from
both gravity and short-ranged self-interactions) attractive to repulsive self-
interactions.

Because of the presence of gravitational self-interaction, kinetic relaxation
is generally accompanied with nucleation of spatially localized clumps/Bose
stars, with their nucleation times dependent on the nature of the short-ranged
self-interactions – attractive or repulsive. For $\lambda\gtrsim\lambda_{\rm
cr}$, the net typical self-interaction is attractive, and nucleation happens
quickly after relaxation. On the other hand for $\lambda\lesssim\lambda_{\rm
cr}$, the net typical self-interaction is repulsive and nucleation gets
delayed. We will study relaxation and nucleation of Bose stars, and also
discuss their eventual fate.111Following conventional nomenclature, we shall
use the words relaxation and condensation interchangeably, but it is to be
stressed that nucleation (of a bound state) is not always equivalent to
relaxation/condensation. As we shall see, it is equivalent to the other two in
the net attractive regime, whereas different in the net repulsive regime.
However we will not dwell into a careful analysis of the growth rate of these
nucleated stars. See [43, 44, 45] for the gravity only ($\lambda=0$) case.

The rest of the paper is organized as follows: Starting with the basic model
of fuzzy scalar DM carrying both GSI and point-like nGSI in sec. II, we
describe the associated wave kinetic Boltzmann equation for the evolution of
the occupation number function in sec. III. Highlighting the presence of the
cross term (that gives rise to $\lambda_{\rm cr}$), we estimate the total rate
of kinetic relaxation/condensation. In sec. IV we discuss the two cases of
$\lambda\gtrsim\lambda_{\rm cr}$ and $-\lambda\gtrsim-\lambda_{\rm cr}$ and
write down the associated nucleation time scales of spatially localized bound
objects. In sec. V we discuss our 3D simulations and compare our analytical
estimates with them. We also discuss eventual behavior of Bose clumps observed
in simulations. Finally in sec.VI, we summarize our work and also compare our
results with the existing literature on this subject. In appendix A we discuss
statistical convergence of our simulations, and in appendix B we discuss a
peculiarity observed in the case of repulsive short-ranged self-interactions,
over longer time scales as compared to nucleation.

Conventions: Unless stated otherwise, we will work in units where $\hbar=c=1$.

##  II Model

Ignoring Hubble flow (for we are interested in sufficiently sub-horizon
dynamics), the evolution of the cold/non-relativistic fuzzy scalar dark matter
with both GSI and short-ranged quartic nGSI, can be described using mean field
theory. The dark matter field $\psi$ obeys the following non-linear
Schrödinger (Gross-Pitaevskii) equation:

$\displaystyle i\frac{\partial}{\partial
t}\psi=-\frac{1}{2m}\nabla^{2}\psi+\psi\Bigl{(}4\pi
Gm^{2}\nabla^{-2}_{/0}-\frac{\lambda}{m^{2}}\Bigr{)}\psi^{\ast}\psi\,.$ (1)

Here $G$ is the Newton’s constant, and $\lambda$ is the point like self-
interaction strength. In our convention, $\lambda>0$ and $\lambda<0$ dictate
attractive and repulsive self-interaction respectively. To obtain the form Eq.
(1), we have plugged the self-gravitational potential, $\Phi=4\pi
G\,\nabla^{-2}(m\psi^{\ast}\psi-\bar{\rho})\equiv 4\pi
Gm\nabla^{-2}_{/0}\psi^{\ast}\psi$, in the usual Schrödinger-Poisson system of
equations. The $\nabla^{-2}_{/0}$ denotes exclusion of the homogeneous part of
the number density field $\psi^{\ast}\psi$. In Fourier space with the
decomposition
$\psi({\bm{x}},t)=(2\pi)^{-3}\int\mathrm{d}{\bm{k}}\,e^{-i{\bm{k}}\cdot{\bm{x}}}\,\Psi_{\bm{k}}(t)$,
the Schrödinger equation becomes

$\displaystyle i\dot{\Psi}_{\bm{k}}$
$\displaystyle=\frac{k^{2}}{2m}\Psi_{\bm{k}}+\int\frac{\mathrm{d}{\bm{p}}}{(2\pi)^{3}}\frac{\mathrm{d}{\bm{q}}}{(2\pi)^{3}}\frac{\mathrm{d}{\bm{\ell}}}{(2\pi)^{3}}\mathcal{T}_{\bm{k},\bm{p},\bm{q},\bm{\ell}}\,\Psi^{\ast}_{\bm{p}}\Psi_{\bm{q}}\Psi_{\bm{\ell}}$
$\displaystyle\qquad\qquad\qquad\qquad\times\,(2\pi)^{3}\delta^{(3)}(\bm{k}+\bm{p}-\bm{q}-\bm{\ell})\,,$
(2)

where

$\displaystyle\mathcal{T}_{\bm{k},\bm{p},\bm{q},\bm{\ell}}=$
$\displaystyle-\frac{4\pi
Gm^{2}}{|{\bm{k}}-{\bm{\ell}}|^{2}}-\frac{\lambda}{m^{2}}\,,$ (3)

and it is understood that ${\bm{k}}\neq{\bm{\ell}}\neq 0$ in the above. For
later convenience, it is also useful to write down the Hamiltonian density (in
physical space) for the mean field $\psi$:

$\displaystyle\mathcal{H}=\frac{1}{2m}|\nabla\psi|^{2}+m\Phi|\psi|^{2}-\frac{\lambda}{2m^{2}}|\psi|^{4}\,.$
(4)

Here the different terms in the above can be attributed to the wave-pressure
$\mathcal{H}_{\rm wp}=|\nabla\psi|^{2}/2m$, gravitational self-interaction
$\mathcal{H}_{\rm gr}=m\Phi|\psi|^{2}$, and short-ranged self-interaction
$\mathcal{H}_{\rm self}=-\lambda|\psi|^{4}/2m^{2}$.

The Gross-Pitaevskii (GP) equation, being non-linear, renders it difficult to
analyze and study the evolution of the $\psi$ field in generality. However for
the purposes of kinetic relaxation leading to nucleation of localized Bose
clumps however, wave kinetic Boltzmann analysis can be performed which we
discuss next. To test and verify our analytical understanding, we perform 3D
field simulations which we discuss in a later section.

##  III Wave kinetics and relaxation

For kinetic relaxation in wave dynamics, we can study the evolution of the
mode occupation number function $f_{\bm{k}}=|\Psi_{\bm{k}}|^{2}/V$ ($V$ is the
volume), which is nothing but the Fourier transform of the $2$-point volume
averaged field correlator
$\zeta({\bm{x}},t)=V^{-1}\int\mathrm{d}{\bm{y}}\,\psi^{\ast}(\bm{y},t)\,\psi(\bm{y}+\bm{x},t)$.
Under random phase approximation with weak interactions, the relevant wave-
kinetic Boltzmann equation can be derived. See for instance [46]. For a
derivation for the general case of arbitrary number of fields and $2$ body
interactions, see [41]. For the scalar case at hand, characterizing the
dependence of the occupation number functions on wavenumbers as
$f_{\bm{k}/m}$, the wave-kinetic equation takes the familiar form

$\displaystyle\frac{\partial f_{\bm{k}/m}}{\partial
t}=\int\frac{\mathrm{d}{\bm{p}}}{(2\pi)^{3}}\,\mathrm{d}\sigma_{{\bm{k}}+{\bm{p}}\rightarrow{\bm{q}}+{\bm{\ell}}}\,|{\bm{v}}-\tilde{\bm{v}}|\,\Biggl{[}(f_{\bm{k}/m}+f_{\bm{p}/m})f_{\bm{q}/m}f_{\bm{\ell}/m}-(f_{\bm{q}/m}+f_{\bm{\ell}/m})f_{\bm{k}/m}f_{\bm{p}/m}\Biggr{]}\,,$
$\displaystyle{\rm
where}\quad\mathrm{d}\sigma_{{\bm{k}}+{\bm{p}}\rightarrow{\bm{q}}+{\bm{\ell}}}=\frac{1}{2|{\bm{v}}-\tilde{\bm{v}}|}\frac{\mathrm{d}{\bm{q}}}{(2\pi)^{3}}\frac{\mathrm{d}{\bm{\ell}}}{(2\pi)^{3}}\Bigl{(}\mathcal{T}_{\bm{k},\bm{p},\bm{q},\bm{\ell}}+\mathcal{T}_{\bm{k},\bm{p},\bm{\ell},\bm{q}}\Bigr{)}\Bigl{(}\mathcal{T}_{\bm{k},\bm{p},\bm{q},\bm{\ell}}+\mathcal{T}_{\bm{k},\bm{p},\bm{\ell},\bm{q}}\Bigr{)}^{\ast}\,\times$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad(2\pi)^{4}\delta^{(3)}(\bm{k}+\bm{p}-\bm{q}-\bm{\ell})\,\delta(E_{\bm{k}}+E_{\bm{p}}-E_{\bm{q}}-E_{\bm{\ell}})\,.$
(5)

Here ${\bm{v}}={\bm{k}}/m$ and $\tilde{\bm{v}}={\bm{p}}/m$ are the incoming
“velocities” in the $2$-wave interaction, and the quantities in the
$1$-dimensional Dirac delta function are the free wave energies
$E_{\bm{k}}=k^{2}/2m$. The quantity $\mathrm{d}\sigma$ is the effective
differential cross section. The cubic nature of the terms in the right hand
side bracket ($\sim f_{\bm{x}}f_{\bm{y}}f_{\bm{z}}$), usually understood as
Boltzmann enhancement terms, arise due to the wave-mechanical nature of the
system (1) and are crucial for the phenomenon of Bose condensation. Last but
not the least, it is the form of the differential cross section that appears
in the wave-kinetic equation, $\sim|\mathcal{T}|^{2}$, that is of utmost
importance for our discussion. The scattering amplitudes due to the different
kinds of $2$-body interactions (here gravity and point-like self-
interactions), are added first and then squared: What appears in the
differential cross section is $|\mathcal{T}|^{2}$ where
$\mathcal{T}=\mathcal{T}_{G}+\mathcal{T}_{\lambda}$ (c.f. Eq. (3)), and
$|\mathcal{T}_{G}+\mathcal{T}_{\lambda}|^{2}\neq|\mathcal{T}_{G}|^{2}+|\mathcal{T}_{\lambda}|^{2}$
(since both $\mathcal{T}_{G}\propto-4\pi Gm^{2}$ and
$\mathcal{T}_{\lambda}=-\lambda/m^{2}$ are real). This can be attributed to
the wave dynamical nature of the GP system. The above equation (III), after
integration over the Dirac deltas, can be re-written in terms of the incoming
and outgoing relative velocities ${\bm{u}}=u\hat{\bm{n}}$ and
${\bm{u}}^{\prime}=u\hat{\bm{n}}^{\prime}$ respectively222Note that the
magnitude of the relative velocity does not change in an elastic collision,
i.e. $|{\bm{u}}|=|{\bm{u}}^{\prime}|\equiv u$., by redefining
${\bm{p}}/m={\bm{k}}/m-u\hat{{\bm{n}}}$ and
${\bm{q}}/m={\bm{\ell}}/m-u\hat{\bm{n}}^{\prime}$:

$\displaystyle\frac{\partial f_{\bm{v}}}{\partial t}$
$\displaystyle=m^{3}\int\frac{\mathrm{d}{\bm{u}}}{(2\pi)^{3}}\,\mathrm{d}\sigma\,u\Biggl{[}(f_{\bm{v}}+f_{\tilde{\bm{v}}})f_{\tilde{\bm{v}}-{\bm{w}}}f_{\bm{v}+\bm{w}}-(f_{\tilde{\bm{v}}-\bm{w}}+f_{\bm{v}+\bm{w}})f_{\bm{v}}f_{\tilde{\bm{v}}}\Biggr{]}\,,$
$\displaystyle{\rm where}\quad\mathrm{d}\sigma$
$\displaystyle=\frac{\mathrm{d}\Omega_{n^{\prime}}}{32\pi^{2}m^{2}}\Biggl{[}\Bigl{(}\frac{16\pi
m^{2}G}{u^{2}|\hat{\bm{n}}^{\prime}-\hat{\bm{n}}|^{2}}+\lambda\Bigr{)}^{2}+\Bigl{(}\frac{16\pi
m^{2}G}{u^{2}|\hat{\bm{n}}^{\prime}+\hat{\bm{n}}|^{2}}+\lambda\Bigr{)}^{2}+2\Bigl{(}\frac{16\pi
m^{2}G}{u^{2}|\hat{\bm{n}}^{\prime}-\hat{\bm{n}}|^{2}}+\lambda\Bigr{)}\Bigl{(}\frac{16\pi
m^{2}G}{u^{2}|\hat{\bm{n}}^{\prime}+\hat{\bm{n}}|^{2}}+\lambda\Bigr{)}\Biggr{]}\,.$
(6)

Let us briefly discuss the different terms in the differential cross section
explicitly. Broadly speaking, there are two types of interference terms that
arise. One is the interference between the $t$ and $u$ channels (relevant
mainly for the gravitational interaction), and the second is the interference
between the two different types of interactions (gravitational and short-
ranged). See fig. 1 for a pictorial representation.

For GSI only ($\lambda=0$) case, the contribution from the $t$ and $u$
channels are the first two terms $\propto
G^{2}\,|\hat{\bm{n}}^{\prime}-\hat{\bm{n}}|^{-4}$ and $\propto
G^{2}\,|\hat{\bm{n}}^{\prime}+\hat{\bm{n}}|^{-4}$, whereas the second term
$\propto
G^{2}\,|\hat{\bm{n}}^{\prime}-\hat{\bm{n}}|^{-2}\,|\hat{\bm{n}}^{\prime}+\hat{\bm{n}}|^{-2}$
is due to their mutual interference (as also discussed in [41]). Note that the
sole contributions from the $t$ and $u$ channels are identical: The full
integral with the $|\hat{\bm{n}}^{\prime}+\hat{\bm{n}}|^{-4}$ term is
identical to that with the $|\hat{\bm{n}}^{\prime}-\hat{\bm{n}}|^{-4}$ term.
The sole contributions give rise to the Rutherford scattering cross section,
carrying a logarithmic IR divergence (aka the Coulomb logarithm), while the
interference term becomes sub-dominant in the large log limit and can be
omitted.

For nGSI only ($G=0$) case, contributions from $t$ and $u$ channels are
identical to their mutual interference one, and goes as $\lambda^{2}$. This is
simply due to the interaction being a contact/point interaction.

Importantly when both of the interactions are present, their respective
scattering amplitudes (for either of the two channels) are added first and
then squared. All the terms $\propto G\lambda$, while giving identical
contributions, characterize the interference between the two types of
interactions. Splitting the contributions from GSI, nGSI, and their
interference, we have the following wave-kinetic Boltzmann equation

$\displaystyle\frac{\partial f_{\bm{v}}}{\partial t}$
$\displaystyle=\mathcal{C}_{\rm GSI}+\mathcal{C}_{\rm cross}+\mathcal{C}_{\rm
nGSI}$ $\displaystyle\qquad\qquad{\rm where}\quad\mathcal{C}_{\rm GSI}$
$\displaystyle=\frac{\Lambda(4\pi
G)^{2}m^{5}}{4\pi}\nabla_{{v}^{i}}\Biggl{[}\frac{1}{2}\nabla_{v^{j}}f_{{\bm{v}}}\int\frac{\mathrm{d}\tilde{\bm{v}}}{(2\pi)^{3}}\,f_{\tilde{\bm{v}}}\,\frac{\delta_{ij}-\hat{u}_{i}\hat{u}_{j}}{u}\,f_{\tilde{\bm{v}}}+f_{{\bm{v}}}\,f_{{\bm{v}}}\int\frac{\mathrm{d}\tilde{\bm{v}}}{(2\pi)^{3}}\,\frac{\hat{u}_{i}}{u^{2}}\,f_{\tilde{\bm{v}}}\Biggr{]},$
$\displaystyle\mathcal{C}_{\rm cross}$ $\displaystyle=\frac{(4\pi G)\lambda
m^{3}}{4\pi}\int\frac{\mathrm{d}\Omega_{n}}{4\pi}\frac{\mathrm{d}u\,u^{2}}{2\pi^{2}}\frac{\mathrm{d}\Omega_{n^{\prime}}}{4\pi}\frac{u}{|{\bm{w}}|^{2}}\Biggl{[}(f_{\bm{v}}+f_{\tilde{\bm{v}}})f_{\bm{v}+{\bm{w}}}f_{\tilde{{\bm{v}}}-{\bm{w}}}-(f_{\bm{v}+{\bm{w}}}+f_{\tilde{\bm{v}}-{\bm{w}}})f_{\bm{v}}f_{\tilde{\bm{v}}}\Biggr{]},$
$\displaystyle\mathcal{C}_{\rm nGSI}$
$\displaystyle=\frac{\lambda^{2}m}{2\pi}\int\frac{\mathrm{d}\Omega_{n}}{4\pi}\frac{\mathrm{d}u\,u^{2}}{2\pi^{2}}\frac{\mathrm{d}\Omega_{n^{\prime}}}{4\pi}\,u\Biggl{[}(f_{\bm{v}}+f_{\tilde{\bm{v}}})f_{\bm{v}+{\bm{w}}}f_{\tilde{{\bm{v}}}-{\bm{w}}}-(f_{\bm{v}+{\bm{w}}}+f_{\tilde{\bm{v}}-{\bm{w}}})f_{\bm{v}}f_{\tilde{\bm{v}}}\Biggr{]}\,.$
(7)

Here ${\bm{w}}=u(\hat{\bm{n}}^{\prime}-\hat{\bm{n}})/2$, and
$\Lambda=\log(mv_{0}L)$ is the aforementioned Coulomb logarithm with $v_{0}$
and $L$ equal to typical velocity and box size (or halo size for physical
considerations) respectively. While the cross term and nGSI term follow
straightforwardly from Eq. (III), the Rutherford scattering collision term
$\mathcal{C}_{\rm GSI}$ is obtained after an eikonal approximation and was
derived explicitly in [41]. Also see [46, 47] for the same equation for a
scalar field.

Now in order to get a typical estimate for the total relaxation rate
$\Gamma_{\rm relax}\equiv\frac{1}{f_{\bm{v}}}\frac{\partial
f_{\bm{v}}}{\partial t}$, we can replace different quantities in the three
collision terms with their appropriate scalings. Replacing angular volume
$\int\mathrm{d}\Omega\rightarrow 4\pi$, typical relative velocity
$|\hat{\bm{n}}^{\prime}-\hat{\bm{n}}|\rightarrow\sqrt{2}$, velocity derivative
$\nabla_{v}\rightarrow 1/v_{0}$, velocity integral
$\int\mathrm{d}u\,u^{n-1}\rightarrow v_{0}^{n}/n$, and finally the occupation
number function
$f_{\bm{v}}\rightarrow(2\pi)^{3/2}\bar{\rho}/(m^{4}v_{0}^{3})$, the total
relaxation rate is parameterized as

$\displaystyle\Gamma_{\rm relax}\simeq\alpha_{1}\frac{(4\pi
G)^{2}\bar{\rho}^{2}\Lambda}{4m^{3}v_{0}^{6}}+\alpha_{12}\frac{(4\pi
G)\lambda\bar{\rho}^{2}}{m^{5}v_{0}^{4}}+\alpha_{2}\frac{\lambda^{2}\bar{\rho}^{2}}{m^{7}v_{0}^{2}}\,.$
(8)

Our scaling of the occupation number is dictated by Gaussian initial condition
(see Eq. (V) ahead) which we shall use to perform simulations, described in
the next section. In general, $\alpha_{1}$, $\alpha_{12}$, and $\alpha_{2}$
are positive $\mathcal{O}(1)$ coefficients that would depend on the specific
initial conditions. Eq. (8) is our master formula for the relaxation rate. The
value of $\lambda$ around which the relaxation rate becomes smallest, is
easily estimated to be

$\displaystyle\lambda_{\rm cr}=-\beta\frac{2\pi G\,m^{2}}{v_{0}^{2}}\sim
10^{-57}\left(\frac{10^{-4}}{v_{0}}\right)^{2}\left(\frac{m}{10^{-5}\,{\rm
eV}}\right)^{2}\,,$ (9)

where $\beta=\alpha_{12}/\alpha_{2}\sim\mathcal{O}(1)$, and the associated
(minimum) rate is333In the large Coulomb limit (relevant for realistic
scenarios), $\Lambda=\log(mv_{0}L)\gg 1$, and the rate is always positive.

$\displaystyle\Gamma_{\rm relax}\left(\lambda_{\rm cr}\right)\simeq\frac{(4\pi
G)^{2}\bar{\rho}^{2}}{4m^{3}v_{0}^{6}}\left(\alpha_{1}\Lambda-\frac{\alpha_{12}^{2}}{\alpha_{2}}\right)\,.$

Notice that this critical value of $\lambda_{\rm cr}$, can also be obtained
from the GP equation (1) by balancing the gravitational term with the self-
interaction term together with replacing the exchange momenta by its typical
value $|{\bm{k}}-{\bm{\ell}}|^{2}\sim 2(mv_{0})^{2}$. This criticality marks
the transition point from attractive to repulsive net typical self-
interactions: For $\lambda\gtrsim\lambda_{\rm cr}$, typical interactions
within the bath of DM waves are attractive since typical $\mathcal{T}$ is
negative, whereas for $-\lambda\gtrsim-\lambda_{\rm cr}$ they are repulsive
since typical $\mathcal{T}$ is positive.

Figure 1: Pictorial / Feynman graph representation of all the terms appearing
in the differential cross section in Eq. (III). The total contribution to the
interaction rate (left hand side of the equality), is the square of the sum of
both gravitational amplitude (top two graphs) and point self-interaction
amplitude (bottom graph $\times 2$). For gravitational interaction, there are
two distinct channels ($t$ and $u$). Their mutual interference, as compared to
their sole contributions, becomes subdominant in the large log limit (leading
to Rutherford scattering result). More importantly, interference between
scattering amplitudes of the two different interactions matters. See main text
for details.

##  IV Nucleation and behavior of solitons

###  IV.1 Nucleation

In general, the process of kinetic relaxation is characterized by an
increasing support of the occupation number function $f_{\bm{k}}$ at vanishing
wavenumber. (For instance see [34] and [41] for discussions of nGSI and GSI
cases respectively). This implies increasing field correlation over larger
length scales with diminishing density fluctuations, i.e. field
homogenization. A heuristic understanding of the subsequent nucleation of a
spatially localized and bound clump, can perhaps be gained most easily from a
particle physics perspective, together with recalling that $\lambda_{\rm cr}$
also marks the transition from typical net attractive self-interaction to
typical net repulsive self-interaction. As particles lose kinetic energy on
account of self-interactions and move towards smaller momenta (condensate
state), there comes a time when within some region, the collective net
potential (due to both self-gravitational and short-ranged interactions)
becomes comparable to wave pressure. The time scale of this process is nothing
but the inverse relaxation rate Eq. (8), which in the case of net attractive
self-potential $\lambda\gtrsim\lambda_{\rm cr}$, leads to ‘immediate locking’
of such a region into a bound clump (having negative energy). That is,
$\tau_{\rm nuc}\simeq\Gamma_{\rm relax}^{-1}$. Strictly speaking, this can be
taken as a definition of $\tau_{\rm nuc}$ with $\Gamma_{\rm nuc}=\Gamma_{\rm
relax}$, in which case the different $\alpha$ constant coefficients in the
rate Eq. (8), are understood as such.

On the other hand for $-\lambda\gtrsim-\lambda_{\rm cr}$, relaxation cannot
immediately lead to nucleation of a localized bound clump. This is because the
net typical interaction is repulsive: the collective self-potential within
density fluctuation regions is not binding yet. Over time though, more
particles get driven towards the condensate phase, and eventually there arises
a potential for a bound object to nucleate (within which net gravity can now
compensate for both repulsive short-ranged interaction and wave-pressure).
This gives $\tau_{\rm nuc}>\Gamma_{\rm relax}^{-1}$. In general, we can
therefore write the following

$\displaystyle\tau_{\rm nuc}\simeq\frac{1}{\Gamma_{\rm
relax}}\begin{cases}1\qquad&\lambda\gtrsim\lambda_{\rm cr}\\\
h(\lambda)\qquad&-\lambda\gtrsim-\lambda_{\rm cr}\end{cases}$ (10)

where $h(\lambda)$ is a threshold function (or the delay factor), that relates
nucleation times to relaxation rates. As mentioned earlier, relaxation means
field homogenization, and we expect the rate at which the system relaxes to be
comparable to the rate at which density fluctuations decrease. The delay
factor can then be estimated as the ratio of typical density fluctuation at
relaxation, to that at nucleation, $h\sim\delta\rho_{\rm
relax}/\delta\rho_{\rm nuc}$. While we expect this to be order unity for net
attractive case (first case of Eq. 10), for large repulsive strengths it
should increase with increasing $-\lambda$. Below we estimate this scaling.

Consider a region of typical size $\sim(mv_{0})^{-1}$ where the field would
have ‘locked’ itself into a bound configuration upon relaxation/condensation,
were the net potential was binding. However this is not the case yet, and we
may balance the typical self-interaction energy density (mostly due to short-
ranged interactions) $\mathcal{H}_{\rm self}\sim-\lambda\,\delta\rho_{\rm
relax}^{2}/2m^{4}$, with the wave pressure within $\mathcal{H}_{\rm wp}\sim
v_{0}^{2}\delta\rho_{\rm relax}/2$. This gives $\delta\rho_{\rm relax}\sim
m^{4}v_{0}^{2}/\lambda$. As relaxation continues (meaning more particles are
driven towards low momenta state), the value of both density fluctuations
$\delta\rho$ and typical size of fluctuation regions $(mv)^{-1}$ change. The
former decreases and the latter increases so as to maintain $\mathcal{H}_{\rm
self}\sim\mathcal{H}_{\rm wp}$. Then, nucleation is expected to occur when
gravity can compensate for both the wave pressure and repulsive short-ranged
self-interaction. That is, we may balance (the magnitudes of) all the three
energy densities, $\mathcal{H}_{\rm wp}\sim v^{2}_{\rm nuc}\delta\rho_{\rm
nuc}/2$, $|\mathcal{H}_{\rm gr}|\sim 2\pi G\delta\rho^{2}_{\rm nuc}/(mv_{\rm
nuc})^{2}$, and $\mathcal{H}_{\rm self}\sim-\lambda\delta\rho^{2}_{\rm
nuc}/2m^{4}$, to give $\delta\rho_{\rm nuc}\sim m^{2}v^{4}_{\rm nuc}/(4\pi G)$
and $v_{\rm nuc}\sim(4\pi Gm^{2}/(-\lambda))^{1/2}$. Eliminating $v_{\rm nuc}$
from $\delta\rho_{\rm nuc}$ gives the following estimate for the delay factor

$\displaystyle h(\lambda)\sim\frac{\delta\rho_{\rm relax}}{\delta\rho_{\rm
nuc}}\rightarrow\alpha_{3}\left(\frac{\lambda}{\lambda_{\rm cr}}\right)\,.$
(11)

Here we have inserted another constant coefficient $\alpha_{3}$ that depends
upon the initial conditions. Through simulations, we will confirm our estimate
Eq. (10) (together with Eq. (8) and Eq. (11)), and also extract the different
$\alpha$ coefficients for Gaussian initial conditions.

###  IV.2 Eventual behavior

Once a spatially-localized Bose clump/soliton emerges, its subsequent
evolution and long term dynamics depends on whether the short-ranged self-
interactions are attractive or repulsive. The full spectrum of such solitons
is extensively discussed in the literature. See e.g. [20, 48, 49, 50]. To
recapitulate some of the basic points that may suffice for our purposes,
consider the energy landscape for objects of radius $r_{s}$ and mass $M_{s}$
in the theory. Using Eq. (4), the wave pressure energy, self-gravitational
potential energy, and short-ranged self-interaction potential energy are
$H_{\rm wp}=aM_{s}/(m^{2}r_{s}^{2})$, $H_{\rm gr}=-b(4\pi
G)M_{s}^{2}/(r_{s})$, and $H_{\rm self}=-c\lambda M_{s}^{2}/(m^{4}r_{s}^{3})$
respectively, with $a$, $b$, and $c$ some positive coefficients that depend
upon the exact profile of the object. The total energy is the sum of all
three.

####  IV.2.1 Attractive short-ranged interactions $\lambda>0$

For this case, the energy vs radius curve (for a given mass $M_{s}$) has a
local minima that corresponds to quasi stable negative energy (bound) states /
solitons. It is separated from the runaway behavior towards small radii,
$\sim-\lambda/r_{s}^{3}$, by a barrier whose height decreases with increasing
$M_{s}$. The barrier disappears at a critical mass $M_{s,\rm
crit}\propto(\lambda G)^{-1/2}$, beyond which the theory does not admit any
quasi-stable bound states anymore. Starting in the kinetic regime and upon
relaxation, a quasi-stable Bose clump nucleates and starts to accrete mass
from its surroundings. Ultimately once it accumulates enough mass such that
the energy barrier gets sufficiently low, and/or it ‘breathes’ rapidly enough
so as to be able to probe beyond the energy barrier, it ‘collapses’. This is
because the region now prefers to lower its energy by transitioning on the
runaway $\sim r_{s}^{-3}$ curve. This is sometimes referred to as “Bosenova”
(owing to its analogy with a type II supernova). While subsequent evolution of
the object beyond this criticality requires a fully relativistic analysis and
has been pursued in the literature [51] (also see [52] for an associated
astrophysical phenomenology), the evolution leading upto this criticality (and
even beyond until the wave pressure starts to become comparable to rest mass
energy) is well captured by the non-relativistic treatment. For large
attractive strengths $\lambda\gtrsim-\lambda_{\rm cr}$, the barrier is less
significant and the objects quickly collapses upon nucleation. In our
simulations we indeed observe this phenomenon (see fig. 3 in sec. V ahead).

####  IV.2.2 Repulsive short-ranged interactions $\lambda<0$

In the repulsive scenario there is no runaway domain since the energy for low
radii is now $-\lambda/r_{s}^{3}>0$. This renders the previous local minima
stable (hence now global minima), corresponding to bound soliton states. The
critical mass $M_{s,\rm crit}\propto(-\lambda G)^{-1/2}$ serves as the
transition point into the Thomas-Fermi regime [53, 20, 54]. This is where the
mass of solitons gets large enough to admit comparable amounts of self-
gravitational and short-ranged interaction energy densities, with gradient
pressure becoming sub-dominant. As a result, the radius starts to approach a
constant $r_{\rm s}\sim\sqrt{-\lambda}/(m^{2}\sqrt{4\pi G})$ (with the mass
dependent correction term dying out as $\sim M_{s}^{-1}$). Up until the mass
becomes sufficiently large where $GM_{s,\rm relv}\sim r_{s}$ and relativistic
effects start to become important (see [55, 56]), the theory then admits a set
of “Chandrasekhar solitons” with masses ranging anywhere between $M_{s,\rm
crit}$ and $M_{s,\rm relv}$, and radii approximately around $r_{\rm
s}\sim\sqrt{-\lambda}/(m^{2}\sqrt{4\pi G})$.444The reason we call them
“Chandrasekhar” solitons (also see [57]) is because of the scaling of their
maximum mass $M_{s,\rm relv}$. It behaves similarly to that of the
Chandrasekhar limit for degenerate stars $M\propto G^{-3/2}m^{-2}$, and can be
attributed to the fact that the Fermi pressure essentially gets replaced with
repulsive short-ranged self-interactions [53, 20, 54, 58].

Though the theory admits these stable Chandrasekhar solitons, understanding
their evolution and long term behavior within the bath of DM waves is crucial,
and has been extensively studied in the literature. See [59, 60, 61, 62, 63]
for simulation setups using the fluid/Madelung equations instead of the
Schrödinger field equation for scalar wave dark matter (with repulsive short-
ranged self-interaction). For our Fourier split simulation technique (which we
discuss in the next section), we find that over longer time scales (after the
nucleation of Bose clumps), the system reaches some sort of criticality when
high frequency modes (near cutoff) start to appear in the simulation box. This
leads to a breakdown of the simulation (along with the disruption of the
clump), visible in the form of a checkerboard-like pattern. We present this
peculiar artifact from our simulations in appendix B, although a detailed
investigation of it is left for future work.

##  V Field simulations

To verify our analytical understanding of kinetic relaxation and associated
nucleation of bound Bose stars, we have carried out a large suite ($\sim 500$)
of 3D simulations of the GP system (1) with varying values of the nGSI
strength $\lambda$. We evolve the GP system (1) with the following initial
Gaussian function for the ${\bm{k}}$-space Schrödinger field

$\displaystyle V^{-1/2}\,\Psi_{\bm{k}/m}\Bigl{|}_{t=0}$
$\displaystyle=e^{i\theta_{\bm{k}/m}}\sqrt{f_{\bm{v}}}\Bigl{|}_{t=0}$
$\displaystyle=e^{i\theta_{\bm{k}/m}}\left[\frac{(2\pi)^{3/2}\bar{\rho}}{m(mv_{0})^{3}}\,e^{-\frac{v^{2}}{2v_{0}^{2}}}\right]^{1/2},$
(12)

where $\theta_{\bm{k}/m}$ are random phases, uniformly distributed in
$(0,2\pi)$, for every wavenumber ${\bm{k}}$. Our numerical algorithm is based
on the well known split Fourier technique/pseudo-spectral method [64, 15, 65,
66, 67, 68], and we have used both Python based and Matlab based codes to
generate our simulation data.

As mentioned earlier, in order to be in the kinetic regime we require (i)
interactions to be tiny as compared with the typical free wave evolution
(occurring over time scales $\sim 2/mv_{0}^{2}$), and (ii) the box size to be
larger than the typical field fluctuation scale $\pi(mv_{0})^{-1}$.
Furthermore, we also impose the box size to be smaller than the gravitational
Jeans scale associated with a incoherent bound halo $\ell_{J}\sim
v_{0}(\pi/G\bar{\rho})^{1/2}$, in order to avoid its formation within our
simulation box. In this sense, our simulation box of a collection of DM waves
with typical fluctuation scale $\sim\pi(mv_{0})^{-1}$, may be regarded as a
region within a DM halo. In summary we require the following to hold true

$\displaystyle{\rm Kinetic\,regime:}\quad L$
$\displaystyle\gg\pi(mv_{0})^{-1}\quad{\&}\quad\Gamma_{\rm relax}\ll
mv_{0}^{2}/2\,,$ $\displaystyle{\rm sub\,Jeans\,scale:}\quad L$
$\displaystyle<\ell_{J}\sim v_{0}(\pi/G\bar{\rho})^{1/2}\,.$ (13)

In our simulations, we work with dimensionless quantities, for which purpose
we set $G=1/(8\pi)$ and $m=1$. More explicitly, one can rescale different
quantities in the fashion $t\rightarrow t/\mathcal{E}$,
${\bm{x}}\rightarrow{\bm{x}}/\sqrt{m\,\mathcal{E}}$,
$\psi\rightarrow\psi\,\mathcal{E}/\sqrt{8\pi Gm}$ and
$\lambda\rightarrow\lambda\mathcal{E}/(8\pi Gm^{3})$ to get Eq. (1) with both
$8\pi G$ and $m$ replaced by unity. Here $\mathcal{E}$ is any reference energy
scale in the system (for instance $\mathcal{E}=mv_{0}^{2}/2$). The
discretization in space is simply $\Delta x=L/(N_{x}-1)$ where $L$ and
$N_{x}^{3}$ are the box size and number of grid points respectively, and the
time discretization is $\Delta t=2\pi(\Delta x)^{2}m/(3\eta)$ with $\eta\geq
1$.555The $\eta\geq 1$ makes sure that there is at-least one time point in
between the full $2\pi$ rotation of the fastest oscillating mode $k_{\rm
max}\sim 2\pi/\Delta x$. Since any faithful dynamics of the system should not
be sensitive to high frequencies (corresponding to the box discretization
scale), $\eta$ can even be smaller than unity. For all our simulations, $\eta$
is at-least as big as unity. In the split Fourier technique, the field
evolution is split into a drift part where it is evolved solely due to the
gradient term (free field evolution), and a kick part where it is evolved
solely due to interactions. The Courant-Friedrichs-Lewy (CFL) condition
ensures that the fastest process in the dynamics is captured appropriately.
Hence, the fastest amongst the kick and drift processes, at any time
iteration, sets the time discretization $\Delta t$ (e.g. see [67, 68] for
details). In the kinetic regime, the time discretization is always set by the
free field evolution term $\sim(\Delta x)^{-2}/2m$, and hence by space
discretization as given above.

In all of our simulations we set $v_{0}=1/\sqrt{2}$, and choose the box size
and average mass density such that we are deep in the kinetic regime. Most of
the simulations were performed with $L=40$, $L=45$, and $L=50$ box sizes, and
the average mass densities were chosen to be small enough such that the factor
$\Gamma^{-1}_{\rm relax}\,mv_{0}^{2}/2$ was at-least as large as $\sim 250$,
going all the way up-to even $\sim 4500$. For robustness, we have performed
simulations with different grid sizes $N_{x}=\\{192,216,256\\}$, scanning over
different $\lambda/\lambda_{\rm cr}$ values. (We also performed simulations
with $N_{x}=150$ and $N_{x}=300$ to test convergence of our results (see
Appendix A).

To capture the formation of localized Bose clumps, we keep track of the mass
density in the box, radially averaged (in $\bm{k}$ space) occupation number
function $f_{k}$, the associated volume averaged correlation function
$\zeta(r)$, and the maximum mass density in the box $\rho_{\rm max}$.666In all
of our simulations, we confirm the behavior of $f_{k}$, in that it develops
increasing support towards smaller $k$ values, at-least up until nucleation.
Nucleation of a localized clump can be characterized by a change of trend of
$\rho_{\rm max}$, wherein it starts to monotonically increase beyond just the
statistical fluctuations that happen over short time scales. We record the
corresponding times in all of our simulations, both by direct inspection and
statistical methods such as moving average.777We note that this is not the
only way to know whether a bound clump has formed or not. For instance one can
alternatively construct an energy spectral function as in [35], to extract the
time scale when the function develops a support towards negative $\omega$.

In order to gauge the validity of our analytical estimate of nucleation times
(c.f. Eq (10) with Eq. (8)), and to extract the different $\alpha$
coefficients, we split the data set into two, with $\lambda=-2\pi
Gm^{2}/v_{0}^{2}$ being the splitting point. Below we elaborate on the
statistical analysis we performed in the two regimes.

Figure 2: Upper Panel: Our main figure showing the nucleation time $\tau_{\rm
nuc}$ (normalized by the gravity only case) as a function of short-ranged
self-interaction strength $\lambda$ (normalized by the critical factor $2\pi
Gm^{2}/v_{0}^{2}$). Solid gray curves are from the theory estimate Eq. (16)
(c.f. Eq. (10) with Eq. (8) and Eq. (11)), where the different $\alpha$
coefficients are obtained from least square fitting as described in the main
text. The different colored $1\sigma$ bars are from simulations (performed
with Gaussian initial conditions (V)), with box sizes $L=36$ (brown), $L=40$
(pink), $L=45$ (magenta), and $L=50$ (blue), and varying average densities
$\bar{\rho}$. Here we have only plotted one theory curve for $L=50$ (all
curves for the four different box sizes lie practically on top of each other
since the $L$ dependence is quite mild). To show the effect of the
interference term in the relaxation rate and the delay factor in the
nucleation time, we have also plotted dotted and dashed gray curves.
Respectively, these correspond to when the interference term from the
relaxation rate is set to zero, and the delay factor in the nucleation time
scale is set to unity. This delay factor is only relevant in the
$\lambda\lesssim 2\pi Gm^{2}/v_{0}^{2}$ case, and the dashed gray curve in the
left panel is simply the extension of the main solid gray curve in the right
panel. Bottom panel: Normalized $\rho_{\rm max}$ (by their respective initial
values) vs time curves for six different $\lambda$ values, highlighted by
colored points in the upper panel. The points of ‘sudden’ rise correspond to
nucleation of respective localized Bose clumps.

####  V.0.1 Net attractive interactions ($\lambda\gtrsim\lambda_{\rm cr}$).

In order to test the $\lambda$ dependence of our estimate, we construct the
quantity $r(\lambda)=\log(mv_{0}L)(\tau_{\rm nuc}(0)-\tau_{\rm
nuc}(\lambda))/2\tau_{\rm nuc}(\lambda)$ using Eq. (8) and Eq. (10). This gets
rid of the $\bar{\rho}$ and $L$ dependence, giving

$\displaystyle
r(\lambda)=\frac{\alpha_{12}}{\alpha_{1}}\left(\frac{\lambda\,v_{0}^{2}}{2\pi
Gm^{2}}\right)+\frac{\alpha_{2}}{2\alpha_{1}}\left(\frac{\lambda\,v_{0}^{2}}{2\pi
Gm^{2}}\right)^{2}\,\,.$

Not only is the curve simple enough to do statistics with, this way we can
also combine all of our simulation data (with different $\bar{\rho}$ and $L$).
The analogous quantity for simulations is

$\displaystyle\hat{r}(\lambda)=\Biggl{[}\frac{\langle\hat{\tau}_{\rm
nuc}(0)\rangle-\hat{\tau}_{\rm nuc}(\lambda)}{2\langle\hat{\tau}_{\rm
nuc}(\lambda)\rangle}\,\Biggr{]}\log(mv_{0}L)\,,$

where hats denote simulation data and angle brackets denote averaging over all
of the data (for a given $\lambda$ value). To extract the ratios
$\alpha_{12}/\alpha_{1}$ and $\alpha_{2}/\alpha_{1}$ for the theory curve, we
construct the cost function

$\displaystyle{\rm
cost}\left(\frac{\alpha_{12}}{\alpha_{1}},\frac{\alpha_{2}}{\alpha_{1}}\right)=\sum^{\sim\lambda_{\rm
cr}}_{\lambda}\frac{1}{N_{\lambda}}\sum^{N_{\lambda}}_{i=1}\Biggl{[}\frac{\hat{r}_{i}(\lambda)-r(\lambda)}{r(\lambda)}\Biggr{]}^{2}$
(14)

for least square fitting. Here $N_{\lambda}$ is the number of different
simulations performed for a given $\lambda$ value. Minimizing this cost
function then fetches the optimal values for $\alpha_{12}/\alpha_{1}$, and
$\alpha_{2}/\alpha_{1}$. For $\alpha_{1}$, we simply find the average of
$\hat{\tau}_{\rm nuc}(0)/\tau_{\rm nuc}(0)$, which we then use to get
$\alpha_{12}$ and $\alpha_{2}$ from the previous two ratios. For our Gaussian
initial condition (V), we found $\alpha_{1}\simeq 0.8$, $\alpha_{12}\simeq
1.2$, and $\alpha_{2}\simeq 1.2$.

####  V.0.2 Net repulsive interactions ($-\lambda\gtrsim-\lambda_{\rm cr}$)

In this case, we expect nucleation to happen later than relaxation, given by
Eq. (10) with the delay factor $h(\lambda)$ in Eq. (11). We can use the
previous case relationship $\tau_{\rm nuc}(\lambda\gtrsim\lambda_{\rm
cr})\simeq\Gamma^{-1}_{\rm relax}$, to test the scaling of $h(\lambda)$ for
the $-\lambda\gtrsim-\lambda_{\rm cr}$ case. From simulations, we construct
$\hat{\tau}_{\rm nuc}\Gamma_{\rm relax}$ with the three $\alpha$s in the
relaxation rate set to the ones obtained above. We then perform least square
fitting by constructing the cost function similar to the previous case

$\displaystyle{\rm cost}(\alpha_{3})=\sum^{\sim\lambda_{\rm
cr}}_{\lambda}\frac{1}{N_{\lambda}}\sum^{N_{\lambda}}_{i=1}\Biggl{[}\frac{\hat{\tau}_{\rm
nuc}\Gamma_{\rm relax}-h(\lambda)}{h(\lambda)}\Biggr{]}^{2}\,,$ (15)

and minimizing it. For Gaussian initial conditions, we found $\alpha_{3}\simeq
1$.

With the above analysis and all the four $\alpha$ values obtained, upper panel
of fig. 2 shows our main plot. We plot $\tau_{\rm nuc}(\lambda)$ (normalized
by $\tau_{\rm nuc}(0)$) as a function of $\lambda(v_{0}^{2}/2\pi Gm^{2})$:

$\displaystyle\frac{\tau_{\rm nuc}}{\tau_{\rm
nuc,0}}(x)=\frac{2\alpha_{1}\Lambda\,h(x)}{2\alpha_{1}\Lambda+4\alpha_{12}x+\alpha_{2}x^{2}}\,;\quad
x\equiv\frac{\lambda v_{0}^{2}}{2\pi Gm^{2}}\,,$ (16)

along with our simulation data. Here $\Lambda=\log(mv_{0}L)$ is the Coulomb
logarithm, and $h(x)$ is unity for $x\gtrsim-1$ (right upper panel of fig 2)
while linearly increasing as $-x\gtrsim-1$ (left upper panel of fig 2). Note
that we have only plotted one curve for $L=50$ (solid gray), since the
dependence on $L$ is very mild and renders different curves for different
values of $L$ practically on top of each other. The error bars correspond to
$1$-$\sigma$ fluctuations (owing to random and different initial condition for
every simulation seed), with different colors corresponding to three different
box sizes considered. In general, the agreement between analytical estimates
and simulations is evident. Let us highlight our two main results: (a) The
rising feature as $\lambda$ goes from positive to negative, with a peak
occurring around $\lambda\simeq-2\pi Gm^{2}/v_{0}^{2}\simeq\lambda_{\rm cr}$,
is a clear evidence of the interference term in the relaxation rate. To
represent the effect of the interference term visually, we have also plotted a
dotted gray curve (in the upper right panel of fig. 2), which is equal to
inverse of the relaxation rate with the interference term dropped. That is,
inverse of Eq. (8) with the term $\propto G\lambda$ set to zero; (b) To the
left of the peak and increasing $-\lambda$, nucleation happens later than just
the inverse relaxation rate. The delay factor $h$ and the relaxation rate
$\Gamma_{\rm relax}$ scale as $\sim-\lambda$ and $\lambda^{2}$ (to leading
order) respectively, resulting in the scaling of the nucleation time as
$(-\lambda)^{-1}$ (and not $\lambda^{-2}$) to leading order. To highlight
this, we have augmented the upper left panel of fig. 2 with just the
relaxation time curve, i.e. $\Gamma_{\rm rel}^{-1}$, shown in dashed gray.
(This is nothing but the solid gray curve on the right upper panel, extended
towards the left upper panel).

Figure 3: Density projection snapshots for 6 different $\lambda$ values, at
different times $\hat{t}$ in the respective simulations. Upper Panel:
Snapshots for three $\lambda$ values in the typical net attractive regime
$\lambda\gtrsim\lambda_{\rm cr}\simeq-2\pi Gm^{2}/v_{0}^{2}$. In the rightmost
snapshot, for $\lambda\approx 3\,\lambda_{\rm cr}$, the nucleated Bose clump
quickly collapses (within $5dt$), shown in the smaller right corner image.
Bottom panel: Snapshots for three different $\lambda$ values for the other
case of typical net repulsive self-interactions $-\lambda\gtrsim-\lambda_{\rm
cr}\simeq 2\pi Gm^{2}/v_{0}^{2}$.

The lower panel of fig. 2 shows moving averaged $\rho_{\rm max}$ vs time
curves for six simulations with different $\lambda$ values. The unambiguous
“sudden” rise in $\rho_{\rm max}$ marks the nucleation of a localized object
within which density grows over time.888In all of our simulations, we have
explicitly verified, by visually tracking the simulation box, that this rising
feature indeed corresponds to appearance of an overdense region.

As visual examples, in fig. 3 we also present density projection snapshots for
six different $\lambda$ values at later times, showing the presence of
nucleated Bose stars. For attractive short-ranged self-interaction
$\lambda>0$, nucleated Bose clumps eventually collapse into a Bosenova. This
happens when it reaches the critical mass where it can no longer remain stable
(see section IV.2).

##  VI Summary and Discussion

In this paper we have investigated kinetic relaxation and associated
nucleation times of Bose stars, in scalar fuzzy dark matter with short-ranged
$2$-body self-interactions. Starting with the wave-kinetic Boltzmann equation
for the mode occupation number function (which we derived in an earlier work),
we first highlighted the presence of a cross/interference term $\propto
G\lambda$ in the rate of relaxation $\Gamma_{\rm relax}$, alongside the usual
two terms $\propto G^{2}$ and $\lambda^{2}$ due to both gravitational and
short-ranged self-interaction individually. This is because of the wave-
mechanical nature of the system: The rate depends on the total cross section,
which is not just the sum of individual cross sections due to the different
processes. Rather the scattering amplitudes due to all the processes must be
added first, and then use its absolute square to get the cross section and
associated rate of relaxation/condensation.

The presence of this cross term gives rise to a critical repulsive self-
interaction strength $\lambda_{\rm cr}\approx-2\pi Gm^{2}/v_{0}^{2}$, around
which the typical net self-interaction (due to both gravitational and short-
ranged self-interaction), transitions between being attractive and repulsive,
and the relaxation rate becomes smallest. Here $k_{0}=mv_{0}$ is the typical
wave-mode present in the system initially.

For nucleation times as a function of $\lambda$, we found that for the case of
net attractive self-interaction $\lambda\gtrsim\lambda_{\rm cr}$, nucleation
happens quickly upon relaxation, giving rise to the relationship $\tau_{\rm
nuc}\simeq\Gamma^{-1}_{\rm relax}$. One the other hand for net repulsive self-
interaction $-\lambda\gtrsim-\lambda_{\rm cr}$, nucleation is delayed. This is
because upon relaxation, short-ranged self-interaction dominate over
gravitational self-interaction, preventing nucleation of a bound object. Over
time as more particles are driven towards the condensate phase (equivalently,
as field correlation length scale increases along with diminishing density
fluctuations), a potential arises for the formation of a bound region where
gravity can now overcome both the wave pressure and short-ranged self-
interaction. The associated delay factor rises linearly with $-\lambda$,
giving the nucleation time scale as $\tau_{\rm nuc}\simeq(\lambda/\lambda_{\rm
cr})\Gamma^{-1}_{\rm relax}$. In summary, Eq. (10) along with Eq. (8) (with
the delay factor give in Eq. (11)) is our main analytical estimate for the
nucleation timescale of Bose stars, as a function of the short-ranged self-
interaction strength $\lambda$.

To analyze this, we performed a large suit of $3$+$1$ dimensional simulations
of the Schrödinger -Poisson / Gross-Pitaevskii system (Eq. (1)), for many
different values of $\lambda$ and different parameters such as the box size
$L$ and average mass density $\bar{\rho}$. All of our simulations were carried
out with Maxwell-Boltzmann distribution, with random phases for each value of
the wavemode ${\bm{k}}$ (Eq. (V)). Throughout most of our simulations, we kept
track of the max density in the box, occupation number function $f_{k}$
(radially averaged $f_{\bm{k}}$ in ${\bm{k}}$ space), the associated
correlation function $\zeta(r)$, and projected mass density along some
direction. By reading the times at which $\rho_{\rm max}$ starts to
monotonically rise beyond just the statistical fluctuations (together with
making sure that a localized over dense region does appear in the simulation
box around this time), we record the times of nucleation. Upper panel of Fig.
2 presents the comparison between simulations and analytical estimate. As
examples, the figure is also appended (lower panel) with $\rho_{\rm max}$ vs
time curves for six different $\lambda$ values.

While in this paper we have not analyzed our simulation data for the rate at
which Bose stars accrete mass from their surroundings, we kept track of the
eventual behavior of these objects (post nucleation), for many of our
simulations. For the attractive case ($\lambda>0$), we confirmed that the
nucleated Bose stars eventually decay away. This is expected since there
exists a maximum critical mass beyond which the star becomes unstable and
collapses into a Bosenova. For instance see upper right snapshot in Fig. 3,
when the nucleated star ‘immediately’ collapses. While the study of eventual
dynamics and fate of such regions requires a full relativistic treatment,
field dynamics up-to this point is well described by the non-relativistic GP
equation (e.g. see [51]).

For the repulsive case $\lambda<0$, we found a peculiar decay behavior. We
find that the nucleated clump eventually (over time scales longer than the
nucleation time) reaches a type of criticality at which point very high
frequency modes, passing through the clump and travelling along the three
directions of the simulation box, appear in the system. See appendix B for
some discussion. This could be an artifact of the periodic boundary conditions
of the split-Fourier simulation setup, and if so, bringing into question its
use to study long term dynamics of fuzzy dark matter with repulsive short-
ranged self-interactions via such simulation setups. We leave a detailed
investigation of this behavior for a separate work.

###  VI.1 Comparison with earlier work

Let us now compare our results with some of the earlier work on the subject of
kinetic nucleation of Bose stars. First, our results encompass the result of
[35] for the gravity only ($\lambda=0$) case, and is even in very good
agreement with the order unity coefficient $\alpha_{1}$ in the rate expression
(besides the overall scaling with $\bar{\rho}$, $m$, $v_{0}$, $L$ and $G$),
obtained for Gaussian initial condition. Upon inclusion of short-ranged self-
interaction ($\lambda\neq 0$ case), our results differ significantly from the
existing literature [37, 38, 39]. First, we find that there exists an
interference term $\propto G\lambda$ in the relaxation rate, which in fact is
the leading order $\lambda$ dependent term when short-ranged self-interaction
is not dominating over gravitational self-interaction. Only in the scenario
when the former is dominant, does the relaxation rate goes as $\lambda^{2}$ to
leading order. Secondly, the nucleation time scale is not always equal to the
inverse relaxation rate. While for $\lambda\gtrsim\lambda_{\rm cr}$,
nucleation time scale is just the inverse relaxation rate, for the strong
repulsive self-interaction $-\lambda\gtrsim-\lambda_{\rm cr}$, nucleation time
is delayed by an extra factor of $(\lambda/\lambda_{\rm cr})$. Therefore for
the purposes of nucleation of Bose stars, only in the case of strong
attractive short-ranged self-interaction, $\lambda\gtrsim-\lambda_{\rm cr}$,
is it true that the nucleation time goes as $\lambda^{-2}$ to leading order.
For in the opposite case of strong repulsive short-ranged self-interaction,
$-\lambda\gtrsim-\lambda_{\rm cr}$, the nucleation time scale goes as
$\lambda^{-1}$ to leading order instead.

###  VI.2 Implications

Our results could have important implications in the context of self-
interacting fuzzy dark matter and various interesting phenomenon that it
entails. The appearance of the interference term in the relaxation rate, and
hence in the nucleation time scale of Boson stars, may modify results for some
of the phenomenon such as recurrent axinovae [52], de-stabilization of
gravitational atoms [69], among others.

In general, irrespective of the nature (attractive or repulsive) of point-like
self-interaction, the interference term becomes the leading order $\lambda$
dependent term (and hence extremely important), when $|\lambda|$ is at best
comparable to the critical value $|\lambda_{\rm cr}|$. As an example, even for
the QCD axion we have $\lambda_{\rm qcd}/|\lambda_{\rm cr}|\simeq
1.3(v_{0}^{2}m_{\mathrm{pl}}^{2}/f_{a}^{2})$, hence becoming comparable to, or
less than $|\lambda_{\rm cr}|$, in cosmological environments with
$v_{0}\lesssim(f_{a}/m_{\mathrm{pl}})\sim 10^{-5}$. For instance this could be
important in the study of axion miniclusters [70].

In this paper we have focused on kinetic nucleation via both gravitational and
short-ranged self-interactions for a single scalar field. A natural
generalization is to include multiple scalar fields with naturally different
masses and 4-point interactions, or a single spin-$1$ field including density-
density and spin-spin interactions [71, 57], or even multiple spin-$1$ fields
with extra Yang-Mills interactions [57], or a combination thereof. While there
would necessarily be interference terms $\propto G\lambda$, and we expect
similar scaling of nucleation time scales as presented in this work (as a
function of $\lambda$), a detailed analysis of such cases is left for future
work.999We thank Benjamin Schussler for carrying out some preliminary
simulations for the self-interacting vector case, confirming the presence of
the interference term in the relaxation rate.

## Acknowledgements

We thank Mustafa Amin, Mark Hertzberg, Andrew Long, and David J.E. Marsh for
many helpful discussions and also their comments on this manuscript. MJ is
partially supported by a DOE grant DE-SC0021619, and partly supported by a
Leverhulme Trust Research Project (RPG-2022-145). JT and WW acknowledge
undergraduate summer support from the Department of Physics and Astronomy at
Rice University.

## References

  * Arvanitaki  _et al._ [2010] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, String Axiverse, Phys. Rev. D 81, 123530 (2010), arXiv:0905.4720 [hep-th] .
  * Acharya  _et al._ [2010] B. S. Acharya, K. Bobkov, and P. Kumar, An M Theory Solution to the Strong CP Problem and Constraints on the Axiverse, JHEP 11, 105, arXiv:1004.5138 [hep-th] .
  * Cicoli  _et al._ [2012] M. Cicoli, M. Goodsell, and A. Ringwald, The type IIB string axiverse and its low-energy phenomenology, JHEP 10, 146, arXiv:1206.0819 [hep-th] .
  * Giedt [2001] J. Giedt, Completion of standard model like embeddings, Annals Phys. 289, 251 (2001), arXiv:hep-th/0009104 .
  * Taylor and Wang [2016] W. Taylor and Y.-N. Wang, A Monte Carlo exploration of threefold base geometries for 4d F-theory vacua, JHEP 01, 137, arXiv:1510.04978 [hep-th] .
  * Alexander  _et al._ [2023] S. Alexander, H. Gilmer, T. Manton, and E. McDonough, The $\pi$-axion and $\pi$-axiverse of dark QCD (2023), arXiv:2304.11176 [hep-ph] .
  * Peccei and Quinn [1977] R. D. Peccei and H. R. Quinn, CP Conservation in the Presence of Instantons, Phys. Rev. Lett. 38, 1440 (1977).
  * Weinberg [1978] S. Weinberg, A New Light Boson?, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek [1978] F. Wilczek, Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons, Phys. Rev. Lett. 40, 279 (1978).
  * Preskill  _et al._ [1983] J. Preskill, M. B. Wise, and F. Wilczek, Cosmology of the Invisible Axion, Phys. Lett. B 120, 127 (1983).
  * Abbott and Sikivie [1983] L. Abbott and P. Sikivie, A Cosmological Bound on the Invisible Axion, Phys. Lett. B 120, 133 (1983).
  * Dine and Fischler [1983] M. Dine and W. Fischler, The Not So Harmless Axion, Phys. Lett. B 120, 137 (1983).
  * Matos and Urena-Lopez [2001] T. Matos and L. A. Urena-Lopez, A Further analysis of a cosmological model of quintessence and scalar dark matter, Phys. Rev. D 63, 063506 (2001), arXiv:astro-ph/0006024 .
  * Hu  _et al._ [2000] W. Hu, R. Barkana, and A. Gruzinov, Cold and fuzzy dark matter, Phys. Rev. Lett. 85, 1158 (2000), arXiv:astro-ph/0003365 .
  * Mocz  _et al._ [2017] P. Mocz, M. Vogelsberger, V. H. Robles, J. Zavala, M. Boylan-Kolchin, A. Fialkov, and L. Hernquist, Galaxy formation with BECDM – I. Turbulence and relaxation of idealized haloes, Mon. Not. Roy. Astron. Soc. 471, 4559 (2017), arXiv:1705.05845 [astro-ph.CO] .
  * Brito  _et al._ [2015] R. Brito, V. Cardoso, and P. Pani, Superradiance: New Frontiers in Black Hole Physics, Lect. Notes Phys. 906, pp.1 (2015), arXiv:1501.06570 [gr-qc] .
  * Schive  _et al._ [2014] H.-Y. Schive, T. Chiueh, and T. Broadhurst, Cosmic Structure as the Quantum Interference of a Coherent Dark Wave, Nature Phys. 10, 496 (2014), arXiv:1406.6586 [astro-ph.GA] .
  * Ruffini and Bonazzola [1969] R. Ruffini and S. Bonazzola, Systems of selfgravitating particles in general relativity and the concept of an equation of state, Phys. Rev. 187, 1767 (1969).
  * Kolb and Tkachev [1994] E. W. Kolb and I. I. Tkachev, Nonlinear axion dynamics and formation of cosmological pseudosolitons, Phys. Rev. D 49, 5040 (1994), arXiv:astro-ph/9311037 .
  * Chavanis [2011] P.-H. Chavanis, Mass-radius relation of Newtonian self-gravitating Bose-Einstein condensates with short-range interactions: I. Analytical results, Phys. Rev. D 84, 043531 (2011), arXiv:1103.2050 [astro-ph.CO] .
  * Amin  _et al._ [2012] M. A. Amin, R. Easther, H. Finkel, R. Flauger, and M. P. Hertzberg, Oscillons After Inflation, Phys. Rev. Lett. 108, 241302 (2012), arXiv:1106.3335 [astro-ph.CO] .
  * Brito  _et al._ [2016] R. Brito, V. Cardoso, C. A. R. Herdeiro, and E. Radu, Proca stars: Gravitating Bose–Einstein condensates of massive spin 1 particles, Phys. Lett. B 752, 291 (2016), arXiv:1508.05395 [gr-qc] .
  * Aoki  _et al._ [2018] K. Aoki, K.-i. Maeda, Y. Misonoh, and H. Okawa, Massive Graviton Geons, Phys. Rev. D 97, 044005 (2018), arXiv:1710.05606 [gr-qc] .
  * Amin and Mocz [2019] M. A. Amin and P. Mocz, Formation, gravitational clustering, and interactions of nonrelativistic solitons in an expanding universe, Phys. Rev. D 100, 063507 (2019), arXiv:1902.07261 [astro-ph.CO] .
  * Adshead and Lozanov [2021] P. Adshead and K. D. Lozanov, Self-gravitating Vector Dark Matter, Phys. Rev. D 103, 103501 (2021), arXiv:2101.07265 [gr-qc] .
  * Jain and Amin [2022a] M. Jain and M. A. Amin, Polarized solitons in higher-spin wave dark matter, Phys. Rev. D 105, 056019 (2022a), arXiv:2109.04892 [hep-th] .
  * Gorghetto  _et al._ [2022] M. Gorghetto, E. Hardy, J. March-Russell, N. Song, and S. M. West, Dark Photon Stars: Formation and Role as Dark Matter Substructure (2022), arXiv:2203.10100 [hep-ph] .
  * Amin  _et al._ [2022] M. A. Amin, M. Jain, R. Karur, and P. Mocz, Small-scale structure in vector dark matter (2022), arXiv:2203.11935 [astro-ph.CO] .
  * Hui  _et al._ [2021] L. Hui, A. Joyce, M. J. Landry, and X. Li, Vortices and waves in light dark matter, JCAP 01, 011, arXiv:2004.01188 [astro-ph.CO] .
  * Gosenca  _et al._ [2023] M. Gosenca, A. Eberhardt, Y. Wang, B. Eggemeier, E. Kendall, J. L. Zagorac, and R. Easther, Multifield ultralight dark matter, Phys. Rev. D 107, 083014 (2023), arXiv:2301.07114 [astro-ph.CO] .
  * Guth  _et al._ [2015] A. H. Guth, M. P. Hertzberg, and C. Prescod-Weinstein, Do Dark Matter Axions Form a Condensate with Long-Range Correlation?, Phys. Rev. D 92, 103513 (2015), arXiv:1412.5930 [astro-ph.CO] .
  * Ferreira [2021] E. G. M. Ferreira, Ultra-light dark matter, Astron. Astrophys. Rev. 29, 7 (2021), arXiv:2005.03254 [astro-ph.CO] .
  * Hui [2021] L. Hui, Wave Dark Matter, Ann. Rev. Astron. Astrophys. 59, 247 (2021), arXiv:2101.11735 [astro-ph.CO] .
  * Semikoz and Tkachev [1995] D. V. Semikoz and I. I. Tkachev, Kinetics of Bose condensation, Phys. Rev. Lett. 74, 3093 (1995), arXiv:hep-ph/9409202 .
  * Levkov  _et al._ [2018] D. G. Levkov, A. G. Panin, and I. I. Tkachev, Gravitational Bose-Einstein condensation in the kinetic regime, Phys. Rev. Lett. 121, 151301 (2018), arXiv:1804.05857 [astro-ph.CO] .
  * Eggemeier and Niemeyer [2019] B. Eggemeier and J. C. Niemeyer, Formation and mass growth of axion stars in axion miniclusters, Phys. Rev. D 100, 063528 (2019), arXiv:1906.01348 [astro-ph.CO] .
  * Kirkpatrick  _et al._ [2020] K. Kirkpatrick, A. E. Mirasola, and C. Prescod-Weinstein, Relaxation times for Bose-Einstein condensation in axion miniclusters, Phys. Rev. D 102, 103012 (2020), arXiv:2007.07438 [hep-ph] .
  * Chen  _et al._ [2022] J. Chen, X. Du, E. W. Lentz, and D. J. E. Marsh, Relaxation times for Bose-Einstein condensation by self-interaction and gravity, Phys. Rev. D 106, 023009 (2022), arXiv:2109.11474 [astro-ph.CO] .
  * Kirkpatrick  _et al._ [2022] K. Kirkpatrick, A. E. Mirasola, and C. Prescod-Weinstein, Analysis of Bose-Einstein condensation times for self-interacting scalar dark matter, Phys. Rev. D 106, 043512 (2022), arXiv:2110.08921 [hep-ph] .
  * Purohit  _et al._ [2021] K. J. Purohit, P. K. Natwariya, J. R. Bhatt, and P. K. Mehta, Rotating Scalar Field and Formation of Bose Stars (2021), arXiv:2109.02601 [astro-ph.CO] .
  * Jain  _et al._ [2023a] M. Jain, M. A. Amin, J. Thomas, and W. Wanichwecharungruang, Kinetic relaxation and Bose-star formation in multicomponent dark matter- I (2023a), arXiv:2304.01985 [astro-ph.CO] .
  * Chen  _et al._ [2023] J. Chen, X. Du, M. Zhou, A. Benson, and D. J. E. Marsh, Gravitational Bose-Einstein Condensation of Vector/Hidden Photon Dark Matter (2023), arXiv:2304.01965 [astro-ph.HE] .
  * Chen  _et al._ [2021] J. Chen, X. Du, E. W. Lentz, D. J. E. Marsh, and J. C. Niemeyer, New insights into the formation and growth of boson stars in dark matter halos, Phys. Rev. D 104, 083022 (2021), arXiv:2011.01333 [astro-ph.CO] .
  * Chan  _et al._ [2022] J. H.-H. Chan, S. Sibiryakov, and W. Xue, Condensation and Evaporation of Boson Stars (2022), arXiv:2207.04057 [astro-ph.CO] .
  * Dmitriev  _et al._ [2023] A. S. Dmitriev, D. G. Levkov, A. G. Panin, and I. I. Tkachev, Self-similar growth of Bose stars (2023), arXiv:2305.01005 [astro-ph.CO] .
  * Zakharov  _et al._ [1992] V. E. Zakharov, V. S. L’vov, and G. Falkovich, Kolmogorov spectra of turbulence i (1992).
  * Chavanis [2021] P.-H. Chavanis, Landau equation for self-gravitating classical and quantum particles: application to dark matter, Eur. Phys. J. Plus 136, 703 (2021), arXiv:2012.12858 [astro-ph.GA] .
  * Chavanis and Delfini [2011] P. H. Chavanis and L. Delfini, Mass-radius relation of Newtonian self-gravitating Bose-Einstein condensates with short-range interactions: II. Numerical results, Phys. Rev. D 84, 043532 (2011), arXiv:1103.2054 [astro-ph.CO] .
  * Chavanis [2022] P.-H. Chavanis, Maximum mass of relativistic self-gravitating Bose-Einstein condensates with repulsive or attractive $|\varphi|^{4}$ self-interaction (2022), arXiv:2211.13237 [gr-qc] .
  * Schiappacasse and Hertzberg [2018] E. D. Schiappacasse and M. P. Hertzberg, Analysis of Dark Matter Axion Clumps with Spherical Symmetry, JCAP 01, 037, [Erratum: JCAP 03, E01 (2018)], arXiv:1710.04729 [hep-ph] .
  * Levkov  _et al._ [2017] D. G. Levkov, A. G. Panin, and I. I. Tkachev, Relativistic axions from collapsing Bose stars, Phys. Rev. Lett. 118, 011301 (2017), arXiv:1609.03611 [astro-ph.CO] .
  * Fox  _et al._ [2023] P. J. Fox, N. Weiner, and H. Xiao, Recurrent Axinovae and their Cosmological Constraints (2023), arXiv:2302.00685 [hep-ph] .
  * Boehmer and Harko [2007] C. G. Boehmer and T. Harko, Can dark matter be a Bose-Einstein condensate?, JCAP 06, 025, arXiv:0705.4158 [astro-ph] .
  * Magana and Matos [2012] J. Magana and T. Matos, A brief Review of the Scalar Field Dark Matter model, J. Phys. Conf. Ser. 378, 012012 (2012), arXiv:1201.6107 [astro-ph.CO] .
  * Salehian  _et al._ [2021] B. Salehian, H.-Y. Zhang, M. A. Amin, D. I. Kaiser, and M. H. Namjoo, Beyond Schrödinger-Poisson: nonrelativistic effective field theory for scalar dark matter, JHEP 09, 050, arXiv:2104.10128 [astro-ph.CO] .
  * Croon  _et al._ [2019] D. Croon, J. Fan, and C. Sun, Boson Star from Repulsive Light Scalars and Gravitational Waves, JCAP 04, 008, arXiv:1810.01420 [hep-ph] .
  * Jain [2022] M. Jain, Soliton stars in Yang-Mills-Higgs theories, Phys. Rev. D 106, 085011 (2022), arXiv:2205.03418 [hep-ph] .
  * Hertzberg  _et al._ [2021] M. P. Hertzberg, F. Rompineve, and J. Yang, Decay of Boson Stars with Application to Glueballs and Other Real Scalars, Phys. Rev. D 103, 023536 (2021), arXiv:2010.07927 [hep-ph] .
  * Dawoodbhoy  _et al._ [2021] T. Dawoodbhoy, P. R. Shapiro, and T. Rindler-Daller, Core-envelope haloes in scalar field dark matter with repulsive self-interaction: fluid dynamics beyond the de Broglie wavelength, Mon. Not. Roy. Astron. Soc. 506, 2418 (2021), arXiv:2104.07043 [astro-ph.CO] .
  * Shapiro  _et al._ [2021] P. R. Shapiro, T. Dawoodbhoy, and T. Rindler-Daller, Cosmological structure formation in scalar field dark matter with repulsive self-interaction: the incredible shrinking Jeans mass, Mon. Not. Roy. Astron. Soc. 509, 145 (2021), arXiv:2106.13244 [astro-ph.CO] .
  * Hartman  _et al._ [2022a] S. T. H. Hartman, H. A. Winther, and D. F. Mota, Constraints on self-interacting Bose-Einstein condensate dark matter using large-scale observables, JCAP 02 (02), 005, arXiv:2108.07496 [astro-ph.CO] .
  * Hartman  _et al._ [2022b] S. T. H. Hartman, H. A. Winther, and D. F. Mota, Cosmological simulations of self-interacting Bose-Einstein condensate dark matter, Astron. Astrophys. 666, A95 (2022b), arXiv:2203.03946 [astro-ph.CO] .
  * Foidl  _et al._ [2023] H. Foidl, T. Rindler-Daller, and W. W. Zeilinger, Halo formation and evolution in scalar field dark matter and cold dark matter: New insights from the fluid approach, Phys. Rev. D 108, 043012 (2023), arXiv:2305.12982 [astro-ph.GA] .
  * Springel [2005] V. Springel, The Cosmological simulation code GADGET-2, Mon. Not. Roy. Astron. Soc. 364, 1105 (2005), arXiv:astro-ph/0505010 .
  * Edwards  _et al._ [2018] F. Edwards, E. Kendall, S. Hotchkiss, and R. Easther, PyUltraLight: A Pseudo-Spectral Solver for Ultralight Dark Matter Dynamics, JCAP 10, 027, arXiv:1807.04037 [astro-ph.CO] .
  * Glennon and Prescod-Weinstein [2021] N. Glennon and C. Prescod-Weinstein, Modifying PyUltraLight to model scalar dark matter with self-interactions, Phys. Rev. D 104, 083532 (2021), arXiv:2011.09510 [astro-ph.CO] .
  * Jain and Amin [2022b] M. Jain and M. A. Amin, i-SPin: An integrator for multicomponent Schrödinger-Poisson systems with self-interactions (2022b), arXiv:2211.08433 [astro-ph.CO] .
  * Jain  _et al._ [2023b] M. Jain, M. A. Amin, and H. Pu, i-SPin 2: An integrator for general spin-s Gross-Pitaevskii systems (2023b), arXiv:2305.01675 [cond-mat.quant-gas] .
  * Budker  _et al._ [2023] D. Budker, J. Eby, M. Gorghetto, M. Jiang, and G. Perez, A Generic Formation Mechanism of Ultralight Dark Matter Solar Halos (2023), arXiv:2306.12477 [hep-ph] .
  * Kolb and Tkachev [1993] E. W. Kolb and I. I. Tkachev, Axion miniclusters and Bose stars, Phys. Rev. Lett. 71, 3051 (1993), arXiv:hep-ph/9303313 .
  * Zhang  _et al._ [2021] H.-Y. Zhang, M. Jain, and M. A. Amin, Polarized Vector Oscillons (2021), arXiv:2111.08700 [astro-ph.CO] .

##  Appendix A Statistical convergence

Here we show convergence and reliability of our simulation results. Figure 4
compares the main theory curve (in solid gray), with data points for the lower
($N_{x}=150$) and higher ($N_{x}=300$) resolution grids, compared to
$N_{x}=192$, $N_{x}=216$ and $N_{x}=256$ used for our results presented in the
main text. See caption for details. Convergence of our results is evident from
this plot. Note that this is not the usual convergence, where two or more
simulations with similar initial conditions are performed, with different
values of $\Delta x$ and $\Delta t$. Rather, here we show a ‘statistical
convergence’ of sorts. Also, since $\Delta t\propto(\Delta x)^{2}$, increasing
grid size reduces both $\Delta x$ and $\Delta t$.

Figure 4: Similar to the upper panel of figure 2, for simulations performed
with lower and higher resolution grids compared to the ones used in the main
text. With box size $L=40$, the $1\sigma$ bars represent simulations performed
using $150^{3}$ grid, while the solid points are from simulations using
$300^{3}$ grid. The solid gray curves correspond to the analytical estimate
Eq. (16) (c.f. Eq. (10) with Eq. (8) and Eq. (11)), with the $\mathcal{O}(1)$
$\alpha$ coefficients obtained using simulations performed with grid sizes
$192^{3}$, $216^{3}$, and $256^{3}$ as in the main text.

##  Appendix B Peculiar appearance of high frequency modes for the repulsive
case

In the case of repulsive short-ranged self interaction $\lambda<0$, we find a
peculiar behavior over long time scales (later than nucleation). There starts
to appear high frequency modes (near cutoff of the simulation box) that pass
through the clump and in all the three perpendicular directions of the box.
This leads to disruption of the simulation (manifesting in the form of
checkerboard like pattern) and of the clump, as these waves circulate within
the periodic simulation box. An artifact of this is wiping out of density
fluctuations and eventual homogenization of the simulation box.

In fig. 5 we provide some simulation snapshots of this peculiarity, for three
different values of $\lambda$. Notice the appearance of checkerboard like
pattern in the right hand side panel.

Figure 5: Density projection snapshots for 3 values or repulsive short-ranged
self-interaction strength $\lambda$, at two different times in the respective
simulations. The grid size is $216^{3}$, and box length is $50$. The left hand
side snapshots are the same as presented in the main text in fig. 3. The right
hand side snapshots are at later times, when the simulation has now been
rendered ‘unfaithful’. Notice the checkerboard like pattern in all of the
snapshots in the right panel.

In order to check if the phenomenon is an artifact of finite discretization,
we simulated lower and higher resolution grids with the same initial
conditions for long times. We found no clues if this is the case or not. As an
example, for $\lambda\approx-4.2\lambda_{\rm cr}$, with $L=40$ and
$\bar{\rho}=0.01$, in the lower resolution ($128^{3}$) and higher resolution
grids ($256^{3}$) the nucleation times were $\sim 3600$ and $\sim 3800$
respectively. This shows decent convergence of the nucleation time. At the
same time however, the onset of these high frequency waves for both the grids
were also similar ($\sim 5300$ and $\sim 5700$ respectively).

We also performed two other tests to see if something can be learned about
this phenomenon. In one test we put a pre-computed soliton in a bath of DM
waves. We observed the same appearance of high frequency modes and
checkerboard pattern developing, leading to disruption of the soliton. The
times at which this happens depended upon the average mass density and total
mass of the soliton. In yet another test, we performed a few runs with larger
average mass densities, such that the associated gravitational Jeans scale
$\ell_{J}\sim v_{0}(\pi/G\bar{\rho})^{1/2}$ became smaller than the box size.
Even in this case, upon formation of a halo, the same peculiar phenomenon
appears.

While we have confirmed the presence of this feature in our (split-Fourier
technique based-)simulations through multiple tests, we leave a detailed
analysis of this peculiarity for future work.


|  |  |  |  15NH3(1,1) | 56 |  |  |  |  |  15NH3(1,1) | 66
|  |  |  |  15NH3(2,2) | 54 |  |  |  |  |  15NH3(2,2) | 98
|  |  |  |  15NH3(3,3) | 54 |  |  |  |  |  15NH3(3,3) | 70
G183.72 | TMRT | 05:40:24.23 +23:50:54.7 | 39 |  14NH3(1,1) | 68 | J024929.80 | TMRT | 02:49:29.80 60:47:29.50 | 8 |  14NH3(1,1) | 79
|  |  |  |  14NH3(2,2) | 88 |  |  |  |  |  14NH3(2,2) | 98
|  |  |  |  14NH3(3,3) | 80 |  |  |  |  |  14NH3(3,3) | 83
|  |  |  |  15NH3(1,1) | 104 |  |  |  |  |  15NH3(1,1) | 82
|  |  |  |  15NH3(2,2) | 99 |  |  |  |  |  15NH3(2,2) | 85
|  |  |  |  15NH3(3,3) | 100 |  |  |  |  |  15NH3(3,3) | 77
G188.79 | TMRT | 06:09:06.97 +21:50:41.4 | 30 |  14NH3(1,1) | 63 | J040136.49 | TMRT | 04:01:36.49 53:19:41.40 | 9 |  14NH3(1,1) | 60
|  |  |  |  14NH3(2,2) | 63 |  |  |  |  |  14NH3(2,2) | 73
|  |  |  |  14NH3(3,3) | 62 |  |  |  |  |  14NH3(3,3) | 66
|  |  |  |  15NH3(1,1) | 85 |  |  |  |  |  15NH3(1,1) | 81
|  |  |  |  15NH3(2,2) | 75 |  |  |  |  |  15NH3(2,2) | 77
|  |  |  |  15NH3(3,3) | 74 |  |  |  |  |  15NH3(3,3) | 69
G188.79 | TMRT | 06:09:06.96 21:50:41.30 | 9 |  14NH3(1,1) | 76 | J051713.30 | TMRT | 05:17:13.30 39:22:14.00 | 9 |  14NH3(1,1) | 85
|  |  |  |  14NH3(2,2) | 65 |  |  |  |  |  14NH3(2,2) | 85
|  |  |  |  14NH3(3,3) | 88 |  |  |  |  |  14NH3(3,3) | 100
|  |  |  |  15NH3(1,1) | 82 |  |  |  |  |  15NH3(1,1) | 73
|  |  |  |  15NH3(2,2) | 81 |  |  |  |  |  15NH3(2,2) | 98
|  |  |  |  15NH3(3,3) | 78 |  |  |  |  |  15NH3(3,3) | 98
G188.94 | TMRT | 06:08:53.3 +21:38:28.7 | 176 |  14NH3(1,1) | 84 | J053927.60 | TMRT | 05:39:27.60 35:30:58.00 | 9 |  14NH3(1,1) | 68
|  |  |  |  14NH3(2,2) | 12 |  |  |  |  |  14NH3(2,2) | 88
|  |  |  |  14NH3(3,3) | 11 |  |  |  |  |  14NH3(3,3) | 69
|  |  |  |  15NH3(1,1) | 97 |  |  |  |  |  15NH3(1,1) | 62
|  |  |  |  15NH3(2,2) | 12 |  |  |  |  |  15NH3(2,2) | 91
|  |  |  |  15NH3(3,3) | 12 |  |  |  |  |  15NH3(3,3) | 90
G188.94 | TMRT | 06:08:53.32 21:38:29.10 | 9 |  14NH3(1,1) | 75 | J054137.40 | TMRT | 05:41:37.40 35:48:49.00 | 9 |  14NH3(1,1) | 70
|  |  |  |  14NH3(2,2) | 76 |  |  |  |  |  14NH3(2,2) | 83
|  |  |  |  14NH3(3,3) | 74 |  |  |  |  |  14NH3(3,3) | 97
|  |  |  |  15NH3(1,1) | 67 |  |  |  |  |  15NH3(1,1) | 80
|  |  |  |  15NH3(2,2) | 75 |  |  |  |  |  15NH3(2,2) | 99
|  |  |  |  15NH3(3,3) | 67 |  |  |  |  |  15NH3(3,3) | 70
G189.03 | TMRT | 06:08:40.65 21:31:07.00 | 9 |  14NH3(1,1) | 64 | J213213.00 | TMRT | 21:32:13.00 55:52:56.00 | 9 |  14NH3(1,1) | 88
|  |  |  |  14NH3(2,2) | 65 |  |  |  |  |  14NH3(2,2) | 85
|  |  |  |  14NH3(3,3) | 64 |  |  |  |  |  14NH3(3,3) | 84
|  |  |  |  15NH3(1,1) | 70 |  |  |  |  |  15NH3(1,1) | 90
|  |  |  |  15NH3(2,2) | 70 |  |  |  |  |  15NH3(2,2) | 77
|  |  |  |  15NH3(3,3) | 66 |  |  |  |  |  15NH3(3,3) | 73
G189.47 | TMRT | 06:02:08.37 20:09:20.10 | 9 |  14NH3(1,1) | 56 | J214301.50 | TMRT | 21:43:01.50 54:56:18.50 | 9 |  14NH3(1,1) | 67
|  |  |  |  14NH3(2,2) | 70 |  |  |  |  |  14NH3(2,2) | 73
|  |  |  |  14NH3(3,3) | 66 |  |  |  |  |  14NH3(3,3) | 88
|  |  |  |  15NH3(1,1) | 74 |  |  |  |  |  15NH3(1,1) | 73
|  |  |  |  15NH3(2,2) | 65 |  |  |  |  |  15NH3(2,2) | 77
|  |  |  |  15NH3(3,3) | 73 |  |  |  |  |  15NH3(3,3) | 71
G189.77 | TMRT | 06:08:35.28 20:39:06.70 | 9 |  14NH3(1,1) | 68 | J214940.65 | TMRT | 21:49:40.65 55:24:51.80 | 9 |  14NH3(1,1) | 62
|  |  |  |  14NH3(2,2) | 65 |  |  |  |  |  14NH3(2,2) | 91
|  |  |  |  14NH3(3,3) | 69 |  |  |  |  |  14NH3(3,3) | 77
|  |  |  |  15NH3(1,1) | 74 |  |  |  |  |  15NH3(1,1) | 96
|  |  |  |  15NH3(2,2) | 75 |  |  |  |  |  15NH3(2,2) | 94
|  |  |  |  15NH3(3,3) | 72 |  |  |  |  |  15NH3(3,3) | 95
G19.88 | TMRT | 18:29:14.25 -11:50:25.55 | 9 |  14NH3(1,1) | 273 | J222122.50 | TMRT | 22:21:22.50 63:51:13.00 | 9 |  14NH3(1,1) | 98
|  |  |  |  14NH3(2,2) | 292 |  |  |  |  |  14NH3(2,2) | 88
|  |  |  |  14NH3(3,3) | 4090 |  |  |  |  |  14NH3(3,3) | 84
|  |  |  |  15NH3(1,1) | 405 |  |  |  |  |  15NH3(1,1) | 71
|  |  |  |  15NH3(2,2) | 356 |  |  |  |  |  15NH3(2,2) | 68
|  |  |  |  15NH3(3,3) | 264 |  |  |  |  |  15NH3(3,3) | 85
G19.92 | TMRT | 18:28:19.17 -11:40:36.38 | 9 |  14NH3(1,1) | 70 | J222851.41 | TMRT | 22:28:51.41 64:13:41.31 | 9 |  14NH3(1,1) | 68
|  |  |  |  14NH3(2,2) | 74 |  |  |  |  |  14NH3(2,2) | 63
|  |  |  |  14NH3(3,3) | 71 |  |  |  |  |  14NH3(3,3) | 98
|  |  |  |  15NH3(1,1) | 70 |  |  |  |  |  15NH3(1,1) | 79
|  |  |  |  15NH3(2,2) | 77 |  |  |  |  |  15NH3(2,2) | 78
|  |  |  |  15NH3(3,3) | 60 |  |  |  |  |  15NH3(3,3) | 95