solar dynamics of e.g. photons, electrons, or other SM particles.

##  7 Bounds and impact on structure formation

Sufficiently strong DM self-interactions alter the early-universe evolution of
DM perturbations, and are constrained by measurements e.g. of the matter power
spectrum, as discussed in [70, 71, 72]. Specifically, from Eq. (12), the
Fourier transform of the DM overdensity field
$\delta\equiv(\rho-\overline{\rho})/\overline{\rho}$ follows [126, 127, 128,
129, 130, 131]

$\ddot{\delta}_{\mathbf{k}}+2H\dot{\delta}_{\mathbf{k}}-4\pi
G\,\overline{\rho}\left[1\pm\left(\frac{k}{k_{\lambda}}\right)^{2}-\left(\frac{k}{k_{J}}\right)^{4}\right]\delta_{\mathbf{k}}=0\,,$
(62)

where $\overline{\rho}\propto a^{-3}$ is the average DM energy density, $k$ is
the comoving momentum, $H\equiv\dot{a}/a$ is the Hubble parameter and $a$ is
the scale factor. The second and third terms in the square bracket of Eq. (62)
modify the standard Meszaros equation and are a consequence of the self-
interactions and the small ULDM mass, and can be interpreted as arising from
an effective sound speed (in the second term, $\pm$ refers to
attractive/repulsive self-interactions). The two comoving momentum scales
$k_{\lambda}(a)$ and $k_{J}(a)$ in Eq. (62) depend on time and are associated
to the self-interactions and to the so-called ‘quantum pressure’ of the ULDM
field. They are defined by

$k_{\lambda}\equiv a\left(\frac{16\pi
Gm^{4}}{|\lambda|}\right)^{\frac{1}{2}}\simeq 2.7{\,{\rm
Mpc}}^{-1}\left[\frac{a}{a_{\rm eq}}\right]\left[\frac{m}{10^{-14}{\,{\rm
eV}}}\right]\left[\frac{f_{a}}{10^{7}{\,{\rm GeV}}}\right]\,,$ (63)

and

$k_{J}\equiv a\left(16\pi
G\overline{\rho}m^{2}\right)^{\frac{1}{4}}\simeq\,10{\,{\rm
Mpc}}^{-1}\left[\frac{a}{a_{\rm
eq}}\right]^{\frac{1}{4}}\left[\frac{m}{10^{-22}{\,{\rm
eV}}}\right]^{\frac{1}{2}}\,,$ (64)

where $a_{\rm eq}$ is the scale factor at the time of matter-radiation
equality. Note that Eq. (62) holds only when the mode $\delta_{\mathbf{k}}$ is
subhorizon, i.e. $k/a>H$.

The self-interaction and quantum pressure terms become dominant and change the
standard growth of the primordial perturbations (which is
$\delta_{\mathbf{k}}\propto 1+(3/2)a/a_{\rm eq}$) for $k$ larger than the
critical momenta $k_{\lambda}$ and $k_{J}$ respectively, while smaller
momentum scales evolve as for standard cold DM. In particular, attractive
self-interactions lead to an exponential increase of $\delta_{\mathbf{k}}$
(until $\delta_{\mathbf{k}}\simeq 1$, at which point non-linear collapse into
bound objects is expected to occur), while repulsive self-interactions and
quantum pressure prevent their growth by making $\delta_{\mathbf{k}}$
oscillate.

As described in Refs. [70, 71], a rough conservative bound is obtained by
requiring that both of these effects are negligible at the largest comoving
momentum at which the matter power spectrum has been measured precisely (from
CMB, large scale structure or Lyman-$\alpha$ data), i.e. approximately
$k\simeq 1{\,{\rm Mpc}}^{-1}$ (see e.g. [132, 133, 134]). Given that
$k_{\lambda}\propto a$ and $k_{J}\propto a^{1/4}$, see Eqs. (63-64), the
corresponding bound on $m$ and $\lambda$ is the strongest at the earliest
times at which Eq. (62) is valid (i.e. the field behaves as DM), which we
assume to be at matter-radiation equality. The blue region of Figure 14 shows
the constraint $k_{\lambda}^{\rm eq}\equiv k_{\lambda}|_{a=a_{\rm eq}}\gtrsim
1{\,{\rm Mpc}}^{-1}$ for ULDM mass and self-coupling in the range of interest
for the formation of the solar halo, while the equivalent constraint from
$k_{J}$ (not shown) only bounds much smaller $m$ (around $10^{-22}{\,{\rm
eV}}$, see Eq. (64)). The self-interactions also affect larger momenta, though
these are poorly constrained by direct observations; should such measurements
improve, the corresponding lower bound on $f_{a}$ in Figure 14 would be
stronger and scale proportionally to $k$. Note that ref. [72] carried out a
much more detailed analysis for repulsive self-interactions, including the
precise evolution of perturbations via a Boltzmann code and comparison with
CMB measurements and large scale structure data. The resulting bound matches
(within $10\%$) with that obtained in our rough estimate.

Figure 14:  Blue region: Approximate constraint on the size of the self-
coupling $\lambda$ (parameterized for convenience by $f_{a}$) from the matter
power spectrum at largest measured momenta, $k\simeq 1{\,{\rm Mpc}}^{-1}$.
This bound applies both for $\lambda>0$ and $\lambda<0$, and considers the
effect of the self-interactions only after matter-radiation equality. The
dotted line is an indication of an upper bound that could be set on $\lambda$
(for attractive self-interactions) from the matter power spectrum due to the
exponential enhancement of the perturbations during radiation domination, see
Eq. (66). Black lines show the values of $m$ and $\lambda$ for which the
typical formation time $|\Gamma|^{-1}$ of the atom in the Solar System equals
$5{\,{\rm Gyr}}$ or $150\,{\rm Myr}$. Over the gray region, the self-coupling
is too weak for the formation to take place within $5{\,{\rm Gyr}}$. We show
in the axis on the right the corresponding value of $\lambda$ parameterized by
$\lambda_{m}\equiv|\lambda|(10^{-14}{\,{\rm eV}}/m)^{2}$. Although the
misalignment mechanism for a conventional axion potential underproduces the DM
abundance for the plotted values of $m$ and $f_{a}$, other scalar potentials
or different production mechanisms could provide the correct abundance.

The self-interactions could in principle modify the evolution of the
perturbations also during radiation domination (this is the case for the so-
called large-misalignment mechanism [131]), potentially strengthening the
bound on $\lambda$. It is known that the gravitational potential induced by
the perturbations in the radiation bath during radiation domination introduces
a source term into Eq. (62), which provides a (small) non-zero value for the
time-derivative of $\delta_{\mathbf{k}}$ once the modes enter the horizon, so
that $\delta_{\mathbf{k}}$ increases logarithmically from the time when the
mode enters the horizon until matter-radiation equality. In particular, the
mode at $k=1{\,{\rm Mpc}}^{-1}$ enters the horizon at redshift $z_{i}\simeq
5\cdot 10^{5}$, corresponding to the temperature $T_{i}\simeq 100{\,{\rm
eV}}$. In the following we briefly discuss how the self-interactions could
affect this behavior, and give a crude estimate of the corresponding bound
from the matter power spectrum.

For attractive self-interactions, Eq. (62) can be rewritten as

$\delta_{\mathbf{k}}^{\prime\prime}+\frac{3y+2}{2y(y+1)}\delta_{\mathbf{k}}^{\prime}-\frac{3}{2y(y+1)}\left[1+\left(\frac{k}{k_{\lambda}^{\rm
eq}}\right)^{2}\frac{1}{y^{2}}\right]\delta_{\mathbf{k}}=0\,,$ (65)

where $y\equiv a/a_{\rm eq}$, $\delta^{\prime}\equiv d\delta/dy$ and we
neglected the quantum pressure term, which is an excellent approximation for
$m\gtrsim 10^{-21}{\,{\rm eV}}$. The solution of Eq. (65) deep into radiation
domination starting from $a=a_{i}$ (when the perturbation enters the horizon)
with initial conditions $\delta_{\mathbf{k}}(a_{i})=1$ and
$\dot{\delta}_{\mathbf{k}}(a_{i})=0$ is

$\delta_{\mathbf{k}}(a)=\left[\frac{a}{a_{i}}\right]^{\frac{1}{4}}\cosh\left[\sqrt{6\frac{a_{\rm
eq}}{a_{i}}}\frac{k}{k_{\lambda}^{\rm
eq}}\left(1-\sqrt{\frac{a_{i}}{a}}\right)\right]\,.$ (66)

Eq. (66) shows that, soon after the time when $a=a_{i}$, perturbations with
momentum $k\gtrsim k_{\lambda}^{\rm eq}\sqrt{a_{i}/a_{\rm eq}}$ increase
substantially during radiation domination, receiving an enhancement factor
much larger than one, in particular of order $(a_{\rm
eq}/a_{i})^{1/4}\exp{[\sqrt{6a_{\rm eq}/a_{i}}k/k_{\lambda}^{\rm eq}]}$ by the
time of matter-radiation equality.292929The expression in Eq. (66) is a
solution only in the limit $k\gtrsim k_{\lambda}^{\rm eq}\sqrt{a/a_{\rm eq}}$,
and overestimates the true solution otherwise. Eq. (66) has been obtained by
considering the leading $1/y$ terms in Eq. (65) and, in particular, only the
second term in the square bracket, which is a good approximation for $k\gtrsim
k_{\lambda}^{\rm eq}\sqrt{a_{i}/a_{\rm eq}}$.  Measurements of the matter
power spectrum can rule out such an increase and provide a (stronger) bound on
the attractive self-coupling. This can be estimated by requiring
$k_{\lambda}^{\rm eq}\sqrt{a_{i}/a_{\rm eq}}\gtrsim 1{\,{\rm Mpc}}^{-1}$, with
$a_{i}=1/(1+z_{i})$. The corresponding upper limit on $f_{a}$, shown with
dotted lines in Figure 14, is stronger by a factor $1/\sqrt{a_{i}/a_{\rm
eq}}\simeq 10$ than the one discussed previously. Given the roughness of this
analysis, we only take this bound as indicative, and leave a full analysis
(which could differ from this estimate) for future work.

On the other hand, if the self-interactions are repulsive, the solution of Eq.
(62) does not increase, but oscillates at momenta $k\gtrsim k_{\lambda}^{\rm
eq}\sqrt{a_{i}/a_{\rm eq}}$. Therefore, in this case, the effect of the self-
interactions before matter-radiation equality is much less dramatic. Indeed
the analysis of ref. [72], carried out from very high redshift, finds a bound
comparable to that at matter-radiation equality, as discussed above.

Finally, black and red lines in Figure 14 represent the values of $m$ and
$f_{a}$ for which the formation time scale $|\Gamma|^{-1}$ is $5{\,{\rm Gyr}}$
and $150\,{\rm Myr}$ respectively: For $m\lesssim 10^{-14}{\,{\rm eV}}$,
$\Gamma$ is given by Eq. (8), otherwise by Eq. (7); the transition region is
uncertain. In particular, in the gray region the self-coupling is too small
for the formation to occur within $5{\,{\rm Gyr}}$. We conclude that
gravitational atoms can form within the age of the Solar System consistently
with direct bounds from structure formation, for about one or two orders of
magnitude in $f_{a}$, in the range $10^{6}\div 10^{8}{\,{\rm GeV}}$.

#### Other limits on ULDM self-interactions

Reference [135] obtained a bound on the self-interactions that is stronger
than that of Figure 14 by requiring that the ULDM equation of state is
$w=p/\rho<10^{-3}$ at matter-radiation equality, where $p$ is the pressure,
though the choice of this threshold is artificial. Additionally, any
effectively massless scalar field, with mass $m<H$, gets isocurvature
perturbations during inflation, which are strongly constrained by cosmic
microwave background (CMB) observations. However, these depend on the initial
field value during inflation, and are evaded e.g. if the initial field value
is large (e.g. of order $M_{\rm Pl}$), if the theory is a different phase or
the scalar field gets a large effective mass via couplings to the inflation
field; see also [70].

Finally, the existence of ultralight particles, regardless of them being DM,
induces the spinning down of rotating black holes – a phenomenon known as
superradiance – and is constrained by the observation of such spinning objects
[136, 137]. However, the self-interactions suppress the superradiance process
for mass and couplings in Figure 14 which are relevant for the formation of
the gravitational atom [138, 139].

##  8 Conclusions and outlook

We have shown that a halo bound to massive astrophysical bodies, resembling a
‘gravitational atom’, forms from the capture of ultralight (pseudo-)scalar
dark matter (DM) particles, as a consequence of their (weak) quartic self-
interactions. This formation takes place also around our Sun, leading to a
_solar halo_ whose density can be orders of magnitude larger than the
background DM density at the position of the Earth. More generally, the
formation of gravitational atoms around astrophysical objects in our Galaxy
could provide novel ULDM signatures.

The efficiency of the capture process is controlled by the amount of
gravitational focusing of the galactic DM waves around the body, which is
measured by the ratio $\xi_{\rm foc}\equiv\lambda_{\rm
dB}/R_{\star}=2\pi\alpha/v_{\rm dm}$ between their de Broglie wavelength,
$\lambda_{\rm dB}$, and the gravitational Bohr radius, $R_{\star}$ in Eq. (2).
In particular, the bound state either saturates to an (underdense) steady-
state configuration (when $\xi_{\rm foc}\lesssim 1$, i.e. $m\lesssim
10^{-14}{\,{\rm eV}}$ for the Sun) or it exponentially grows in density (when
$\xi_{\rm foc}\gtrsim 1$). We have confirmed this behavior through numerical
simulations. On the other hand, the size of the self-interaction coupling
$\lambda$ only dictates the formation time scale, which in the regime
$\xi_{\rm foc}\gtrsim 1$ is of the order of the relaxation time via self-
interactions, defined in Eq. (47). The macroscopic properties of the resulting
bound state, e.g. its radius and maximum density after the exponential
increase (in Eq. (9)), depend on the particle mass $m$, the central (for
example, solar) mass $M$, the properties of the galactic DM near the system
(through $\rho_{\rm dm}$ and $v_{\rm dm}$), and the magnitude of $\lambda$.
Note that the basic mechanism depends neither on the sign of $\lambda$ nor the
parity of the field.

The most striking implications occur in the exponential-growth regime, when
$\xi_{\rm foc}\gtrsim 1$. For instance, for axion DM with mass $m\simeq{\rm
few}\cdot 10^{-14}$ eV and decay constants $f_{a}\simeq 10^{7}\div 10^{8}$
GeV, on Gyr timescales the solar halo density grows to be larger than the
background DM density by a factor as large as $10^{4}$ at a distance
$r=1\,{\rm AU}$ from the Sun; see Figures 4 and 6. The density at $r\ll$ AU
can be even larger for smaller masses, $m\simeq 10^{-13}$ eV, (see Figure 13)
providing a clear target for proposed searches in space [117]. In the presence
of low-velocity and low-dispersion DM substructure, like a dark disk, the
density could be further enhanced.

The exponential growth stops once the density of the bound states approaches
the critical density in Eq. (56). For repulsive self-interactions
($\lambda>0$), at this point the gravitational atom becomes a stable bound
object, whose density can be large compared to that of the background DM.
Instead, for attractive self-interactions ($\lambda<0$), the atom is unstable
and collapses, possibly triggering a Bosenova emission of relativistic
particles. Both the large density enhancements and a possible nearby Bosenova
can provide novel signals in ULDM detection experiments.

In this work we assumed that the gravitational potential is proportional to
$1/r$, which implies that $R_{\star}$ is larger than the radius $R$ of the
astrophysical body that captures the ULDM. For large enough axion masses
($m\gtrsim 2\cdot 10^{-13}{\,{\rm eV}}$ for the solar halo, see Figure 3),
this assumption is no longer satisfied. In this case, assuming a uniform
density inside the body, the potential is proportional to $r^{2}$ rather than
$1/r$, and, if captured, the particles occupy bound states _inside_ the bulk
of the astrophysical body. A full analysis of the formation of bound states in
this case is worthy of further study.

We emphasize that the formation mechanism introduced in this paper is generic,
requiring only a light DM (pseudo-)scalar field, with sufficiently strong
quartic self-interactions for the formation to occur within Gyr time scales.
In particular, this leads to conclude that gravitational atoms form around any
sufficiently compact object (i.e. with radius smaller than $R_{\star}$). Here
we merely note that when $\xi_{\rm foc}\gtrsim 1$, the Bohr radius $R_{\star}$
of the atom bound to a few different astrophysical systems is given by:

$R_{\star}\lesssim\begin{cases}\displaystyle\,\,\,1.3\,R_{J}\left[\frac{1.7\cdot
10^{-11}{\,{\rm
eV}}}{m}\right]^{2}\left[\frac{10^{-3}M_{\odot}}{M}\right]&\text{(Jupiter-
like\,\,planet)}\\\ \displaystyle\,\,\,\,1\,{\,{\rm AU}}\,\left[\frac{1.7\cdot
10^{-14}{\,{\rm
eV}}}{m}\right]^{2}\left[\frac{M_{\odot}}{M}\right]&{\rm(Star\,\,or\,\,neutron\,\,star)}\\\
\displaystyle\,\,\,30\,\,{\rm pc}\,\left[\frac{1.7\cdot 10^{-21}{\,{\rm
eV}}}{m}\right]^{2}\left[\frac{10^{7}M_{\odot}}{M}\right]&\text{(Supermassive
black hole)}\end{cases}\,$ (67)

where for reference $R_{J}=7\cdot 10^{4}{\,{\rm km}}$ is the radius of
Jupiter. Therefore, neutron stars, black holes (both solar-mass and
supermassive), and planets are among promising targets for a dedicated study.
Note that the gravitational atom bound to the Earth would have a radius larger
than that of the Earth only for $m\lesssim 10^{-9}$ eV [44]. However, for
these masses and $M=M_{\oplus}$ (the mass of the Earth), $\xi_{\rm
foc}\lesssim 0.1$ and therefore the exponential growth does not occur. Using
Eq. (9) one finds that $\rho_{\rm eq}\ll\rho_{\rm dm}$ for all relevant masses
$m$.

It is worthwhile to compare our results to other known configurations of
(light) particles bound to astrophysical objects. The superradiant instability
of ultralight fields around rapidly rotating black holes can populate $l>0$
modes by extracting the angular momentum of the host (see, for example, [136,
137, 138, 139]). Contrary to the mechanism discussed in this paper,
superradiance does not require the field to be the DM. Additionally, the
process is effective when $m$ is of the order of the black hole radius
$1/(2GM)$, and therefore, for fixed $M$, is active for larger boson masses (by
a factor $2\pi/v_{\rm dm}$) than those for which DM capture starts to be
effective, $\xi_{\rm foc}\simeq 1$.303030The analysis of solar-mass and
supermassive spinning black holes led to constraints in the mass range
$m\simeq 10^{-13}\div 10^{-11}$ eV [140, 138] and $m\simeq 10^{-21}\div
10^{-17}{\,{\rm eV}}$ [141], respectively. Similar processes can occur in
dense stellar objects like neutron stars, but require additional SM couplings
beyond gravity and are not relevant for less dense objects like our Sun [142].
Finally, the particles produced by the Sun that lie in the low-velocity tail
of the distribution can get captured into bound orbits around the Sun, giving
rise to yet-another bound configuration, a _solar basin_ [124, 125, 143]. This
is however most relevant for particle masses $m\sim\mathcal{O}({\,{\rm
keV}})$, much larger than those considered here. Note that as with
superradiance, the solar basin will form even if the field constitutes a
negligible fraction of the DM.

Self-gravitating bound configurations of ultralight bosons, known as boson
stars [27, 28, 35], can form on astrophysical timescales (see e.g. [30, 31]).
After formation, their mass $M_{\rm bs}$ changes in time through processes
mediated by gravitational interactions (see e.g. [34, 144]), similar in form
to those that occur in the formation of the gravitational atom discussed in
this paper. The authors of ref. [69, 145] analyzed the mass growth of boson
stars defining the characteristic parameter $\nu\equiv v_{\rm dm}/v_{b}$
($\nu$ is analogous to $2\pi/\xi_{\rm foc}$): Boson stars with $\nu\gg 1$
decrease in mass as a result of processes analogous to _stripping_ , see
Figure 8(c). In the other limit, $\nu\ll 1$, their mass increases, but at a
rate that decreases in time, as $M_{\rm bs}$ becomes larger. While there is a
close connection between the conditions for the growth of the gravitational
atom and the boson stars, the dynamical evolution of the system is different
in the two cases. Vector boson stars, similar configurations for spin-$1$ ULDM
particles, can also form during the evolution of the Universe [39, 40, 41, 42,
43]. Finally, we observe that the gravitational atoms could form also from the
capture of such spin-1 DM, although the capture process could be different
given the derivative form of the self-interactions.

## Acknowledgements

The authors are grateful to Hyungjin Kim and Christopher McKee for crucial
contributions in the early stages of this work. The authors thank Kfir Blum,
Elisa Ferreira, Ricardo Ferreira, Ed Hardy, Oleksii Matsedonskyi, Mehrdad
Mirbabayi, Wolfram Ratzinger, Giovanni Villadoro, and Wei Xue for valuable
discussions. The authors also thank Kfir Blum and Ed Hardy for insightful
feedbacks on a draft. MG is grateful to Ed Hardy for substantial contributions
to the Schrödinger–Poisson code used for the numerical results of this paper.
The work of DB was supported in part by the Cluster of Excellence “Precision
Physics, Fundamental Interactions, and Structure of Matter” (PRISMA+ EXC
2118/1) and the German-Israeli Foundation (GIF). The work of JE was supported
by the World Premier International Research Center Initiative (WPI), MEXT,
Japan and by the JSPS KAKENHI Grant Numbers 21H05451 and 21K20366. The work of
MJ is supported by a research grant from the Musk Foundation. The work of GP
is supported by grants from BSF-NSF, Friedrich Wilhelm Bessel research award
of the Alexander von Humboldt Foundation, GIF, ISF, Minerva, SABRA - Yeda-Sela
- WRC Program, the Estate of Emile Mimran, and the Maurice and Vivienne Wohl
Endowment.

##  Appendix A Energy density and bound states

In this Appendix we derive the energy density $\epsilon$ in Eq. (16) and
provide the expression of excited states $\psi_{nlm}$.

##### Energy density.

As mentioned in Section 4.2, in the nonrelativistic limit the action of the
scalar field reduces to $S\equiv\int dtd^{3}x\mathcal{L}$ in Eq. (42). If the
self-gravitational potential $\Phi_{\rm se}$ (neglected in the main text) is
relevant, $\mathcal{L}$ contains the additional terms $-(8\pi
G)^{-1}|\nabla\Phi_{\rm se}|^{2}-m\Phi_{\rm se}|\psi|^{2}$. The invariance of
$S$ under time translations implies that the time-like component ‘00’ of the
(nonrelativistic) energy momentum tensor, which is given by

$\frac{\partial\mathcal{L}}{\partial\dot{\psi}}\dot{\psi}+\frac{\partial\mathcal{L}}{\partial\dot{\psi}^{*}}\dot{\psi}^{*}+\frac{\partial\mathcal{L}}{\partial\dot{\Phi}}\dot{\Phi}-\mathcal{L}=\frac{1}{2m}|\nabla\psi|^{2}+\frac{1}{8\pi
G}|\nabla\Phi_{\rm se}|^{2}\,+m(\Phi_{\rm ex}+\Phi_{\rm
se})|\psi|^{2}+\frac{1}{2}g|\psi|^{4}\,,$ (68)

(with $\Phi=\Phi_{\rm se}+\Phi_{\rm ex}$) is conserved on the EoM, when
integrated over the whole space. Eq. (68) can be simplified using
$|\nabla\Phi_{\rm se}|^{2}=\nabla\cdot(\Phi_{\rm se}\nabla\Phi_{\rm
se})-\Phi_{\rm se}\nabla^{2}\Phi_{\rm se}$: The first term vanishes upon
spatial integration, while the second can be rewritten using the Poisson
equation (see footnote 6). This leads to the conclusion that the energy
density

$\epsilon=\frac{|\nabla\psi|^{2}}{2m}+m\left(\Phi_{\rm ex}+\frac{\Phi_{\rm
se}}{2}\right)|\psi|^{2}+\frac{1}{2}g|\psi|^{4}\,\,,$ (69)

is globally conserved.

##### Bound states.

The quasi-stationary solutions $\psi=e^{-i\omega t}\Psi(\mathbf{x})$ of the
EoM in Eq. (15) satisfy the eigenvalue equation
$(-\nabla^{2}/2m-\alpha/r)\Psi=\omega\Psi$, whose most general solutions with
negative energy ($\omega<0$) are the well-known hydrogen-like bound states

$\Psi=\psi_{nlm}=R_{nl}(r)Y_{l}^{m}(\theta,\varphi)\,,\quad
R_{nl}(r)\equiv\left[\frac{2}{nR_{\star}}\right]^{\frac{3}{2}}\left[\frac{2r}{nR_{\star}}\right]^{l}\sqrt{\frac{(n-l-1)!}{2n(l+n)!}}L_{n-l-1}^{2l+1}\left[\frac{2r}{nR_{\star}}\right]e^{-\frac{r}{nR_{\star}}}\,,$
(70)

with eigenvalues $\omega=\omega_{n}=-m\alpha^{2}/(2n^{2})$. In Eq. (70),
$Y_{l}^{m}(\theta,\varphi)$ are the spherical harmonics and $L_{a}^{b}[x]$ the
generalized Laguerre polynomials of degree $a$ and order $b$. The bound states
$\psi_{nlm}$ and the scattering states $\psi_{\mathbf{k}}$ in Eq. (18) are
orthonormalized, and are all orthogonal to one another, as

$\displaystyle\int
d^{3}x\,\psi^{*}_{nlm}\psi_{n^{\prime}l^{\prime}m^{\prime}}$
$\displaystyle=\delta_{nn^{\prime}}\delta_{ll^{\prime}}\delta_{mm^{\prime}}\,,\quad\int
d^{3}x\psi_{{\bf k}^{\prime}}^{*}\psi_{\bf k}=(2\pi)^{3}\delta^{3}({\bf
k}-{\bf k}^{\prime})\,,\quad\int d^{3}x\psi^{*}_{nlm}\psi_{\mathbf{k}}=0\,.$
(71)

##  Appendix B Perturbative solution of the Gross–Pitaevskii equation

In this Appendix we provide a detailed derivation of the change in the number
of bound particles $\langle\dot{N}_{nlm}\rangle$ in Eq. (41). This is obtained
by solving the equation of motion

$\left(i\partial_{t}+\frac{\nabla^{2}}{2m}+\frac{\alpha}{r}\right)\psi=g|\psi|^{2}\psi\,,$
(72)

with initial condition $\psi(t=0)=(\psi_{w}+\psi_{b})_{t=0}$ in Eqs. (23) and
(35), perturbatively in the (local) effective self-interaction parameter
$g_{\rm eff}=g\rho/(m\omega)$, where $\rho$ and $\omega$ are the typical
values of the local density and frequency of the field. As mentioned in
Section 4.1, in the regime $\rho<\rho_{\rm crit}$, the self-interactions are
small enough that it is possible to approximate the generic solution of Eq.
(72) as a superposition of solutions of the EoM for $\lambda=g=0$, i.e.

$\psi=\sum_{nlm}c_{nlm}(t)e^{-i\tilde{\omega}_{nlm}t}\psi_{nlm}+\int[dk]c_{\mathbf{k}}(t)e^{-i\tilde{\omega}_{\mathbf{k}}t}\psi_{\mathbf{k}}\,,$
(73)

where $\psi_{nlm}$ is given in Eq. (70) and the coefficients $c_{nlm}$ and
$c_{\mathbf{k}}$ depend on $\lambda$ (or $g$) and vary in time more slowly
than the typical field oscillation periods, given by $1/\tilde{\omega}_{nlm}$
and $1/\tilde{\omega}_{\mathbf{k}}$ for bound and scattering states
respectively. We can expand $c_{nlm}$ and $c_{\mathbf{k}}$ in a series in
$g_{\rm eff}$ as

$c_{nlm}(t)=c_{nlm}^{(0)}(t)+c_{nlm}^{(1)}(t)+c_{nlm}^{(2)}(t)+\ldots,\qquad
c_{\mathbf{k}}(t)=c_{\mathbf{k}}^{(0)}(t)+c_{\mathbf{k}}^{(1)}(t)+c_{\mathbf{k}}^{(2)}(t)+\ldots\,,$
(74)

where the superscripts indicate the order in $g_{\rm eff}$ of each term. In
general $\tilde{\omega}_{nlm}$ and $\tilde{\omega}_{\mathbf{k}}$ in Eq. (73)
are different from the original eigenfrequencies $\omega_{n}$ and
$\omega_{k}$, which indeed receive time-independent corrections (or ‘shifts’)
when $\lambda\neq 0$; see for instance [138]. These corrections may be also
expanded as

$\tilde{\omega}_{nlm}=\omega_{n}+\omega_{nlm}^{(1)}+\omega_{nlm}^{(2)}+\ldots,\qquad\tilde{\omega}_{\mathbf{k}}=\omega_{k}+\omega_{\mathbf{k}}^{(1)}+\omega_{\mathbf{k}}^{(2)}+\ldots,$
(75)

where $\omega_{n}=-m\alpha^{2}/2n^{2}$ and $\omega_{k}=k^{2}/2m$ are the
original $\lambda=0$ frequencies, and $|\omega_{nlm}^{(1)}|\ll|\omega_{n}|$,
etc.. As a result, we can organize the full solution $\psi$ in Eq. (73) in the
perturbative series

$\psi=\psi^{(0)}+\psi^{(1)}+\psi^{(2)}+\ldots\,,$

where, for instance,

$\displaystyle\psi^{(0)}$
$\displaystyle=\sum_{nlm}c_{nlm}^{(0)}(t)e^{-i\omega_{n}t}\psi_{nlm}+\int[dk]c_{\mathbf{k}}^{(0)}(t)e^{-i\omega_{k}t}\psi_{\mathbf{k}}\
,$ (76) $\displaystyle\psi^{(1)}$
$\displaystyle=\sum_{nlm}\left[c_{nlm}^{(1)}(t)-i\omega_{nlm}^{(1)}tc_{nlm}^{(0)}(t)\right]e^{-i\omega_{n}t}\psi_{nlm}+\int[dk]\left[c_{\mathbf{k}}^{(1)}(t)-i\omega_{\mathbf{k}}^{(1)}tc_{\mathbf{k}}^{(0)}(t)\right]e^{-i\omega_{k}t}\psi_{\mathbf{k}}\
.$ (77)

The expansions above allow us to easily solve Eq. (72) at every order in the
perturbative expansion. First, at the zeroth order Eq. (72) simply reads

$\left(i\partial_{t}+\frac{\nabla^{2}}{2m}+\frac{\alpha}{r}\right)\psi^{(0)}=0\,,$
(78)

and implies that $c_{nlm}^{(0)}(t)$ and $c_{\mathbf{k}}^{(0)}(t)$ are in fact
time-independent. Imposing that $\psi$ equals the initial condition
$\psi_{w}+\psi_{b}$ at $t=0$, we find $\psi^{(0)}=\psi_{w}+\psi_{b}$, and
$c_{nlm}^{(0)}=\sqrt{N_{nlm}^{(0)}}$ and $c_{\mathbf{k}}^{(0)}=a(\mathbf{k})$.
The zeroth-order solution in fact coincides with the initial conditions, and
the full nontrivial time evolution is determined by the time-dependent
coefficients $c_{nlm}^{(i)}$ and $c_{\mathbf{k}}^{(i)}$.

As discussed in Section 4.1, the number of particles bound to the $nlm$ level
follows Eq. (36), reported here for simplicity:

$N_{nlm}(t)=|c_{nlm}(t)|^{2}=N_{nlm}^{(0)}+2\sqrt{N_{nlm}^{(0)}}{\rm
Re}[c_{nlm}^{(1)}(t)]+\left\\{|c_{nlm}^{(1)}(t)|^{2}+2\sqrt{N_{nlm}^{(0)}}{\rm
Re}[c_{nlm}^{(2)}(t)]\right\\}+\dots\,.$ (79)

In the following we will describe the solution of the EoM (72) up to the
second order, which will be sufficient to calculate the first non-trivial part
of $\langle\dot{N}_{nlm}\rangle$.

#### First-order solution

The first order EoM is

$\left(i\partial_{t}+\frac{\nabla^{2}}{2m}+\frac{\alpha}{r}\right)\psi^{(1)}=g|\psi^{(0)}|^{2}\psi^{(0)}\,,$
(80)

which is a source equation for $\psi^{(1)}$. We solve Eq. (80) with initial
conditions containing a generic configuration of scattering states
$a(\mathbf{k})$, but only particles in the ground state, i.e.
$N_{nlm}^{(0)}=0$ unless $nlm=100$, or equivalently

$\psi^{(0)}=\sqrt{N_{100}^{(0)}}e^{-i\omega_{1}t}\psi_{100}+\int[dk]a({\mathbf{k}})e^{-i\omega_{k}t}\psi_{\mathbf{k}}\,.$
(81)

With these initial conditions, $\langle\dot{N}_{100}\rangle$ will describe the
evolution of the bound-state occupation number both before the ground state is
populated (upon taking $N_{100}^{(0)}=0$), and after a ground-state population
has been captured. Plugging the expression of $\psi^{(1)}$ of Eq. (77) into
Eq. (80), we obtain the time evolution of the first-order coefficients:

$\sum_{nlm}\left[i\dot{c}_{nlm}^{(1)}(t)+\sqrt{N_{nlm}^{(0)}}\omega_{nlm}^{(1)}\right]e^{-i\omega_{n}t}\psi_{nlm}+\int[dk]\left[c_{\mathbf{k}}^{(1)}(t)+a(\mathbf{k})\omega_{\mathbf{k}}^{(1)}\right]e^{-i\omega_{\mathbf{k}}t}\psi_{\mathbf{k}}=g|\psi^{(0)}|^{2}\psi^{(0)}\,.$
(82)

Using the orthonormality relations of Eq. (71) we can single out the
contribution of the individual bound-state levels and waves:

$\displaystyle i\dot{c}_{nlm}^{(1)}(t)+\sqrt{N_{nlm}^{(0)}}\omega_{nlm}^{(1)}\
$ $\displaystyle=ge^{i\omega_{n}t}\int
d^{3}x|\psi^{(0)}|^{2}\psi^{(0)}\psi_{nlm}^{*}\,,$ (83) $\displaystyle
i\dot{c}_{\mathbf{k}}^{(1)}+a(\mathbf{k})\omega_{\mathbf{k}}^{(1)}$
$\displaystyle=ge^{i\omega_{k}t}\int
d^{3}x|\psi^{(0)}|^{2}\psi^{(0)}\psi_{\mathbf{k}}^{*}\,.$ (84)

Substituting the expression of $\psi^{(0)}$ of Eq. (81) into Eqs. (83) and
(84) and expanding the products, we obtain multiple terms on the right hand
side (RHS) of these equations with different physical meaning. Crucially, a
generic such ‘source’ term $S(\mathbf{x})e^{-i\omega_{s}t}\subset
g|\psi^{(0)}|^{2}\psi^{(0)}$ that oscillates with frequency $\omega_{s}$
contributes to a significant change in $c_{nlm}$ or $c_{\mathbf{k}}$ only if
the frequencies in the phases in the RHS of Eqs. (83) and (84) sum up to zero,
in such a way that $c_{nlm}$ or $c_{\mathbf{k}}$ have a chance of growing with
$t$ instead of oscillating. In particular (see also [138]):

  1. 1.

In Eq. (83), terms with $\omega_{s}<0$ and $\omega_{s}=\omega_{n}$ for some
$n$ correspond to a ‘resonant’ increase of the bound modes labeled by $n$
(i.e. production of particles in the state $n$), or, as we will see later, to
a correction of the eigenfrequency $\omega_{nlm}^{(1)}$. If instead
$\omega_{s}\neq\omega_{n}$ for all $n$, the term is a driven off-resonance
oscillator, and represents a (suppressed) off-resonant production of
particles.

  2. 2.

In Eq. (84), terms with $\omega_{s}>0$ correspond to a resonant production of
scattering states or a correction to their frequency
$\omega_{\mathbf{k}}^{(1)}$.

As we will see, only resonant terms will effectively contribute to $c_{nlm}$
and $c_{\mathbf{k}}$. Let us first focus on the ground state, $nlm=100$. There
are five types of terms on the RHS of Eq. (83). As Eq. (34) in the main text,
these can be pictured diagrammatically as in Figure 15, where the states ‘100’
and $\mathbf{k}$ correspond to the factors of $\psi_{100}$ or
$\psi_{\mathbf{k}}$ appearing in the integrand of each term, and are placed on
the right or left depending on whether the factor appears with or without the
complex conjugate. The dashed line represents the $\psi_{100}^{*}$ (or in
general $\psi_{nlm}^{*}$) factor common in all the terms. This representation
makes it clear that the 100 level is populated, together with the other state
that enters with the conjugate. We will now consider each diagram in turn:

$100$$100$$100$$100$$(1)$$\mathbf{k}_{1}$$100$$100$$100$$(2)$$\mathbf{k}_{1}$$100$$\mathbf{k}_{2}$$100$$(3)$$\mathbf{k}_{1}$$\mathbf{k}_{2}$$100$$100$$(4)$$\mathbf{k}_{1}$$\mathbf{k}_{2}$$\mathbf{k}_{3}$$100$$(5)$
Figure 15: Diagrammatic representation of the terms in the right hand side of
Eq. (83) that could contribute to the time-evolution of the coefficient
$c_{100}^{(1)}$. Eventually only diagram (5) leads to a significant change in
$c_{100}^{(1)}$, as all other diagrams do not conserve energy or do not change
the particle number in the 100 state.

  * •

Diagram (4) in Figure 15 corresponds to

$g\sqrt{N_{100}^{(0)}}\int[dk_{1}][dk_{2}]e^{i(2\omega_{1}-\omega_{k_{1}}-\omega_{k_{2}})t}a(\mathbf{k}_{1})a(\mathbf{k}_{2})\int
d^{3}x\psi_{\mathbf{k}_{1}}\psi_{\mathbf{k}_{2}}{\psi_{100}^{*2}}\,$

and leads to an off-resonance oscillation since the argument in the
exponential is always negative. This also happens for diagram (2), whose
integrand is proportional to $e^{i(\omega_{1}-\omega_{k_{1}})t}$. Intuitively,
for these terms the energy in the corresponding processes is not conserved for
all values of $k_{1}$ and $k_{2}$, and the production of bound states is not
effective. Instead, all the other diagrams have a chance of obtaining a
vanishing phase.

  * •

Diagram (1) corresponds to $gN_{100}^{(0)3/2}\int d^{3}x|\psi_{100}|^{4}$ and
is time-independent, and would in principle lead to a linear increase of
$c^{(1)}_{100}$. In reality, this term is cancelled by the term
$\sqrt{N_{nlm}^{(0)}}\omega_{nlm}^{(1)}$ on the left hand side of Eq. (83),
because the frequency shift $\omega_{100}^{(1)}$ contains the contribution
$gN_{100}^{(0)}\int d^{3}x|\psi_{100}|^{4}=gN_{100}^{(0)}/8\pi R_{\star}^{3}$.
Consequently, diagram (1) represents the change in energy of the $100$ state
due to self-interactions among the particles present in the ground state, and
does not actually change their occupation number. Note that the condition
$|\omega_{100}^{(1)}|<|\omega_{1}|$ for our perturbative method to work can be
translated to $N_{100}^{(0)}<4\pi/gm^{2}\alpha$. Equivalently, this can be
expressed in terms of $\rho_{\rm crit}$ defined in Eq. (56) as
$\rho(r=0)<2\rho_{\rm crit}$.

  * •

Diagram (3) corresponds to

$2g\sqrt{N_{100}^{(0)}}\int[dk_{1}][dk_{2}]e^{-i(\omega_{k_{1}}-\omega_{k_{2}})t}a^{*}(\mathbf{k}_{2})a(\mathbf{k}_{1})\int
d^{3}x|\psi_{100}|^{2}\psi_{\mathbf{k}_{1}}\psi^{*}_{\mathbf{k}_{2}}\,.$

Although the phase could vanish when $k_{1}=k_{2}$, this term also does not
lead to a change in the number of particles $N_{100}$. To understand its
physical interpretation, it is convenient to write
$a^{*}(\mathbf{k}_{2})a(\mathbf{k}_{1})$ as $\langle
a^{*}(\mathbf{k}_{2})a(\mathbf{k}_{1})\rangle+[a^{*}(\mathbf{k}_{2})a(\mathbf{k}_{1})-\langle
a^{*}(\mathbf{k}_{2})a(\mathbf{k}_{1})\rangle]$. Using Eq. (24), the first
part becomes

$2g\sqrt{N_{100}^{(0)}}\int[dk]f(\mathbf{k})\int
d^{3}x|\psi_{100}|^{2}|\psi_{\mathbf{k}}|^{2},$

and thus represents a shift of the energy of the ground state due to the dark
matter waves, $\omega_{100}^{(1)}\supset 2g\int[dk]f(\mathbf{k})\int
d^{3}x|\psi_{100}|^{2}|\psi_{\mathbf{k}}|^{2}$. (For instance, in the plane
wave limit of $\psi_{\mathbf{k}}$, the frequency correction is
$\omega_{100}^{(1)}\supset 2gn_{\rm dm}$, with $n_{\rm dm}=\rho_{\rm dm}/m$
being the number density of the waves. Thus the perturbative condition
$|\omega_{100}^{(1)}|<|\omega_{1}|$ is equivalent $\rho_{\rm dm}<\rho_{\rm
crit}/8$.) The second part is proportional to
$a^{*}(\mathbf{k}_{2})a(\mathbf{k}_{1})-\langle
a^{*}(\mathbf{k}_{2})a(\mathbf{k}_{1})\rangle$, and describes the elastic
scattering of the waves with the bound state, and vanishes upon taking the
average over times much longer than the coherence time, so does not contribute
to $\langle\dot{N}_{100}\rangle$ at this order. Although it actually affects
$\dot{c}_{100}$, it can be verified that its contribution to $N_{100}(t)$ at
second order in $g_{\rm eff}$ is canceled between
$\langle|c_{100}^{(1)}|^{2}\rangle$ and $2\sqrt{N_{100}^{(0)}}{\rm Re}[\langle
c_{100}^{(2)}\rangle]$ in Eq. (89). We thus ignore this term in the following.

  * •

Finally, diagram (5) corresponds to a resonant production of the $100$ state
via the ‘capture’ process ${k}_{1}+{k}_{2}\rightarrow{k}_{3}+100$, and is the
only diagram in Figure 15 that leads to an effective change in
${c}_{100}^{(1)}$.

As a result the part of the solution of Eq. (83) that grows in time is, for
$nlm=100$,

$\displaystyle c_{100}^{(1)}(t)$
$\displaystyle=-ig\int[dk_{1}][dk_{2}][dk_{3}]a^{*}(\mathbf{k}_{3})a(\mathbf{k}_{1})a(\mathbf{k}_{2})\mathcal{M}_{{k}_{1}+{k}_{2}\rightarrow{k}_{3}+100}\int_{0}^{t}e^{i(\omega_{1}+\omega_{k_{3}}-\omega_{k_{1}}-\omega_{k_{2}})t^{\prime}}dt^{\prime}\,,$
(85)

where

$\mathcal{M}_{a+b\rightarrow
c+d}\equiv\int\psi_{a}\psi_{b}\psi_{c}^{*}\psi_{d}^{*}d^{3}x\,.$

This reproduces Eq. (30) in the main text, and, correspondingly, diagram (5)
reproduces $c_{nlm}^{(1)}$ in Eq. (34).

$\mathbf{k}_{1}$$\mathbf{k}_{2}$$\mathbf{k}_{3}$$nlm$$(1)$$100$$\mathbf{k}_{1}$$\mathbf{k}_{2}$$nlm$$(2)$$\mathbf{k}_{1}$$\mathbf{k}_{2}$$\mathbf{k}_{3}$$\mathbf{k}$$(1)$$100$$\mathbf{k}_{1}$$\mathbf{k}_{2}$$\mathbf{k}$$(2)$$\mathbf{k}_{1}$$\mathbf{k}_{2}$$100$$\mathbf{k}$$(3)$
Figure 16:  Diagrammatic representation of the terms in the right hand side of
Eqs. (83) and (84) that are relevant for the change in $c_{nlm}^{(1)}$ (left)
and $c_{\mathbf{k}}^{(1)}$ (right). Only diagram (2, left) and diagrams (2-3,
right) end up contributing to $\langle\dot{N}_{100}\rangle$, and give rise to
the processes of excitation, stripping and stimulated capture respectively.

We observe that, in practice, the relevant terms to consider in the RHS of
Eqs. (83) and (84) are those that _(i)_ conserve energy and _(ii)_ change the
number of particles in the state (in this case 100) in the corresponding
diagram. We can use this information to determine $c^{(1)}_{nlm}$ for $nlm\neq
100$ and $c^{(1)}_{\mathbf{k}}$, whose relevant diagrams are shown in Figure
16 (left and right, respectively). The result for $c^{(1)}_{nlm}$ is

$\displaystyle c_{nlm}^{(1)}(t)$
$\displaystyle=-ig\int[dk_{1}][dk_{2}][dk_{3}]a^{*}(\mathbf{k}_{3})a(\mathbf{k}_{1})a(\mathbf{k}_{2})\mathcal{M}_{{k}_{1}+{k}_{2}\rightarrow{k}_{3}+nlm}\int_{0}^{t}e^{i(\omega_{n}+\omega_{k_{3}}-\omega_{k_{1}}-\omega_{k_{2}})t^{\prime}}dt^{\prime}\,$
$\displaystyle-2ig\sqrt{N_{100}^{(0)}}\int[dk_{1}][dk_{2}]a(\mathbf{k}_{1})a^{*}(\mathbf{k}_{2})\mathcal{M}_{{100}+{k}_{1}\rightarrow{k}_{2}+nlm}\int_{0}^{t}e^{i(\omega_{n}+\omega_{k_{2}}-\omega_{k_{1}}-\omega_{1})t^{\prime}}\,,$
(86)

where the first term represents the ‘capture’ process in diagram (1) of Figure
16 (left) and is completely analogous to that of the $nlm=100$ case. The
second represents the ‘excitation’ process in diagram (2) where a particle in
the $100$ state gets excited to a higher state ($n>1$) because of the
interactions with the waves. This is absent for $nlm=100$ and can occur
because the $100$ state is initially populated and
$|\omega_{n}|<|\omega_{1}|$.

The result for $c^{(1)}_{\mathbf{k}}$ is

$\displaystyle c_{\mathbf{k}}^{(1)}(t)$
$\displaystyle=-ig\int[dk_{1}][dk_{2}][dk_{3}]a^{*}(\mathbf{k}_{3})a(\mathbf{k}_{1})a(\mathbf{k}_{2})\mathcal{M}_{{k}_{1}+{k}_{2}\rightarrow{k}_{3}+k}\int_{0}^{t}e^{i(\omega_{k}+\omega_{k_{3}}-\omega_{k_{1}}-\omega_{k_{2}})t^{\prime}}dt^{\prime}\,$
$\displaystyle-
ig\sqrt{N_{100}^{(0)}}\int[dk_{1}][dk_{2}]\big{\\{}a(\mathbf{k}_{1})a(\mathbf{k}_{2})\mathcal{M}_{{k}_{1}+{k}_{2}\rightarrow{k}+100}\int_{0}^{t}e^{i(\omega_{1}+\omega_{k}-\omega_{k_{1}}-\omega_{k_{2}})t^{\prime}}dt^{\prime}$
$\displaystyle+2a(\mathbf{k}_{1})a^{*}(\mathbf{k}_{2})\mathcal{M}_{{k}_{1}+100\rightarrow{k}+{k}_{2}}\int_{0}^{t}e^{i(\omega_{k}+\omega_{k_{2}}-\omega_{k_{1}}-\omega_{1})t^{\prime}}dt^{\prime}\big{\\}}\,,$
(87)

where the first term represents the pure scattering of waves in diagram (1) of
Figure 16 (right), while the second and third term represent the diagrams (3)
and (2). These – as we will see next – are at the origin of the stimulated
capture and stripping processes discussed in the main text.

Eqs. (85), (86) and (87) provide solution of the EoM at the first order in
$g_{\rm eff}$ for the initial condition in Eq. (81), containing only particles
in the scattering states and in the ground state. Note that the results above
are also valid if the initial conditions contain particles _only_ at one
particular level $\tilde{n}\tilde{l}\tilde{m}$, i.e.

$\psi^{(0)}=\sqrt{N_{\tilde{n}\tilde{l}\tilde{m}}^{(0)}}e^{-i\omega_{{\tilde{n}}}t}\psi_{\tilde{n}\tilde{l}\tilde{m}}+\int[dk]a({\mathbf{k}})e^{-i\omega_{k}t}\psi_{\mathbf{k}}\,,$
(88)

in which case it is sufficient to substitute
$100\rightarrow\tilde{n}\tilde{l}\tilde{m}$ in Eqs. (85), (86) and (87). In
Appendix C we will further generalize this derivation for a system with _both_
$100$ and $200$ levels populated, i.e. $N_{100}^{(0)}\neq 0$ and
$N_{200}^{(0)}\neq 0$. Although in principle doable, we do not attempt to
write the general solution for $N_{nlm}^{(0)}\neq 0$ in this work.

#### Occupation number and second-order solution

As mentioned in Section 4.1, the contribution at order $g_{\rm eff}$ to
$\langle{N}_{100}\rangle$ (the second term in Eq. (79)) vanishes, since
$\langle c_{100}^{(1)}\rangle\propto\langle
a(\mathbf{k}_{1})a(\mathbf{k}_{2})a^{*}(\mathbf{k}_{3})\rangle=0$, and the
first nontrivial correction to $\langle N_{100}\rangle$ comes at order $g_{\rm
eff}^{2}$ and reads

$\langle|c_{100}^{(1)}|^{2}\rangle+2\sqrt{N_{100}^{(0)}}{\rm Re}[\langle
c_{100}^{(2)}\rangle]\,.$ (89)

As discussed in Section 4.1, the first term of Eq. (89) is simply the square
of Eq. (85) and can be readily calculated using the expression of $\langle
a^{*}(\mathbf{k}_{1})a(\mathbf{k}_{2})\rangle$ in Eq. (24), along with the
identity $\lim\limits_{\omega
t\rightarrow\infty}|\int_{0}^{t}dt^{\prime}e^{i\omega
t^{\prime}}|^{2}=\lim\limits_{\omega
t\rightarrow\infty}4\omega^{-2}\sin^{2}(\frac{\omega
t}{2})=2\pi\delta(\omega)t$. We obtain

$\langle|c_{100}^{(1)}|^{2}\rangle=2g^{2}\int[dk_{1}][dk_{2}][dk_{3}]|\mathcal{M}_{{k}_{1}+{k}_{2}\rightarrow
100+{k}_{3}}|^{2}f(\mathbf{k}_{1})f(\mathbf{k}_{2})f(\mathbf{k}_{3})(2\pi)\delta(\Delta\omega)t\,,$
(90)

with
$\Delta\omega\equiv\omega_{k_{1}}+\omega_{k_{2}}-\omega_{k_{3}}-\omega_{1}$.
On the other hand, the second term in Eq. (89) requires the solution of the
EoM at the second order, which we report here for simplicity:

$\left(i\partial_{t}+\frac{\nabla^{2}}{2m}+\frac{\alpha}{r}\right)\psi^{(2)}=g\left(\psi^{(0)2}\psi^{(1)*}+2|\psi^{(0)}|^{2}\psi^{(1)}\right)\,.$
(91)

Similarly to Eq. (80), this is a source equation for $\psi^{(2)}$, generated
by a combination of $\psi^{(0)}$ and $\psi^{(1)}$. By plugging the explicit
expression of $\psi^{(2)}$ in terms of $c_{nlm}^{(2)}$ and
$c_{\mathbf{k}}^{(2)}$ and using again the orthonormality relations, we obtain
(as in Eq. (85))

$i\dot{c}_{100}^{(2)}(t)+\dots=ge^{i\omega_{1}t}\int
d^{3}x\left(\psi^{(0)2}\psi^{(1)*}+2|\psi^{(0)}|^{2}\psi^{(1)}\right)\psi_{100}^{*}\,,$
(92)

where the dots stand for frequency corrections analogous to those in Eq. (85).
It is via $\psi^{(1)}$ in the RHS of Eq. (92) that the stimulated capture and
stripping diagrams (2) and (3) of Figure 16 (right), and the excitation
process in diagram (2) of Figure 16 (left), affect
$\langle\dot{N}_{100}\rangle$.

We can substitute into Eq. (92) the expressions of $\psi^{(0)}$ in Eq. (81)
and the solution of $\psi^{(1)}$ computed in Eqs. (85), (86) and (87).
However, as for the first-order coefficients, only the terms corresponding to
processes that conserve energy and change the particle number in the 100 state
contribute to a change in $N_{100}$. In particular:

  * •

Only the unbound component of $\psi^{(0)}$, i.e.
$\int[dk]a(\mathbf{k})e^{-i\omega_{\mathbf{k}}t}\psi_{\mathbf{k}}\subset\psi^{(0)}$,
gives rise to processes that modify the particle number.

  * •

Additionally, only the unbound component of $\psi^{(1)}$, i.e.
$\int[dk]c_{\mathbf{k}}^{(1)}(t)e^{-i\omega_{k}t}\psi_{\mathbf{k}}\subset\psi^{(1)}$,
can conserve energy in the first term in the RHS of Eq. (92) (this is the
reason why we considered only the unbound component of the first order field
in the diagrams of Eq. (39) in the main text). In the second term of Eq. (92)
there is also a contribution from its bound component with $nlm\neq 100$, i.e.
Eq. (86), which gives rise to an excitation term not reported in the main text
(see discussion below).

Let us first consider the contributions arising from the unbound component of
$\psi^{(1)}$. The first term of $c_{\mathbf{k}}^{(1)}$ in Eq. (87) does not
contribute to $\langle\dot{c}_{100}^{(2)}\rangle$, because this term contains
three powers of $a$ and always appears in the RHS of Eq. (92) in combination
with two additional powers of $a$ from $\psi^{(0)}$, and $\langle
f^{5}\rangle=0$ as discussed in Section 3.1. As a result, the process of pure
wave scattering in diagram (1) of Figure 16 (right) does not contribute at the
leading order to the increase of bound particles. On the other hand, the other
two terms in $c_{\mathbf{k}}^{(1)}$ (diagrams (2) and (3) of Figure 16
(right)), corresponding to the stripping-like and stimulated capture-like
diagrams, are relevant. Once substituted into Eq. (92) they lead to

$\displaystyle c_{100}^{(2)}(t)$
$\displaystyle\supset-g^{2}\sqrt{N_{100}^{(0)}}\bigg{\\{}2\int[dk][dp_{1}][dp_{2}]a(\mathbf{p}_{1})a^{*}(\mathbf{p}_{2})\mathcal{M}_{{p}_{1}+{k}\rightarrow
100+{p}_{2}}\int_{0}^{t}dt^{\prime}e^{i(\omega_{1}+\omega_{p_{2}}-\omega_{p_{1}}-\omega_{k})t^{\prime}}$
$\displaystyle\int[dk_{1}][dk_{2}]\big{(}a(\mathbf{k}_{1})a(\mathbf{k}_{2})\mathcal{M}_{{k}_{1}+{k}_{2}\rightarrow
100+{k}}\int_{0}^{t^{\prime}}e^{i(\omega_{1}+\omega_{k}-\omega_{k_{1}}-\omega_{k_{2}})t^{\prime\prime}}dt^{\prime\prime}$
$\displaystyle+2a(\mathbf{k}_{1})a^{*}(\mathbf{k}_{2})\mathcal{M}_{{k}_{1}+100\rightarrow{k}+{k}_{2}}\int_{0}^{t^{\prime}}e^{i(\omega_{k}+\omega_{k_{2}}-\omega_{1}-\omega_{k_{1}})t^{\prime\prime}}dt^{\prime\prime}\big{)}$
$\displaystyle-\int[dk][dp_{1}][dp_{2}]a(\mathbf{p}_{1})a(\mathbf{p}_{2})\mathcal{M}_{{p}_{1}+{p}_{2}\rightarrow
100+{k}}\int_{0}^{t}dt^{\prime}e^{i(\omega_{1}+\omega_{k}-\omega_{p_{1}}-\omega_{p_{2}})t^{\prime}}$
$\displaystyle\int[dk_{1}][dk_{2}]\big{(}a^{*}(\mathbf{k}_{1})a^{*}(\mathbf{k}_{2})\mathcal{M}_{100+{k}\rightarrow{k}_{1}+{k}_{2}}\int_{0}^{t^{\prime}}e^{-i(\omega_{1}+\omega_{k}-\omega_{k_{1}}-\omega_{k_{2}})t^{\prime\prime}}dt^{\prime\prime}$
$\displaystyle+2a^{*}(\mathbf{k}_{1})a(\mathbf{k}_{2})\mathcal{M}_{{k}+{k}_{2}\rightarrow
100+{k}_{1}}\int_{0}^{t^{\prime}}e^{-i(\omega_{k}+\omega_{k_{2}}-\omega_{1}-\omega_{k_{1}})t^{\prime\prime}}dt^{\prime\prime}\big{)}\bigg{\\}}\,.$
(93)

The first (last) three lines of Eq. (93) arise from the second (first) term in
the RHS of Eq. (92), where the two contributions in the first order field
$\psi^{(1)}$ come from the second and third terms of Eq. (87) respectively. We
can represent the four addends of Eq. (93) via the diagrams (1-4) in Figure
17. As in the main text, one of the incoming/outgoing legs is substituted with
the first order field corresponding to the diagrams (2-3) of Figure 16
(right), with momenta labelled by $\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}$.
Note that in Eq. (39) we only reported the diagrams (2) and (3), which are the
only that do not vanish in the determination of $\langle\dot{N}_{100}\rangle$,
see below.

$\mathbf{k}_{1}$$\mathbf{k}_{2}$$100$$\mathbf{k}$$\mathbf{p}_{2}$$100$$\mathbf{p}_{1}$$-i$$-i$$\mathbf{k}_{1}$$100$$\mathbf{k}_{2}$$\mathbf{k}$$\mathbf{p}_{2}$$100$$\mathbf{p}_{1}$$-i$$-i$$\mathbf{p}_{1}$$\mathbf{p}_{2}$$100$$\mathbf{k}$$\mathbf{k}_{1}$$\mathbf{k}_{2}$$100$$-i$$i$$\mathbf{p}_{1}$$\mathbf{p}_{2}$$\mathbf{k}_{2}$$\mathbf{k}$$\mathbf{k}_{1}$$100$$100$$-i$$i$
Figure 17:  The diagrams representing each term of $c_{100}^{(2)}$ in Eq.
(93). Only (2) and (3) contribute to $\langle\dot{N}_{100}\rangle$.

Using the identity

$\langle
a(\mathbf{k}_{1})a(\mathbf{k}_{2})a^{*}(\mathbf{k}_{3})a^{*}(\mathbf{k}_{4})\rangle=(2\pi)^{6}f(\mathbf{k}_{1})f(\mathbf{k}_{2})\left[\delta(\mathbf{k}_{1}-\mathbf{k}_{3})\delta(\mathbf{k}_{2}-\mathbf{k}_{4})+\delta(\mathbf{k}_{1}-\mathbf{k}_{4})\delta(\mathbf{k}_{2}-\mathbf{k}_{3})\right]\,,$

which follows from Eq. (24) and Wick’s theorem, Eq. (93) simplifies to

$\displaystyle\langle c_{100}^{(2)}(t)\rangle$
$\displaystyle\supset-4g^{2}\sqrt{N_{100}^{(0)}}\int[dk][dk_{1}][dk_{2}]f(\mathbf{k}_{1})f(\mathbf{k}_{2})|\mathcal{M}_{{k}_{1}+{k}_{2}\rightarrow
100+{k}}|^{2}\int_{0}^{t}dt^{\prime}\int_{0}^{t^{\prime}}dt^{\prime\prime}e^{i(\omega_{1}+\omega_{k_{1}}-\omega_{k_{2}}-\omega_{k})(t^{\prime}-t^{\prime\prime})}$
$\displaystyle+2g^{2}\sqrt{N_{100}^{(0)}}\int[dk][dk_{1}][dk_{2}]f(\mathbf{k}_{1})f(\mathbf{k}_{2})|\mathcal{M}_{100+{k}\rightarrow{k}_{1}+{k}_{2}}|^{2}\int_{0}^{t}dt^{\prime}\int_{0}^{t^{\prime}}dt^{\prime\prime}e^{i(\omega_{1}+\omega-\omega_{k_{1}}-\omega_{k_{2}})(t^{\prime}-t^{\prime\prime})}\,.$
(94)

The time integration can be performed explicitly using
$\int_{0}^{t}dt^{\prime}\int_{0}^{t^{\prime}}dt^{\prime\prime}e^{i\omega(t^{\prime}-t^{\prime\prime})}+c.c=4\omega^{-2}\sin^{2}(\omega
t)\rightarrow 2\pi\delta(\omega)t$, and yields

$\displaystyle\langle c_{100}^{(2)}+c_{100}^{(2)*}\rangle$
$\displaystyle\supset
2g^{2}\sqrt{N_{100}^{(0)}}\int[dk_{1}][dk_{2}][dk_{3}]|\mathcal{M}_{{k}_{1}+{k}_{2}\rightarrow
100+{k}_{3}}|^{2}[f(\mathbf{k}_{1})f(\mathbf{k}_{2})-2f(\mathbf{k}_{1})f(\mathbf{k}_{3})](2\pi)\delta(\Delta\omega)t\,.$
(95)

This equation provides the capture and stripping terms in Eq. (41). As
mentioned in Section 4.1, this corresponds to the diagrams in Eq. (40) of the
main text, which can be thought as arising from the ‘contraction’ of the
diagrams (2) and (3) in Figure 17.

Similarly, only the second term in the bound component $c^{(1)}_{nlm}$ in Eq.
(86), representing the excitation process, contributes to $c_{100}^{(2)}$;
this term gives

$\displaystyle\langle c_{100}^{(2)}+c_{100}^{(2)*}\rangle$
$\displaystyle\supset
2g^{2}\sqrt{N_{100}^{(0)}}\int[dk_{1}][dk_{2}]\sum_{nlm}|\mathcal{M}_{nlm+{k}_{1}\rightarrow
100+{k}_{2}}|^{2}[-2f(\mathbf{k}_{1})f(\mathbf{k}_{2})]2\pi\delta(\omega_{n}+\omega_{k_{1}}-\omega_{1}-\omega_{k_{2}})t\,.$
(96)

Plugging the last two equations into Eq. (89) we obtain the full result:

$\displaystyle\langle\dot{N}_{100}$
$\displaystyle\rangle=\frac{d}{dt}\left\langle|c_{100}^{(1)}|^{2}+\sqrt{N_{100}^{(0)}}(c_{100}^{(2)}+c_{100}^{(2)*})\right\rangle$
$\displaystyle=2g^{2}\int[dk_{1}][dk_{2}][dk_{3}]|\mathcal{M}_{{k}_{1}+{k}_{2}\rightarrow{k}_{3}+100}|^{2}(2\pi)\delta(\Delta\omega)\\{f(\mathbf{k}_{1})f(\mathbf{k}_{2})f(\mathbf{k}_{3})+N_{100}^{(0)}[f(\mathbf{k}_{1})f(\mathbf{k}_{2})-2f(\mathbf{k}_{2})f(\mathbf{k}_{3})]\\}\,$
$\displaystyle+2g^{2}N_{100}^{(0)}\int[dk_{1}][dk_{2}]\sum_{nlm}|\mathcal{M}_{{nlm}+k_{1}\rightarrow
100+k_{2}}|^{2}(2\pi)\delta(\omega_{n}+{\omega}_{k_{1}}-{\omega}_{k_{2}}-{\omega}_{1})\left[-2f(\mathbf{k}_{1})f(\mathbf{k}_{2})\right]\,.$
(97)

This equation coincides with Eq. (41) up to the term in the last line,
neglected in the main text. This describes the depletion of particles from the
100 level because of their excitation to all higher levels via scatterings
with the DM waves, as depicted in diagram (2) in Figure 16 (left). Note the
similarity between the stripping and excitation terms, which have an analogous
form provided $\int[dk]$ is substituted with $\sum_{nlm}$ (this is because the
diagrams (2, left) and (2, right) in Figure 16 are identical upon substituting
$\mathbf{k}\leftrightarrow nlm$).

The excitation term could in principle reduce the number of particles in the
$100$ state (when the captured mass is large enough, since it is proportional
to $N_{100}^{(0)}$). However, at the same time, the occupation number of
excited states would grow. As we will show in Appendix C, for a system
consisting of $100$ and $200$ states, the excited state (if populated) decays
(‘relaxes’ or ‘de-excites’) to the ground state, at least in the regime
$\xi_{\rm foc}=2\pi\alpha/v_{\rm dm}\gg 1$. Consequently, $N_{100}$ increases.
As a result, the excitation terms are not expected to affect the discussion in
the main text.

##  Appendix C Role of excited states and two-level system

In the main text we studied the capture of ULDM to the ground state, but
neglected excited states in most of our discussions. In this Appendix, we show
that these are likely to only enhance the accretion of DM onto the ground
state, and not to change the overall conclusions of Section 4.

As argued in Section 4.3, the matrix element
$\mathcal{M}\equiv\mathcal{M}_{k_{1}+k_{2}\rightarrow k_{3}+nlm}$ in Eq. (31)
is suppressed for $n\gg 1$, in both limits of $\xi_{\rm foc}\gg 1$ and
$\xi_{\rm foc}\ll 1$, which inhibits the population of very high-$n$ states.
However, we also argued that for $\xi_{\rm foc}=2\pi\alpha/v_{\rm dm}\gg 1$
the magnitude of $\mathcal{M}$ for the first few levels (those satisfying
$n\lesssim\xi_{\rm foc}^{1/2}$) could be comparable to that of $nlm=100$ (this
can be seen, for instance for $nlm=200$, by using the explicit expression of
$\psi_{nlm}$ in Eq. (70)). Such levels are therefore populated with a similar
rate as the 100 state, and can affect the density of the gravitational atom.
We now study the evolution of a system where both the 100 and 200 levels are
included (and populated by capture from the DM waves) and show that: For
$\xi_{\rm foc}\ll 1$ the population of the excited state $200$ is comparable
to $100$, but with significantly reduced density due to its large radius; and
in the limit $\xi_{\rm foc}\gg 1$, the excited state is rapidly populated but
transitions to the ground state over a timescale similar to the exponential
increase time, $\tau_{\rm rel}$.

Let us first derive a generalized version of Eq. (41) that describes the
simultaneous evolution of the occupation number of both $100$ and $200$
states, and comment later about the inclusion of all other levels. To do so,
we solve the EoM in Eq. (72) perturbatively, starting from the initial
condition:

$\psi^{(0)}=\sqrt{N_{100}^{(0)}}e^{-i\omega_{1}t}\psi_{100}+\sqrt{N_{200}^{(0)}}e^{-i\omega_{2}t}\psi_{200}+\int[dk]a({\mathbf{k}})e^{-i\omega_{k}t}\psi_{\mathbf{k}}\,.$
(98)

Contrary to Eq. (81), this has both $N_{100}^{(0)}\neq 0$ and
$N_{200}^{(0)}\neq 0$. Of course, $N_{100}^{(0)}=N_{200}^{(0)}=0$ at $t=0$ in
our Solar System, but, as in Section 4.1, keeping them generic will allow us
to study the rates $\langle\dot{N}_{100}\rangle$ and
$\langle\dot{N}_{200}\rangle$ at an arbitrary time after such levels get
populated from the DM capture.

As in Appendix B, the occupation number $\langle{N}_{100}\rangle$ is given in
Eq. (89), and requires the calculation of $c_{100}^{(1)}$ and $c_{100}^{(2)}$.
By plugging Eq. (98) into Eq. (83) and selecting only the terms that conserve
energy and change the particle number in the 100 state, we find

$\displaystyle c_{100}^{(1)}(t)$
$\displaystyle=-ig\bigg{[}\int[dk_{1}][dk_{2}][dk_{3}]a(\mathbf{k}_{1})a(\mathbf{k}_{2})a^{*}(\mathbf{k}_{3})\mathcal{M}_{{k}_{1}+{k}_{2}\rightarrow{k}_{3}+100}\int_{0}^{t}e^{i(\omega_{1}+\omega_{k_{3}}-\omega_{k_{1}}-\omega_{k_{2}})t^{\prime}}dt^{\prime}$
$\displaystyle+2\sqrt{N_{200}^{(0)}}\int[dk_{1}][dk_{2}]a(\mathbf{k}_{1})a^{*}(\mathbf{k}_{2})\mathcal{M}_{{k}_{1}+200\rightarrow{k}_{2}+100}\int_{0}^{t}e^{i(\omega_{1}+\omega_{k_{2}}-\omega_{k_{1}}-\omega_{2})t^{\prime}}dt^{\prime}$
$\displaystyle+N_{200}^{(0)}\int[dk]a(\mathbf{k})\mathcal{M}_{200+200\rightarrow{k}+100}\int_{0}^{t}e^{i(\omega_{1}+\omega_{k}-2\omega_{2})t^{\prime}}dt^{\prime}\bigg{]}\,.$
(99)

The $\langle|c_{100}^{(1)}|^{2}\rangle$ part of $\langle N_{100}\rangle$ is
simply given by the square of the amplitude of each term in Eq. (99) (the
cross terms vanish upon average) and reads

$\displaystyle\langle|c_{100}^{(1)}|^{2}\rangle=$
$\displaystyle\,g^{2}t\Big{[}2\int[dk_{1}][dk_{2}][dk_{3}]|\mathcal{M}_{{k}_{1}+{k}_{2}\rightarrow
100+{k}_{3}}|^{2}(2\pi)\delta(\omega_{1}+\omega_{k_{3}}-\omega_{k_{1}}-\omega_{k_{2}})f(\mathbf{k}_{1})f(\mathbf{k}_{2})f(\mathbf{k}_{3})$
$\displaystyle+4N_{200}^{(0)}\int[dk_{1}][dk_{2}]|\mathcal{M}_{{k}_{1}+200\rightarrow{k}_{2}+100}|^{2}(2\pi)\delta(\omega_{1}+\omega_{k_{2}}-\omega_{k_{1}}-\omega_{2})f(\mathbf{k}_{1})f(\mathbf{k}_{2})$
$\displaystyle+2N_{200}^{(0)2}\int[dk]|\mathcal{M}_{200+200\rightarrow{k}+100}|^{2}(2\pi)\delta(\omega_{1}+\omega_{k}-2\omega_{2})f(\mathbf{k})\Big{]}\,.$
(100)

The first term in Eqs. (99) and (100) describes the (already discussed) direct
capture from the DM waves, and is shown in diagram (1a) in Figure 18, where
the conventions for the lines are as in Section 4.1. In addition to this, two
more terms appeared: They describe the increase in the occupation number from
the processes ${k}_{1}+200\rightarrow 100+{k}_{2}$ and
$200+200\rightarrow{k}+100$. The former, dubbed ‘de-excitation’, is shown in
diagram (1b); the latter, dubbed ‘relaxation’, in diagram (1c). These are only
active when $N_{200}^{(0)}\neq 0$, and increase the occupation number in the
ground state at the expenses of the populated 200 state.

$\mathbf{k}_{1}$$\mathbf{k}_{2}$$\mathbf{k}_{3}$$100$capture(1a)$\mathbf{k}_{1}$$\mathbf{k}_{2}$$100$$\mathbf{k}$stimulated
capture(2a)$100$$\mathbf{k}_{3}$$\mathbf{k}_{2}$$\mathbf{k}$stripping(3a)

$200$$\mathbf{k}_{2}$$\mathbf{k}_{1}$$100$de-
excitation(1b)$100$$\mathbf{k}_{1}$$\mathbf{k}_{2}$$200$excitation(2b)$200$$\mathbf{k}_{2}$$100$$\mathbf{k}_{1}$stimulated
de-excitation(3b)$100$$\mathbf{k}_{1}$$200$$\mathbf{k}_{2}$stimulated
excitation(4b)

$200$$200$$\mathbf{k}$$100$relaxation(1c)$200$$200$$100$$\mathbf{k}$stimulated
relaxation(2c)$100$$\mathbf{k}$$200$$200$de-relaxation(3c)

Figure 18: A schematic representation of the processes that contribute to the
change in the number of bound particles $\langle\dot{N}_{100}\rangle$, see Eq.
(C), in a system where both the 100 and 200 levels are present. As in Figure
8, double solid lines represents the occupation number factors in the
different terms of Eq. (C).

As in Appendix B, one could proceed to calculate the second order coefficient
$c_{100}^{(2)}$, and thus the second part of $\langle N_{100}\rangle$. The
only difference is that in $\psi^{(0)}$ on the RHS of the second order EoM,
Eq. (91), both the 200 state and the unbound component contribute. As before,
the latter provides the stimulated capture and stripping terms in Eqs.
(95-96), represented in the diagrams (2a-3a) in Figure 18, while the former
provides new terms that involve $N_{200}^{(0)}$.

Rather than going through this lengthy derivation, we can obtain
(heuristically) the form of the terms in $\langle\dot{N}_{100}\rangle$ that
depend on $N_{200}^{(0)}$ as we did in Section 4.2, i.e. by demanding that the
dependence on the occupation number of the two new processes
${k}_{1}+200\leftrightarrow 100+{k}_{2}$ and $200+200\leftrightarrow{k}+100$
(which include direct and inverse processes) is consistent with Bose
enhancement. In other words, the contribution to $\langle\dot{N}_{100}\rangle$
from the former should be proportional to

$\displaystyle
f(\mathbf{k}_{1})N_{200}[f(\mathbf{k}_{2})+1][N_{100}+1]-[f(\mathbf{k}_{1})+1][N_{200}+1]f(\mathbf{k}_{2})N_{100}$
$\displaystyle\simeq\ $ $\displaystyle
f(\mathbf{k}_{1})f(\mathbf{k}_{2})(N_{200}-N_{100})+N_{100}N_{200}[f(\mathbf{k}_{2})-f(\mathbf{k}_{1})]\,.$
(101)

while for the latter to

$\displaystyle
N_{200}^{2}[f(\mathbf{k})+1][N_{100}+1]-[N_{200}+1]^{2}f(\mathbf{k})N_{100}$
$\displaystyle\simeq\ $ $\displaystyle
f(\mathbf{k})N_{200}^{2}+N_{100}N_{200}^{2}-2N_{100}N_{200}f(\mathbf{k})\,,$
(102)

where the approximations are valid for $N_{nlm}\gg 1$ (for simplicity of
notation, here and in the following we omit the superscript ‘$(0)$’). Then,
the evolution of the occupation number $\langle\dot{N}_{100}\rangle$ in the
ground state can be obtained by starting from Eq. (100) and adding more terms
in the RHS in such a way that the dependence on the occupation numbers
$N_{100}$ and $N_{200}$ matches that of the two previous equations. These
terms are simply the permutations (2b-4b) and (2c-3c) of the diagrams (1b) and
(1c) respectively, weighted with the appropriate factors. The full result is

$\displaystyle\dot{N}_{100}$
$\displaystyle=2g^{2}\int[dk_{1}][dk_{2}][dk_{3}]|\mathcal{M}_{{k}_{1}+{k}_{2}\rightarrow{k}_{3}+100}|^{2}(2\pi)\delta(\omega_{k_{1}}+\omega_{k_{2}}-\omega_{1}-\omega_{k_{3}})$
$\displaystyle\left\\{f(\mathbf{k}_{1})f(\mathbf{k}_{2})f(\mathbf{k}_{3})+N_{100}[f(\mathbf{k}_{1})f(\mathbf{k}_{2})-2f(\mathbf{k}_{2})f(\mathbf{k}_{3})]\right\\}$
$\displaystyle+4g^{2}\int[dk_{1}][dk_{2}]|\mathcal{M}_{{k}_{1}+100\rightarrow{k}_{2}+200}|^{2}(2\pi)\delta(\omega_{k_{1}}+\omega_{1}-\omega_{2}-\omega_{k_{2}})$
$\displaystyle\\{N_{100}N_{200}[f(\mathbf{k}_{2})-f(\mathbf{k}_{1})]-f(\mathbf{k}_{1})f(\mathbf{k}_{2})[N_{100}-N_{200}]\\}$
$\displaystyle+2g^{2}\int[dk]|\mathcal{M}_{200+200\rightarrow
100+{k}}|^{2}(2\pi)\delta(2\omega_{2}-\omega_{1}-\omega_{k})\\{f(\mathbf{k})N_{200}^{2}+N_{100}N_{200}^{2}-2N_{200}N_{100}f(\mathbf{k})\\}\,.$
(103)

One can similarly obtain the evolution of $\langle N_{200}\rangle$, which
reads

$\displaystyle\dot{N}_{200}$
$\displaystyle=2g^{2}\int[dk_{1}][dk_{2}][dk_{3}]|\mathcal{M}_{{k}_{1}+{k}_{2}\rightarrow{k}_{3}+200}|^{2}(2\pi)\delta(\omega_{k_{1}}+\omega_{k_{2}}-\omega_{2}-\omega_{k_{3}})$
$\displaystyle\\{f(\mathbf{k}_{1})f(\mathbf{k}_{2})f(\mathbf{k}_{3})+N_{200}[f(\mathbf{k}_{1})f(\mathbf{k}_{2})-2f(\mathbf{k}_{2})f(\mathbf{k}_{3})]\\}$
$\displaystyle-4g^{2}\int[dk_{1}][dk_{2}]|\mathcal{M}_{{k}_{1}+100\rightarrow{k}_{2}+200}|^{2}(2\pi)\delta(\omega_{k_{1}}+\omega_{1}-\omega_{2}-\omega_{k_{2}})$
$\displaystyle\\{N_{100}N_{200}[f(\mathbf{k}_{2})-f(\mathbf{k}_{1})]-f(\mathbf{k}_{1})f(\mathbf{k}_{2})[N_{100}-N_{200}]\\}$
$\displaystyle-4g^{2}\int[dk]|\mathcal{M}_{200+200\rightarrow
100+{k}}|^{2}(2\pi)\delta(2\omega_{2}-\omega_{1}-\omega_{k})\\{f(\mathbf{k})N_{200}^{2}+N_{100}N_{200}^{2}-2N_{200}N_{100}f(\mathbf{k})\\}\,,$
(104)

and is analogous to Eq. (C), provided that $N_{100}$ and $N_{200}$ are
exchanged, appropriate minus signs added and a factor of 2 in the last line of
Eq. (C). Note that:

  * •

The first two lines of Eq. (C) reproduce the capture, stimulated capture and
stripping processes, diagrams (1a-3a) of Figure 18.

  * •

The third and fourth lines include the de-excitation process (1b) and its
inverse (‘excitation’, (2b)), as well as the same processes that are
‘stimulated’ by the 200 state being already populated (3b-4b); these appear
thanks to Bose enhancement. Excitations deplete the 100 state given their
minus sign. Observe that the excitation term comes together with a tower of
additional excitations to higher and higher states, similar to the last term
in Eq. (B). We omitted these because not crucial for what follows.

  * •

Finally, the last line includes the relaxation process (1c), as well as its
stimulated version (2c) and its inverse process (‘de-relaxation’, 3c).

Eqs. (C-C) constitute a closed system that determines the evolution of
$N_{100}(t)$ and $N_{200}(t)$. However, as discussed at the beginning of this
Section, higher levels (or levels with nontrivial angular momentum) could be
also populated from DM capture. Thus, in principle, one should derive a more
general set of coupled equations describing also the evolution e.g. of the
$300$ level, by including $N_{300}^{(0)}\neq 0$ in the initial conditions of
Eq. (98). In this way one could study also the population of these levels, and
track their evolution afterwards. Eventually this leads to an infinite (and
intractable) set of coupled equations (similar in form to Eqs. (C-C))
describing the simultaneous evolution of the occupation number of all levels.
However, in the following we will prove that:

  * •

In the limit $2\pi\alpha/v_{\rm dm}\ll 1$, the two-level system in Eqs. (C-C)
tends to an equilibrium configuration where $N_{100}$ and (later) $N_{200}$
have a common value, which turns out to coincide with that obtained from the
single-level equation, Eq. (41) in the main text. Thus, the density at the
center is dominated by the ground state.

  * •

In the limit $2\pi\alpha/v_{\rm dm}\gg 1$ the excited state, initially
exponentially enhanced, quickly decays to the ground state.

We expect this dynamics to be prototypical of the all-level system: Excited
states quickly decay to lower-level states in the second limit, and the
density is well reproduced by the single level system in the first; as a
result, it is likely that the full system behaves, to a good approximation, as
the ground state.

#### Dilute gravitational atoms: $\xi_{\rm foc}=2\pi\alpha/v_{\rm dm}\ll 1$

We first consider the limit $2\pi\alpha/v_{\rm dm}\ll 1$. As mentioned in
Section 4.3, $\mathcal{M}_{{k}_{1}+{k}_{2}\rightarrow{k}_{3}+n00}\propto
n^{-3/2}$; thus, the capture and stripping rates for the $200$ state are
approximately $1/8$ times those for $100$. Energy conservation,
$\omega_{k_{1}}+\omega_{1}-\omega_{2}-\omega_{k_{2}}=0$, forces
${k}_{1}\simeq{k}_{2}$ in the $2\pi\alpha/v_{\rm dm}\ll 1$ limit. As a result,
since the sum of the excitation and de-excitation terms is proportional to
$f(\mathbf{k}_{2})-f(\mathbf{k}_{1})$, these cancel between each other.
Consequently, Eqs. (C) and (C) can be approximated as

$\displaystyle\dot{N}_{100}$
$\displaystyle=a_{1}-a_{2}N_{100}-a_{3}(N_{100}-N_{200})+a_{4}(N_{200}^{2}-2N_{200}N_{100})+a_{5}N_{100}N_{200}^{2}\,,$
$\displaystyle\dot{N}_{200}$
$\displaystyle=a_{1}/8-(a_{2}/8)N_{200}+a_{3}(N_{100}-N_{200})-2a_{4}(N_{200}^{2}-2N_{200}N_{100})-2a_{5}N_{100}N_{200}^{2}\,.$
(105)

As in Section 4.3, the coefficients $a_{1}$ and $a_{2}$ are

$\displaystyle a_{1}=\frac{16384\sqrt{\pi}Bg^{2}\rho_{\rm
dm}^{3}\alpha^{4}}{m^{7}v_{\rm dm}^{9}}~{},\qquad
a_{2}=\frac{A}{2.4}\frac{(2\pi\alpha/v_{\rm dm})^{4}}{\tau_{\rm rel}}~{},$
(106)

with $A\simeq 0.21$ and $B\simeq 0.08$. The other coefficients may be
evaluated in a similar way:

$\displaystyle a_{3}\simeq 0.04\frac{(2\pi\alpha/v_{\rm dm})^{4}}{\tau_{\rm
rel}}~{},\qquad a_{4}=\frac{10816g^{2}\rho_{\rm
dm}\alpha^{3}m}{4782969\sqrt{\pi}v_{\rm dm}^{2}}f_{4}(\alpha/v_{\rm
dm})~{},\qquad a_{5}=\frac{2704\sqrt{2}g^{2}m^{5}\alpha^{4}}{4782969\pi}~{},$
(107)

where $f_{4}(a)=2\sqrt{2}a/e$ for $a\ll 1$.313131The exact form is
$f_{4}(a)\equiv e^{-(1+a/\sqrt{2})^{2}}(e^{2\sqrt{2}a}-1)$. These equations
are greatly simplified when rewritten in terms of $t^{\prime}=t/\tau_{\rm
rel}$ and of the normalized occupation number

$\tilde{N}_{n00}\equiv\frac{m^{4}\alpha^{3}}{\rho_{\rm dm}}N_{n00}=\pi
n^{3}\frac{\rho_{n00}}{\rho_{\rm dm}}~{},\qquad\rho_{n00}\equiv
mN_{n00}|\psi_{n00}(0)|^{2}=\frac{mN_{n00}}{\pi(nR_{\star})^{3}}\,.$ (108)

Note that $\tilde{N}_{n00}$ represents (modulo a $\pi n^{3}$ factor) the
contribution $\rho_{n00}$ from the state $n00$ to the density of the bound
state at the center $r=0$, normalized to $\rho_{\rm dm}$ (i.e.
$\rho(0)/\rho_{\rm dm}$). In particular, for the single-level system in Eq.
(41) at equilibrium we have $\tilde{N}_{100}\sim\rho_{100}/\rho_{\rm
dm}\sim(2\pi\alpha/v_{\rm dm})^{3}\ll 1$, see Eq. (9). Using these variables,
we obtain

$\displaystyle\frac{d\tilde{N}_{100}}{dt^{\prime}}$ $\displaystyle=2\cdot
10^{3}\left[\frac{\alpha}{v_{\rm
dm}}\right]^{7}-10^{2}\left[\frac{\alpha}{v_{\rm
dm}}\right]^{4}\tilde{N}_{100}-60\left[\frac{\alpha}{v_{\rm
dm}}\right]^{4}\left[\tilde{N}_{100}-\tilde{N}_{200}\right]$
$\displaystyle+10^{-3}f_{4}\left(\frac{\alpha}{v_{\rm
dm}}\right)\left[\tilde{N}_{200}^{2}-2\tilde{N}_{100}\tilde{N}_{200}\right]+2.5\cdot
10^{-4}\left[\frac{v_{\rm
dm}}{\alpha}\right]^{2}\tilde{N}_{100}\tilde{N}_{200}^{2}\,,$ (109)
$\displaystyle\frac{d\tilde{N}_{200}}{dt^{\prime}}$ $\displaystyle=3\cdot
10^{2}\left[\frac{\alpha}{v_{\rm dm}}\right]^{7}-20\left[\frac{\alpha}{v_{\rm
dm}}\right]^{4}\tilde{N}_{200}+60\left[\frac{\alpha}{v_{\rm
dm}}\right]^{4}\left[\tilde{N}_{100}-\tilde{N}_{200}\right]$
$\displaystyle-2\cdot 10^{-3}f_{4}\left(\frac{\alpha}{v_{\rm
dm}}\right)\left[\tilde{N}_{200}^{2}-2\tilde{N}_{100}\tilde{N}_{200}\right]-5\cdot
10^{-4}\left[\frac{v_{\rm
dm}}{\alpha}\right]^{2}\tilde{N}_{100}\tilde{N}_{200}^{2}\,.$ (110)

For $\xi_{\rm foc}=2\pi\alpha/v_{\rm dm}\ll 1$ the equilibrium solution (i.e.
$\dot{N}_{nlm}=0$) of this system occurs at $N_{100}=N_{200}\equiv N$, with
equilibrium value for $N_{100}$ coinciding with that of the single-level
system in Eq. (41). This is because for such a configuration the stimulated
excitation and de-excitation terms (proportional to
$\tilde{N}_{100}-\tilde{N}_{200}$ in the previous equations) cancel between
each other, while all the relaxation terms (in last line) are much smaller
than the capture and stimulated capture ones, which therefore dominate.
Indeed, the relaxation terms are proportional to
$f_{4}\tilde{N}^{2}\propto(\alpha/v_{\rm dm})^{7}$ and $(v_{\rm
dm}/\alpha)^{2}\tilde{N}^{3}\propto(\alpha/v_{\rm dm})^{7}$, but their
coefficient is suppressed by at least four orders of magnitudes compared to
that of the capture and stimulated capture, ultimately because $\tilde{N}\ll
1$.

Figure 19 (left) shows the numerical evolution of the system in Eqs. (109-110)
starting from $N_{100}=N_{200}=0$ for $\xi_{\rm foc}=0.1$. Initially the 200
level grows slower than the ground state due to the $1/8$ factor. However,
when the occupation numbers become large enough, the excitation terms become
important and quickly bring the occupation numbers close to each other. In the
end, both states reach the equilibrium, with occupation numbers similar to the
value in the one level system. As mentioned, since the $200$ level has a
radius twice as larger as that of the ground state, the density near of center
of the gravitational atom is dominated by the ground state.

Figure 19: The time evolution of the normalized occupation numbers
$\tilde{N}_{100}$ (blue) and $\tilde{N}_{200}$ (orange), by numerically
solving the two-level system, Eqs. (109-110), with $\xi_{\rm
foc}=2\pi\alpha/v_{\rm dm}=0.1$ (left), and Eqs. (111-112) with $\xi_{\rm
foc}=10$ (right). $\tilde{N}_{n00}/(\pi n^{3})$ represents the contribution
from the state $n00$ to the central density, normalized to $\rho_{\rm dm}$.
For comparison, we also show the evolution of the occupation number of the
one-level system studied in the main text (green). In the second limit, both
the states increases exponentially, but the excited state quickly transitions
to the ground state soon after $t\simeq\tau_{\rm rel}$.

#### Dense gravitational atoms: $\xi_{\rm foc}=2\pi\alpha/v_{\rm dm}\gg 1$

The only processes in Figure 18 that are not suppressed in the regime
$2\pi\alpha/v_{\rm dm}\gg 1$ are _(i)_ stimulated capture (diagram (2a) of
Figure 18, discussed in the main text), _(ii)_ stimulated de-excitation
(diagram (3b)) and _(iii)_ stimulated relaxation (diagram (2c)). All the other
diagrams are exponentially suppressed because do not benefit from Bose
enhancement. In particular:

  * •

For the excitation and de-excitation processes ${k}_{1}+100\leftrightarrow
200+{k}_{2}$, central panel of Figure 18, energy conservation dictates
$k_{1}=\sqrt{k_{2}^{2}+3m^{2}\alpha^{2}/4}\gg mv_{\rm dm}$, so the diagrams
(1b,2b,4b), proportional to $f(\mathbf{k}_{1})$, are exponentially suppressed.

  * •

Similarly, the relaxation and de-relaxation processes $200+200\leftrightarrow
100+{k}$ in diagrams (1c,3c) are proportional to $f(\mathbf{k})$ and by energy
conservation $k=m\alpha/\sqrt{2}\gg mv_{\rm dm}$, so the resulting integrals
are small.

The unsuppressed processes (2a,3b,2c) only induce level transitions from 200
to 100. As we now show in detail, if populated, the 200 state will be
therefore depleted into the 100 state, and will increase its growth. At $t=0$,
$N_{100}=N_{200}=0$ and the only term that contributes to a change in the
number of bound particles is the (suppressed) direct capture, as all the other
terms are proportional to further powers of $N_{100}$ and $N_{200}$. Ignoring
all the terms that are exponentially suppressed except for direct capture, in
the $2\pi\alpha/v_{\rm dm}\gg 1$ limit Eqs. (C-C) are well approximated by:

$\displaystyle\dot{N}_{100}$
$\displaystyle=C_{100}/m+\Gamma_{100}N_{100}+DN_{100}N_{200}+RN_{100}N_{200}^{2}\,,$
(111) $\displaystyle\dot{N}_{200}$
$\displaystyle=C_{200}/m+\Gamma_{200}N_{200}-DN_{100}N_{200}-2RN_{100}N_{200}^{2}\,,$
(112)

where the coefficients can be calculated similarly to Section 4.3,323232The
integral in the matrix element $\mathcal{M}_{100+k_{1}\rightarrow 200+k_{2}}$
is a function only of one variable, $\hat{k}_{1}\cdot\hat{k}_{3}$, and can be
computed numerically; the corresponding integral in
$\mathcal{M}_{200+200\rightarrow 100+k}$ is computable analytically. and read

$\displaystyle C_{n00}\simeq\frac{16g^{2}\rho_{\rm dm}^{3}}{m^{6}v_{\rm
dm}^{5}}s(\alpha^{2}/n^{2}v_{\rm dm}^{2})\,,\ \ \
n=1,2\,,\qquad\quad\Gamma_{100}\simeq\Gamma_{200}\simeq\frac{1}{0.34\tau_{\rm
rel}}\simeq\frac{3g^{2}\rho_{\rm dm}^{2}}{m^{3}v_{\rm dm}^{2}}\,,$ (113)

where to a good approximation $s(x)\simeq e^{-x+\sqrt{4+2x}-2}$, and

$\displaystyle D\simeq\frac{0.035g^{2}\alpha^{2}m\rho_{\rm dm}}{v_{\rm
dm}}\,,\qquad\quad R\simeq 1.3\cdot 10^{-4}g^{2}\alpha^{4}m^{5}\,.$ (114)

The solution of Eqs. (111-112) with initial conditions $N_{100}=N_{200}=0$ can
be understood qualitatively. First, notice that $\Gamma_{100}$ and
$\Gamma_{200}$ are of the same order (because the matrix elements
$M_{k_{1}+k_{2}\rightarrow k_{3}+n00}$ with $n=1,2$ are similar), so we expect
also $N_{200}$ to increase exponentially, with similar timescale $\tau_{\rm
rel}$. The capture coefficients $C_{100}$ and $C_{200}$ are comparable, but
since $s(x)$ is exponentially decreasing function of $\alpha/(nv_{\rm dm})$
with $n=1,2$, the latter is less suppressed at large $2\pi\alpha/v_{\rm dm}$.
Ultimately, this is because the exponential suppression from
$f(\mathbf{k}_{3})$, with
$k_{3}=\sqrt{k_{1}^{2}+k_{2}^{2}+\alpha^{2}m^{2}/n^{2}}$, is smaller for $n=2$
than $n=1$. Thus, as mentioned in Section 4.3, for $2\pi\alpha/v_{\rm dm}\gg
1$ the mass captured until the time $1/\Gamma$ is larger for the 200 level.

As a result, initially both levels grow; while $N_{100}$ and $N_{200}$ are
small, the growth is driven by the capture terms and is linear. As soon as
$N_{100}$ and $N_{200}$ are large enough, the stimulated capture terms drive
the exponential increase of the levels at $t\simeq\tau_{\rm rel}$. Soon after,
the stimulated de-excitation and stimulated relaxation terms (last and second-
to-last terms in Eq. (112)) become as large as the stimulated capture term,
because they are proportional to higher powers of $N_{nlm}$.333333On the other
hand, the exponentially-suppressed terms of Figure 19, not reported in Eq.
(112), will be suppressed with respect to these, and therefore not considered
in the evolution. By comparing the capture term with the stimulated de-
excitation and the relaxation ones, and using the fact that (before these
become relevant) $N_{100}\propto N_{200}\propto e^{t/\tau_{\rm rel}}$, we can
estimate that the latter start dominating over the former, inducing the decay
of the 200 level, at times $t_{c,D}\simeq\tau_{\rm rel}\log_{D}$ and
$t_{c,R}\simeq 2\tau_{\rm rel}\log_{R}$, respectively. The factors $\log_{D}$
and $\log_{R}$ are of order one and can be estimated from Eq. (112) as
$\log_{D}\simeq\log(15v_{\rm dm}^{2}/f_{2}\alpha^{2})$ and $\log_{R}\simeq
2\log(27v_{\rm dm}^{2}/f_{2}\alpha^{2})$; here $f_{2}$ is the density of the
200 state in units of the DM density at time $t\simeq\tau_{\rm rel}$. We
conclude that, soon after the exponential increase starts, 200 decays into
100. Figure 19 shows the time evolution of $\tilde{N}_{100}$ and
$\tilde{N}_{200}$ for $2\pi\alpha/v_{\rm dm}=10$, from the numerical solution
of Eqs. (111-112) with initial condition $N_{100}=N_{200}=0$, where the main
features described above are reproduced. In particular, the particle number in
the $200$ state becomes negligible at $t\gtrsim 2\tau_{\rm rel}$. The dynamics
of this simplified system should be prototypical of the full case including
additional excited states. It could also possibly describe states with
nontrivial angular momentum, though a more detailed analysis is required.

##  Appendix D Details on Gross–Pitaevskii simulations

In this Appendix we give a brief description of the simulations presented in
Section 5 and provide additional results.

##### Equations of motion and initial conditions.

As mentioned in Section 5, the combination of the EoM in Eq. (15) and initial
conditions in Eq. (23) depends only on the three dimensionless parameters
$\\{\alpha/v_{\rm dm},\tilde{g}\equiv g\rho_{\rm dm}/(m^{2}v_{\rm
dm}^{2}),\tilde{R}_{s}\equiv mv_{\rm dm}R_{s}\\}$. This can be easily seen by
rewriting $(i\partial_{t}+\nabla^{2}/2m-m\Phi)\psi=g|\psi|^{2}\psi$ in terms
of the rescaled variables $\tilde{\mathbf{x}}\equiv mv_{\rm dm}\mathbf{x}$,
$\tilde{t}\equiv mv_{\rm dm}^{2}t$ and $\tilde{\psi}\equiv\psi/\sqrt{\rho_{\rm
dm}/m}$ (note that space and time are normalized to the typical wavelength and
period of the DM waves, and $m\psi^{2}$ to their average density; the
relaxation time is simply $mv_{\rm dm}^{2}\tau_{\rm rel}=1/\tilde{g}^{2}$). In
this way, the EoM become

$\left(i\partial_{\tilde{t}}+\frac{\tilde{\nabla}^{2}}{2}-\tilde{\Phi}\right)\tilde{\psi}=\tilde{g}|\tilde{\psi}|^{2}\tilde{\psi}\,,$
(115)

where $\tilde{\Phi}=-\frac{\alpha/v_{\rm dm}}{\tilde{r}}$ for
$\tilde{r}>\tilde{R}_{s}$ (i.e. away from the Sun) and
$\tilde{\Phi}=-\frac{\alpha/v_{\rm
dm}}{2\tilde{R}_{s}}(3-\tilde{r}^{2}/\tilde{R}_{s}^{2})$ for
$\tilde{r}<\tilde{R}_{s}$ (inside the Sun). These redefinitions are useful
because the initial conditions only depend on $\alpha/v_{\rm dm}$, and read

$\tilde{\psi}_{w}=\int\frac{d^{3}\tilde{k}}{(2\pi)^{3}}\tilde{a}(\tilde{\mathbf{k}})e^{-i\frac{\tilde{k}^{2}\tilde{t}}{2}}\psi_{\tilde{\mathbf{k}}}(\tilde{\mathbf{x}})~{},\qquad\langle\tilde{a}^{*}(\tilde{\mathbf{k}})\tilde{a}(\tilde{\mathbf{k}}^{\prime})\rangle=\left[8\pi^{3/2}e^{-(\tilde{\mathbf{k}}-\hat{z})^{2}}\right](2\pi)^{3}\delta(\tilde{\mathbf{k}}-\tilde{\mathbf{k}}^{\prime})\,,$
(116)

where we took $\sigma=v_{\rm dm}/\sqrt{2}$ and $\mathbf{k}_{\rm dm}$ parallel
to $\hat{z}$, and it is understood that $R_{\star}\rightarrow\frac{v_{\rm
dm}}{\alpha}$ in $\psi_{\tilde{\mathbf{k}}}(\tilde{\mathbf{x}})$. In practice,
a realization of these initial conditions is (at $t=0$)

$\tilde{\psi}_{w,0}=A\int\frac{d^{3}\tilde{k}}{(2\pi)^{3}}e^{-\frac{1}{2}{(\tilde{\mathbf{k}}-\hat{z})^{2}}+i\gamma(\tilde{\mathbf{k}})}\psi_{\tilde{\mathbf{k}}}(\tilde{\mathbf{x}})\,,$
(117)

where $\gamma(\mathbf{k})$ is a random function of $\mathbf{k}$ taking values
in $[0,2\pi[$. The coefficient $A$ in Eq. (117) is fixed by the condition
$m\langle|\psi_{w,0}|^{2}\rangle=\rho_{\rm dm}$, i.e.
$\langle|\tilde{\psi}_{w,0}|^{2}\rangle=1$.

##### Evolution and systematics.

We solve Eq. (115) with initial conditions in Eq. (117) in a cubic box of
length $L$ with $N_{x}^{3}=64^{3}\div 128^{3}$ points.343434In practice, to
compute the initial conditions we used the plane wave form
$\psi_{\tilde{\mathbf{k}}}(\tilde{\mathbf{x}})=e^{i\tilde{\mathbf{k}}\cdot\tilde{\mathbf{x}}}$,
which is correct except close to $\tilde{r}=0$. Given that plane waves are not
hydrogen atom eigenstates, this leads to a small population of the bound
states in the initial conditions, which however does not affect the timescale
of the exponential growth, given the results of Section 4. The coefficient $A$
in Eq. (117) is in this case $A=8(\pi\tilde{L})^{3/2}$, where $\tilde{L}\equiv
mv_{\rm dm}L$ is the length in units of the inverse momentum of the waves. The
evolution is carried out with periodic boundary conditions and a standard
second-order pseudo-spectral method, employed e.g. in [31, 34, 57, 146, 58],
to which we refer for all the details. (Note that the two additional
parameters $\\{N_{x},\tilde{L}\\}$ need to be specified in a simulation.)

The EoM imply the conservation of the total mass, $m\int d^{3}x|\psi|^{2}$,
(which is automatic with the pseudo-spectral method we use) and of the
nonrelativistic energy $\mathcal{E}\equiv\int d^{3}x\epsilon$, where
$\epsilon$ is the energy density in Eq. (16). The time-independence of
$\mathcal{E}$ is therefore a nontrivial check of the numerical evolution. We
compute $\mathcal{E}$ throughout all the simulations, check its conservation
and rely on simulations where this is conserved better than
$\delta\mathcal{E}/\mathcal{E}<10^{-4}$ (this requires the time-step to be
small enough, see footnote 35).

As anticipated in Section 5, there are multiple competing requirements for a
simulation to capture the formation of the gravitational atom and be free from
systematic uncertainties.

  1. (a)

The lattice spacing $\Delta\equiv L/N_{x}$ should be small enough to resolve
the wavelength of the DM waves and the Sun radius. We checked that this
requires $\Delta mv_{\rm dm}=\tilde{L}/N_{x}\lesssim 1$ and
$\Delta/R_{s}=\tilde{L}/(N_{x}\tilde{R}_{s})\lesssim 1$ respectively, i.e.
there should be at least approximately one grid point per inverse DM momentum
and per Sun radius.

  2. (b)

The length $L$ should be large enough for the box to be in the infinite volume
limit, which requires $L/(2\pi/mv_{\rm dm})=\tilde{L}/2\pi\gg 1$. In other
words, the number of wavelengths per box side should be large;
$\tilde{L}/2\pi\gtrsim\mathcal{O}(10)$ is enough for this to happen. In
particular, this means that the mass of the gravitational atom is always a
small fraction of the total mass in the box.

  3. (c)

Finally, $R_{s}$ should be much smaller than $R_{\star}$ for the point-like
approximation of the Sun to hold, i.e. $R_{s}/R_{\star}=(\alpha/v_{\rm
dm})\tilde{R}_{s}\ll 1$.

In the simulations of the main text (Figure 11 in Section 5), we used
$|\tilde{g}|=0.006$ ($\tilde{g}>0$ for the two lowest values of
$2\pi\alpha/v_{\rm dm}$ and $\tilde{g}<0$ otherwise) and $\tilde{R}_{s}=0.8$
(except for the two highest values of $2\pi\alpha/v_{\rm dm}$, for which
$\tilde{R}_{s}=0.4$). We chose $\tilde{L}=40$ and performed simulations with
both $N_{x}^{3}=64^{3}$ and $N_{x}^{3}=128^{3}$ to check the lattice spacing
requirement (a). Finally, the function $\gamma(\tilde{\mathbf{k}})$ in the
initial conditions in Eq. (117) has been fixed to be the same (random)
function of $\tilde{\mathbf{k}}$ for all simulations in Figure 11.353535 We
use time-step $\Delta_{t}mv_{\rm dm}^{2}\simeq(\Delta mv_{\rm dm})^{2}/4$,
which is small enough for the energy $\mathcal{E}$ to be conserved as
described above. Note that a much larger time-step is sufficient for
$\mathcal{E}$ to be conserved at the $10^{-4}$ level before the gravitational
atom forms. However, after formation, energy conservation is worsened because
the shorter timescales associated to the ground state need to be resolved for
large $2\pi\alpha/v_{\rm dm}$.

Figure 20: The density profiles of Figure 11 for two subsequent times during
the exponential increase, corresponding to $t/\tau_{\rm
rel}\simeq\\{0.83,0.90\\}$ and $t/\tau_{\rm rel}\simeq\\{1.36,1.50\\}$ for
$\xi_{\rm foc}=8$ and $\xi_{\rm foc}=4$ respectively. Dashed lines represent
the ground-state profile.

Note that the conditions (a) and (b) are both satisfied with this choice of
parameters. As for (c), although $R_{s}<R_{\star}$ for all simulations, for
the largest values of $2\pi\alpha/v_{\rm dm}$ the ratio $R_{s}/R_{\star}$ is
not too much smaller than 1 (the worst case is $R_{s}/R_{\star}\simeq 0.8$,
for $2\pi\alpha/v_{\rm dm}=6.5$). However this does not seem to have a
significant effect on either the exponential-increase time and the
gravitational atom profile; see also Figure 11.

Figure 21:  _Left:_ The time evolution of the dark matter overdensity
$\rho(0)/\rho_{\rm dm}$ at the position of the Sun for $\xi_{\rm
foc}=2\pi\alpha/v_{\rm dm}=\\{6.5,2.5,1\\}$ and repulsive self-interactions
with $\tilde{g}=0.006\cdot\\{1,1,0.5\\}$. After an exponential increase, the
density tends to saturate at the predicted critical value $\rho_{\rm crit}$,
represented by an empty disk. _Right:_ The overdensity profiles during the
exponential increase, at times $t/\tau_{\rm rel}\simeq\\{1.5,1.7,1\\}$; dashed
lines represent the ground-state profile of Eq. (17).

##### Density profiles and repulsive self-interactions.

In Figure 20 we show the time evolution of the overdensity profile for the
simulations of Figure 11 in Section 5, for the two values $\xi_{\rm
foc}=2\pi\alpha/v_{\rm dm}=4,8$. This is obtained by taking the spherical
average of the field $|\tilde{\psi}|^{2}$ around the center of the Sun, $r=0$.
We show two times during the exponential increase; these times can be
extracted from the value of the density in Figure 11 (left). Note that the
profile oscillates in time (especially when the maximum density is not much
larger than the average density). Thus, as in Figure 11 (right), we plot the
profiles averaged over times longer than $1/mv^{2}_{\rm dm}$, so that the
oscillations can average out. From Figure 20, it is clear that the profile
resembles the ground state during the whole increase, in agreement with our
analytic expectations.

Finally, in Figure 21 we study more closely the system for repulsive self-
interactions, $g>0$. We show the overdensity $\rho(0)/\rho_{\rm dm}$ at the
position of the Sun as a function of time (left) and the density profiles
(right) in three simulations with $2\pi\alpha/v_{\rm dm}=\\{6.5,2.5,1\\}$ and
$\tilde{g}=0.006\cdot\\{1,1,0.5\\}$.363636As before, $\tilde{L}=40$ and we use
$\tilde{R}_{s}=\\{0.3,0.3,0.8\\}$ and $\tilde{N}_{x}=\\{128,128,64\\}$. First,
the exponential increase of the density is evident, given the logarithmic
scale in the $y$-axis. Soon after $t\simeq\tau_{\rm rel}$, the overdensity
tends to saturate to the expected critical value $\rho_{\rm crit}/\rho_{\rm
dm}=(2\pi\alpha/v_{\rm dm})^{2}/(2\pi^{2}\tilde{g})$ in Eq. (56), which is
shown (multiplied by $0.7$) by an empty disk. Additionally, the corresponding
density profiles before saturation are well reproduced by the ground state
(instead, the profile is deformed after saturation happens, and is flatter).
These simulations also show that the exponential increase in the density, i.e.
$\Gamma>0$, is likely to occur for values of $2\pi\alpha/v_{\rm dm}$ even as
small as 1, which is outside the regime of validity of our analytic
calculations.

As mentioned in the main text, we tested numerically values of $|\tilde{g}|$
and $\tilde{R}_{s}$ that are much larger than those relevant for our Solar
System. Thus, these cannot be directly used to justify the formation of the
solar halo. However, by performing simulations with various values of
$\tilde{g}$ and $\tilde{R}_{s}$, we checked that the time (relative to
$\tau_{\rm rel}$) when the density increases exponentially is independent of
these parameters, as long as $R_{s}\lesssim R_{\star}$, which confirms our
analytic analysis of Section 4.3. An example of this is the red curve of
Figure 21, for which $\tilde{g}$ is half in size than all the other
simulations. The results of this Appendix and Section 4 are therefore a
convincing evidence of the analytic picture developed in the main text, which
can be used in our Solar System for the relevant values of $|{g}|$ and
${R}_{s}$.

## References

  * [1] R. Peccei and H. R. Quinn, _CP Conservation in the Presence of Instantons_ , _Phys. Rev. Lett._ 38 (1977) 1440–1443.
  * [2] S. Weinberg, _A New Light Boson?_ , _Phys. Rev. Lett._ 40 (1978) 223–226.
  * [3] F. Wilczek, _Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons_, _Phys. Rev. Lett._ 40 (1978) 279–282.
  * [4] J. E. Kim, _Weak Interaction Singlet and Strong CP Invariance_ , _Phys. Rev. Lett._ 43 (1979) 103.
  * [5] M. A. Shifman, A. Vainshtein and V. I. Zakharov, _Can Confinement Ensure Natural CP Invariance of Strong Interactions?_ , _Nucl. Phys. B_ 166 (1980) 493–506.
  * [6] M. Dine, W. Fischler and M. Srednicki, _A Simple Solution to the Strong CP Problem with a Harmless Axion_ , _Phys. Lett. B_ 104 (1981) 199–202.
  * [7] A. Zhitnitsky, _On Possible Suppression of the Axion Hadron Interactions. (In Russian)_ , _Sov. J. Nucl. Phys._ 31 (1980) 260.
  * [8] P. W. Graham, D. E. Kaplan and S. Rajendran, _Cosmological Relaxation of the Electroweak Scale_ , _Phys. Rev. Lett._ 115 (2015) 221801, [1504.07551].
  * [9] K. Choi, H. Kim and T. Sekiguchi, _Dynamics of the cosmological relaxation after reheating_ , _Phys. Rev. D_ 95 (2017) 075008, [1611.08569].
  * [10] T. Flacke, C. Frugiuele, E. Fuchs, R. S. Gupta and G. Perez, _Phenomenology of relaxion-Higgs mixing_ , _JHEP_ 06 (2017) 050, [1610.02025].
  * [11] A. Banerjee, H. Kim and G. Perez, _Coherent relaxion dark matter_ , _Phys. Rev. D_ 100 (2019) 115026, [1810.01889].
  * [12] P. Svrcek and E. Witten, _Axions In String Theory_ , _JHEP_ 06 (2006) 051, [hep-th/0605206].
  * [13] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper and J. March-Russell, _String Axiverse_ , _Phys. Rev. D_ 81 (2010) 123530, [0905.4720].
  * [14] J. Preskill, M. B. Wise and F. Wilczek, _Cosmology of the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 127–132.
  * [15] L. F. Abbott and P. Sikivie, _A Cosmological Bound on the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 133–136.
  * [16] M. Dine and W. Fischler, _The Not So Harmless Axion_ , _Phys. Lett. B_ 120 (1983) 137–141.
  * [17] E. G. M. Ferreira, _Ultra-light dark matter_ , _Astron. Astrophys. Rev._ 29 (2021) 7, [2005.03254].
  * [18] M. Weber and W. de Boer, _Determination of the local dark matter density in our galaxy_ , _Astronomy and Astrophysics_ 509 (jan, 2010) A25.
  * [19] F. Nesti and P. Salucci, _The Local Dark Matter Density_ , _PoS_ DSU2012 (2012) 041, [1212.3670].
  * [20] J. Bovy and S. Tremaine, _On the local dark matter density_ , _Astrophys. J._ 756 (2012) 89, [1205.4033].
  * [21] J. I. Read, _The Local Dark Matter Density_ , _J. Phys. G_ 41 (2014) 063101, [1404.1938].
  * [22] N. W. Evans, C. A. J. O’Hare and C. McCabe, _Refinement of the standard halo model for dark matter searches in light of the Gaia Sausage_ , _Phys. Rev. D_ 99 (2019) 023012, [1810.11468].
  * [23] L. Necib, M. Lisanti, S. Garrison-Kimmel, A. Wetzel, R. Sanderson, P. F. Hopkins et al., _Under the FIRElight: Stellar Tracers of the Local Dark Matter Velocity Distribution in the Milky Way_ , _The Astrophysical Journal_ 883 (10, 2019) , [1810.12301].
  * [24] J. W. Foster, N. L. Rodd and B. R. Safdi, _Revealing the Dark Matter Halo with Axion Direct Detection_ , _Phys. Rev. D_ 97 (2018) 123006, [1711.10489].
  * [25] G. P. Centers et al., _Stochastic fluctuations of bosonic dark matter_ , _Nature Commun._ 12 (2021) 7321, [1905.13650].
  * [26] M. Lisanti, M. Moschella and W. Terrano, _Stochastic properties of ultralight scalar field gradients_ , _Phys. Rev. D_ 104 (2021) 055037, [2107.10260].
  * [27] D. J. Kaup, _Klein-Gordon Geon_ , _Phys. Rev._ 172 (1968) 1331–1342.
  * [28] R. Ruffini and S. Bonazzola, _Systems of selfgravitating particles in general relativity and the concept of an equation of state_ , _Phys. Rev._ 187 (1969) 1767–1783.
  * [29] E. W. Kolb and I. I. Tkachev, _Axion miniclusters and Bose stars_ , _Phys. Rev. Lett._ 71 (1993) 3051–3054, [hep-ph/9303313].
  * [30] H.-Y. Schive, T. Chiueh and T. Broadhurst, _Cosmic Structure as the Quantum Interference of a Coherent Dark Wave_ , _Nature Phys._ 10 (2014) 496–499, [1406.6586].
  * [31] D. G. Levkov, A. G. Panin and I. I. Tkachev, _Gravitational Bose-Einstein condensation in the kinetic regime_ , _Phys. Rev. Lett._ 121 (2018) 151301, [1804.05857].
  * [32] B. Eggemeier and J. C. Niemeyer, _Formation and mass growth of axion stars in axion miniclusters_ , _Phys. Rev. D_ 100 (2019) 063528, [1906.01348].
  * [33] J. Veltmaat, B. Schwabe and J. C. Niemeyer, _Baryon-driven growth of solitonic cores in fuzzy dark matter halos_ , _Phys. Rev. D_ 101 (2020) 083518, [1911.09614].
  * [34] 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 (2021) 083022, [2011.01333].
  * [35] M. Colpi, S. L. Shapiro and I. Wasserman, _Boson Stars: Gravitational Equilibria of Selfinteracting Scalar Fields_ , _Phys. Rev. Lett._ 57 (1986) 2485–2488.
  * [36] P.-H. Chavanis, _Mass-radius relation of Newtonian self-gravitating Bose-Einstein condensates with short-range interactions: I. Analytical results_ , _Phys. Rev. D_ 84 (2011) 043531, [1103.2050].
  * [37] P. 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 (2011) 043532, [1103.2054].
  * [38] J. Eby, P. Suranyi, C. Vaz and L. Wijewardhana, _Axion Stars in the Infrared Limit_ , _JHEP_ 03 (2015) 080, [1412.3430].
  * [39] 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 (2016) 291–295, [1508.05395].
  * [40] M. Minamitsuji, _Vector boson star solutions with a quartic order self-interaction_ , _Phys. Rev. D_ 97 (2018) 104023, [1805.09867].
  * [41] M. A. Amin and P. Mocz, _Formation, gravitational clustering, and interactions of nonrelativistic solitons in an expanding universe_ , _Phys. Rev. D_ 100 (2019) 063507, [1902.07261].
  * [42] M. Jain and M. A. Amin, _Polarized solitons in higher-spin wave dark matter_ , _Phys. Rev. D_ 105 (2022) 056019, [2109.04892].
  * [43] M. Gorghetto, E. Hardy, J. March-Russell, N. Song and S. M. West, _Dark photon stars: formation and role as dark matter substructure_ , _JCAP_ 08 (2022) 018, [2203.10100].
  * [44] A. Banerjee, D. Budker, J. Eby, H. Kim and G. Perez, _Relaxion Stars and their detection via Atomic Physics_ , _Commun. Phys._ 3 (2020) 1, [1902.08212].
  * [45] A. Banerjee, D. Budker, J. Eby, V. V. Flambaum, H. Kim, O. Matsedonskyi et al., _Searching for Earth/Solar Axion Halos_ , _JHEP_ 09 (2020) 004, [1912.04295].
  * [46] A. Kryemadhi, M. Maroudas, A. Mastronikolis and K. Zioutas, _Gravitational focusing effects on streaming dark matter as a new detection concept_ , 2210.07367.
  * [47] D. H. H. Hoffmann, J. Jacoby and K. Zioutas, _Gravitational lensing by the Sun of non-relativistic penetrating particles_ , _Astropart. Phys._ 20 (2003) 73–78.
  * [48] M. S. Alenazi and P. Gondolo, _Phase-space distribution of unbound dark matter near the Sun_ , _Phys. Rev. D_ 74 (2006) 083518, [astro-ph/0608390].
  * [49] B. R. Patla, R. J. Nemiroff, D. H. H. Hoffmann and K. Zioutas, _Flux Enhancement of Slow-moving Particles by Sun or Jupiter: Can they be Detected on Earth?_ , _Astrophys. J._ 780 (2014) 158, [1305.2454].
  * [50] G. Prézeau, _Dense Dark Matter Hairs Spreading Out from Earth, Jupiter and Other Compact Bodies_ , _Astrophys. J._ 814 (2015) 122, [1507.07009].
  * [51] Y. Sofue, _Gravitational Focusing of Low-Velocity Dark Matter on the Earth’s Surface_ , _Galaxies_ 8 (2020) 42, [2005.08252].
  * [52] S. K. Lee, M. Lisanti, A. H. G. Peter and B. R. Safdi, _Effect of Gravitational Focusing on Annual Modulation in Dark-Matter Direct-Detection Experiments_ , _Phys. Rev. Lett._ 112 (2014) 011301, [1308.1953].
  * [53] H. Kim and A. Lenoci, _Gravitational focusing of wave dark matter_ , _Phys. Rev. D_ 105 (2022) 063032, [2112.05718].
  * [54] R. J. Crewther, P. Di Vecchia, G. Veneziano and E. Witten, _Chiral Estimate of the Electric Dipole Moment of the Neutron in Quantum Chromodynamics_ , _Phys. Lett. B_ 88 (1979) 123.
  * [55] G. Grilli di Cortona, E. Hardy, J. Pardo Vega and G. Villadoro, _The QCD axion, precisely_ , _JHEP_ 01 (2016) 034, [1511.02867].
  * [56] A. Sommerfeld, _Über die beugung und bremsung der elektronen_ , _Annalen der Physik_ 403 (1931) 257–330, [https://onlinelibrary.wiley.com/doi/pdf/10.1002/andp.19314030302].
  * [57] 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 (2022) 023009, [2109.11474].
  * [58] 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 (2022) 043512, [2110.08921].
  * [59] J. Eby, M. Leembruggen, P. Suranyi and L. C. R. Wijewardhana, _Collapse of Axion Stars_ , _JHEP_ 12 (2016) 066, [1608.06911].
  * [60] D. G. Levkov, A. G. Panin and I. I. Tkachev, _Relativistic axions from collapsing Bose stars_ , _Phys. Rev. Lett._ 118 (2017) 011301, [1609.03611].
  * [61] N. P. Pitjev and E. V. Pitjeva, _Constraints on dark matter in the solar system_ , _Astron. Lett._ 39 (2013) 141–149, [1306.5534].
  * [62] Y.-D. Tsai, J. Eby, J. Arakawa, D. Farnocchia and M. S. Safronova, _New Constraints on Dark Matter and Cosmic Neutrino Profiles through Gravity_ , 2210.03749.
  * [63] J. I. Read, G. Lake, O. Agertz and V. P. Debattista, _Thin, thick and dark discs in LCDM_ , _Mon. Not. Roy. Astron. Soc._ 389 (2008) 1041–1057, [0803.2714].
  * [64] C. W. Purcell, J. S. Bullock and M. Kaplinghat, _The Dark Disk of the Milky Way_ , _apj_ 703 (Oct., 2009) 2275–2284, [0906.5348].
  * [65] D. Budker, P. W. Graham, M. Ledbetter, S. Rajendran and A. Sushkov, _Proposal for a Cosmic Axion Spin Precession Experiment (CASPEr)_ , _Phys. Rev. X_ 4 (2014) 021030, [1306.6089].
  * [66] J. A. Dror, S. Gori, J. M. Leedom and N. L. Rodd, _Sensitivity of Spin-Precession Axion Experiments_ , _Phys. Rev. Lett._ 130 (2023) 181801, [2210.06481].
  * [67] N. Bar, K. Blum, T. Lacroix and P. Panci, _Looking for ultralight dark matter near supermassive black holes_ , _JCAP_ 07 (2019) 045, [1905.11745].
  * [68] E. Y. Davies and P. Mocz, _Fuzzy Dark Matter Soliton Cores around Supermassive Black Holes_ , _Mon. Not. Roy. Astron. Soc._ 492 (2020) 5721–5729, [1908.04790].
  * [69] J. H.-H. Chan, S. Sibiryakov and W. Xue, _Condensation and Evaporation of Boson Stars_ , 2207.04057.
  * [70] A. Arvanitaki, J. Huang and K. Van Tilburg, _Searching for dilaton dark matter with atomic clocks_ , _Phys. Rev. D_ 91 (2015) 015015, [1405.2925].
  * [71] J. Fan, _Ultralight Repulsive Dark Matter and BEC_ , _Phys. Dark Univ._ 14 (2016) 84–94, [1603.06580].
  * [72] J. A. R. Cembranos, A. L. Maroto, S. J. Núñez Jareño and H. Villarrubia-Rojo, _Constraints on anharmonic corrections of Fuzzy Dark Matter_ , _JHEP_ 08 (2018) 073, [1805.08112].
  * [73] J. Jaeckel, V. M. Mehta and L. T. Witkowski, _Monodromy Dark Matter_ , _JCAP_ 01 (2017) 036, [1605.01367].
  * [74] R. T. Co, L. J. Hall and K. Harigaya, _QCD Axion Dark Matter with a Small Decay Constant_ , _Phys. Rev. Lett._ 120 (2018) 211602, [1711.10486].
  * [75] R. T. Co, E. Gonzalez and K. Harigaya, _Axion Misalignment Driven to the Hilltop_ , _JHEP_ 05 (2019) 163, [1812.11192].
  * [76] D. Cyncynates, T. Giurgica-Tiron, O. Simon and J. O. Thompson, _Resonant nonlinear pairs in the axiverse and their late-time direct and astrophysical signatures_ , _Phys. Rev. D_ 105 (2022) 055005, [2109.09755].
  * [77] CAST collaboration, V. Anastassopoulos et al., _New CAST Limit on the Axion-Photon Interaction_ , _Nature Phys._ 13 (2017) 584–590, [1705.02290].
  * [78] Fermi-LAT collaboration, M. Ajello et al., _Search for Spectral Irregularities due to Photon–Axionlike-Particle Oscillations with the Fermi Large Area Telescope_ , _Phys. Rev. Lett._ 116 (2016) 161101, [1603.06978].
  * [79] C. S. Reynolds, M. C. D. Marsh, H. R. Russell, A. C. Fabian, R. Smith, F. Tombesi et al., _Astrophysical limits on very light axion-like particles from Chandra grating spectroscopy of NGC 1275_ , _Astrophys. J._ 890 (2020) 59, [1907.05475].
  * [80] M. J. Dolan, F. J. Hiskens and R. R. Volkas, _Advancing globular cluster constraints on the axion-photon coupling_ , _JCAP_ 10 (2022) 096, [2207.03102].
  * [81] ADMX collaboration, N. Du et al., _A Search for Invisible Axion Dark Matter with the Axion Dark Matter Experiment_ , _Phys. Rev. Lett._ 120 (2018) 151301, [1804.05750].
  * [82] ADMX collaboration, S. J. Asztalos et al., _An Improved RF cavity search for halo axions_ , _Phys. Rev._ D69 (2004) 011101, [astro-ph/0310042].
  * [83] ADMX collaboration, S. J. Asztalos et al., _A SQUID-based microwave cavity search for dark-matter axions_ , _Phys. Rev. Lett._ 104 (2010) 041301, [0910.5914].
  * [84] P. W. Graham and S. Rajendran, _New Observables for Direct Detection of Axion Dark Matter_ , _Phys. Rev. D_ 88 (2013) 035023, [1306.6088].
  * [85] G. Ruoso, A. Lombardi, A. Ortolan, R. Pengo, C. Braggio, G. Carugno et al., _The QUAX proposal: a search of galactic axion with magnetic materials_ , _J. Phys. Conf. Ser._ 718 (2016) 042051, [1511.09461].
  * [86] H. Liu, B. D. Elwood, M. Evans and J. Thaler, _Searching for Axion Dark Matter with Birefringent Cavities_ , _Phys. Rev. D_ 100 (2019) 023548, [1809.01656].
  * [87] K. Nagano, T. Fujita, Y. Michimura and I. Obata, _Axion Dark Matter Search with Interferometric Gravitational Wave Detectors_ , _Phys. Rev. Lett._ 123 (2019) 111301, [1903.02017].
  * [88] I. Obata, T. Fujita and Y. Michimura, _Optical Ring Cavity Search for Axion Dark Matter_ , _Phys. Rev. Lett._ 121 (2018) 161301, [1805.11753].
  * [89] W. DeRocco and A. Hook, _Axion interferometry_ , _Phys. Rev._ D98 (2018) 035021, [1802.07273].
  * [90] P. Sikivie, N. Sullivan and D. B. Tanner, _Proposal for Axion Dark Matter Detection Using an LC Circuit_ , _Phys. Rev. Lett._ 112 (2014) 131301, [1310.8545].
  * [91] S. Chaudhuri, P. W. Graham, K. Irwin, J. Mardon, S. Rajendran and Y. Zhao, _Radio for hidden-photon dark matter detection_ , _Phys. Rev._ D92 (2015) 075012, [1411.7382].
  * [92] Y. Kahn, B. R. Safdi and J. Thaler, _Broadband and Resonant Approaches to Axion Dark Matter Detection_ , _Phys. Rev. Lett._ 117 (2016) 141801, [1602.01086].
  * [93] M. Silva-Feaver et al., _Design Overview of DM Radio Pathfinder Experiment_ , _IEEE Trans. Appl. Supercond._ 27 (2017) 1400204, [1610.09344].
  * [94] J. L. Ouellet et al., _First Results from ABRACADABRA-10 cm: A Search for Sub- $\mu$eV Axion Dark Matter_, _Phys. Rev. Lett._ 122 (2019) 121802, [1810.12257].
  * [95] DMRadio collaboration, A. AlShirawi et al., _Electromagnetic modeling and science reach of DMRadio-m 3_, 2302.14084.
  * [96] P. Sikivie, _Axion Dark Matter Detection using Atomic Transitions_ , _Phys. Rev. Lett._ 113 (2014) 201301, [1409.2806].
  * [97] H. Kim and G. Perez, _Oscillations of atomic energy levels induced by QCD axion dark matter_ , 2205.12988.
  * [98] C. Abel et al., _Search for Axionlike Dark Matter through Nuclear Spin Precession in Electric and Magnetic Fields_ , _Phys. Rev._ X7 (2017) 041034, [1708.06367].
  * [99] J. E. Kim and G. Carosi, _Axions and the Strong CP Problem_ , _Rev. Mod. Phys._ 82 (2010) 557–602, [0807.3125].
  * [100] P. W. Graham, I. G. Irastorza, S. K. Lamoreaux, A. Lindner and K. A. van Bibber, _Experimental Searches for the Axion and Axion-Like Particles_ , _Ann. Rev. Nucl. Part. Sci._ 65 (2015) 485–514, [1602.00039].
  * [101] D. J. Gross, R. D. Pisarski and L. G. Yaffe, _QCD and Instantons at Finite Temperature_ , _Rev. Mod. Phys._ 53 (1981) 43.
  * [102] J. Breit, S. Gupta and A. Zaks, _Cold bose stars_ , _Physics Letters B_ 140 (1984) 329 – 332.
  * [103] L. Hui, J. P. Ostriker, S. Tremaine and E. Witten, _Ultralight scalars as cosmological dark matter_ , _Phys. Rev. D_ 95 (2017) 043541, [1610.08297].
  * [104] G. C. Wick, _The evaluation of the collision matrix_ , _Phys. Rev._ 80 (Oct, 1950) 268–272.
  * [105] R. Schoenrich, J. Binney and W. Dehnen, _Local Kinematics and the Local Standard of Rest_ , _Mon. Not. Roy. Astron. Soc._ 403 (2010) 1829, [0912.3693].
  * [106] A.-C. Eilers, D. W. Hogg, H.-W. Rix and M. K. Ness, _The circular velocity curve of the milky way from 5 to 25 kpc_ , _The Astrophysical Journal_ 871 (jan, 2019) 120.
  * [107] T. Helfer, D. J. E. Marsh, K. Clough, M. Fairbairn, E. A. Lim and R. Becerril, _Black hole formation from axion stars_ , _JCAP_ 03 (2017) 055, [1609.04724].
  * [108] J. Eby, M. Leembruggen, P. Suranyi and L. C. R. Wijewardhana, _QCD Axion Star Collapse with the Chiral Potential_ , _JHEP_ 06 (2017) 014, [1702.05504].
  * [109] F. Michel and I. G. Moss, _Relativistic collapse of axion stars_ , _Phys. Lett. B_ 785 (2018) 9–13, [1802.10085].
  * [110] J. Eby, P. Suranyi and L. C. R. Wijewardhana, _The Lifetime of Axion Stars_ , _Mod. Phys. Lett. A_ 31 (2016) 1650090, [1512.01709].
  * [111] J. Eby, M. Ma, P. Suranyi and L. C. R. Wijewardhana, _Decay of Ultralight Axion Condensates_ , _JHEP_ 01 (2018) 066, [1705.05385].
  * [112] J. Eby, S. Shirai, Y. V. Stadnik and V. Takhistov, _Probing relativistic axions from transient astrophysical sources_ , _Phys. Lett. B_ 825 (2022) 136858, [2106.14893].
  * [113] P. J. Fox, N. Weiner and H. Xiao, _Recurrent Axinovae and their Cosmological Constraints_ , 2302.00685.
  * [114] X. Du, D. J. E. Marsh, M. Escudero, A. Benson, D. Blas, C. K. Pooni et al., _Soliton Merger Rates and Enhanced Axion Dark Matter Decay_ , 2301.09769.
  * [115] M. Escudero, C. K. Pooni, M. Fairbairn, D. Blas, X. Du and D. J. E. Marsh, _Axion Star Explosions: A New Source for Axion Indirect Detection_ , 2302.10206.
  * [116] C. G. Boehmer and T. Harko, _Can dark matter be a Bose-Einstein condensate?_ , _JCAP_ 06 (2007) 025, [0705.4158].
  * [117] Y.-D. Tsai, J. Eby and M. S. Safronova, _Direct detection of ultralight dark matter bound to the Sun with space quantum sensors_ , _Nature Astron._ 7 (2023) 113–121, [2112.07674].
  * [118] A. Caputo, M. Regis, M. Taoso and S. J. Witte, _Detecting the Stimulated Decay of Axions at RadioFrequencies_ , _JCAP_ 03 (2019) 027, [1811.08436].
  * [119] M. P. Hertzberg and E. D. Schiappacasse, _Dark Matter Axion Clump Resonance of Photons_ , _JCAP_ 11 (2018) 004, [1805.00430].
  * [120] M. P. Hertzberg, Y. Li and E. D. Schiappacasse, _Merger of Dark Matter Axion Clumps and Resonant Photon Emission_ , _JCAP_ 07 (2020) 067, [2005.02405].
  * [121] D. G. Levkov, A. G. Panin and I. I. Tkachev, _Radio-emission of axion stars_ , _Phys. Rev. D_ 102 (2020) 023501, [2004.05179].
  * [122] D. Antypas et al., _New Horizons: Scalar and Vector Ultralight Dark Matter_ , 2203.14915.
  * [123] Y. Hochberg, T. Lin and K. M. Zurek, _Absorption of light dark matter in semiconductors_ , _Phys. Rev. D_ 95 (2017) 023013, [1608.01994].
  * [124] K. Van Tilburg, _Stellar basins of gravitationally bound particles_ , _Phys. Rev. D_ 104 (2021) 023019, [2006.12431].
  * [125] R. Lasenby and K. Van Tilburg, _Dark photons in the solar basin_ , _Phys. Rev. D_ 104 (2021) 023020, [2008.08594].
  * [126] P.-H. Chavanis, _Growth of perturbations in an expanding universe with Bose-Einstein condensate dark matter_ , _Astron. Astrophys._ 537 (2012) A127, [1103.2698].
  * [127] A. Suárez and P.-H. Chavanis, _Hydrodynamic representation of the Klein-Gordon-Einstein equations in the weak field limit: General formalism and perturbations analysis_ , _Phys. Rev. D_ 92 (2015) 023510, [1503.07437].
  * [128] J.-c. Hwang and H. Noh, _Axion as a Cold Dark Matter candidate_ , _Phys. Lett. B_ 680 (2009) 1–3, [0902.4738].
  * [129] P.-H. Chavanis, _Collapse of a self-gravitating Bose-Einstein condensate with attractive self-interaction_ , _Phys. Rev. D_ 94 (2016) 083007, [1604.05904].
  * [130] W. Hu, R. Barkana and A. Gruzinov, _Cold and fuzzy dark matter_ , _Phys. Rev. Lett._ 85 (2000) 1158–1161, [astro-ph/0003365].
  * [131] A. Arvanitaki, S. Dimopoulos, M. Galanis, L. Lehner, J. O. Thompson and K. Van Tilburg, _Large-misalignment mechanism for the formation of compact axion structures: Signatures from the QCD axion to fuzzy dark matter_ , _Phys. Rev. D_ 101 (2020) 083014, [1909.11665].
  * [132] SDSS collaboration, M. Tegmark et al., _Cosmological parameters from SDSS and WMAP_ , _Phys. Rev. D_ 69 (2004) 103501, [astro-ph/0310723].
  * [133] M. M. Ivanov, M. Simonović and M. Zaldarriaga, _Cosmological Parameters from the BOSS Galaxy Power Spectrum_ , _JCAP_ 05 (2020) 042, [1909.05277].
  * [134] S. A. Appleby, C. Park, S. E. Hong, H. S. Hwang and J. Kim, _Cosmological Parameter Estimation from the Two-Dimensional Genus Topology – Measuring the Shape of the Matter Power Spectrum_ , _Astrophys. J._ 896 (2020) 145, [2004.01424].
  * [135] B. Li, T. Rindler-Daller and P. R. Shapiro, _Cosmological Constraints on Bose-Einstein-Condensed Scalar Field Dark Matter_ , _Phys. Rev. D_ 89 (2014) 083536, [1310.6061].
  * [136] A. Arvanitaki and S. Dubovsky, _Exploring the String Axiverse with Precision Black Hole Physics_ , _Phys. Rev. D_ 83 (2011) 044026, [1004.3558].
  * [137] A. Arvanitaki, M. Baryakhtar, S. Dimopoulos, S. Dubovsky and R. Lasenby, _Black Hole Mergers and the QCD Axion at Advanced LIGO_ , _Phys. Rev. D_ 95 (2017) 043001, [1604.03958].
  * [138] M. Baryakhtar, M. Galanis, R. Lasenby and O. Simon, _Black hole superradiance of self-interacting scalar fields_ , _Phys. Rev. D_ 103 (2021) 095019, [2011.11646].
  * [139] N. P. Branco, R. Z. Ferreira and J. a. G. Rosa, _Superradiant axion clouds around asteroid-mass primordial black holes_ , _JCAP_ 04 (2023) 003, [2301.01780].
  * [140] A. Arvanitaki, M. Baryakhtar and X. Huang, _Discovering the QCD Axion with Black Holes and Gravitational Waves_ , _Phys. Rev. D_ 91 (2015) 084011, [1411.2263].
  * [141] C. Ünal, F. Pacucci and A. Loeb, _Properties of ultralight bosons from heavy quasar spins via superradiance_ , _JCAP_ 05 (2021) 007, [2012.12790].
  * [142] F. Chadha-Day, B. Garbrecht and J. McDonald, _Superradiance in stars: non-equilibrium approach to damping of fields in stellar media_ , _JCAP_ 12 (2022) 008, [2207.07662].
  * [143] W. DeRocco, S. Wegsman, B. Grefenstette, J. Huang and K. Van Tilburg, _First Indirect Detection Constraints on Axions in the Solar Basin_ , _Phys. Rev. Lett._ 129 (2022) 101101, [2205.05700].
  * [144] A. S. Dmitriev, D. G. Levkov, A. G. Panin and I. I. Tkachev, _Self-similar growth of Bose stars_ , 2305.01005.
  * [145] J. H.-H. Chan, S. Sibiryakov and W. Xue, _Boson Star Normal Modes_ , 2304.13054.
  * [146] F. Edwards, E. Kendall, S. Hotchkiss and R. Easther, _PyUltraLight: A Pseudo-Spectral Solver for Ultralight Dark Matter Dynamics_ , _JCAP_ 10 (2018) 027, [1807.04037].