A generic formation mechanism of ultralight dark matter solar halos

As-yet undiscovered light bosons may constitute all or part of the dark matter (DM) of our Universe, and are expected to have (weak) self-interactions. We show that the quartic self-interactions generically induce the capture of dark matter from the surrounding halo by external gravitational potentials such as those of stars, including the Sun. This leads to the subsequent formation of dark matter bound states supported by such external potentials, resembling gravitational atoms (e.g. a solar halo around our own Sun). Their growth is governed by the ratio ξ foc ≡ λdB/R ⋆ between the de Broglie wavelength of the incoming DM waves, λdB, and the radius of the ground state R ⋆. For ξ foc ≲ 1, the gravitational atom grows to an (underdense) steady state that balances the capture of particles and the inverse (stripping) process. For ξ foc ≳ 1, a significant gravitational-focusing effect leads to exponential accumulation of mass from the galactic DM halo into the gravitational atom. For instance, a dark matter axion with mass of the order of 10-14 eV and decay constant between 107 and 108 GeV would form a dense halo around the Sun on a timescale comparable to the lifetime of the Solar System, leading to a local DM density at the position of the Earth 𝒪(104) times larger than that expected in the standard halo model. For attractive self-interactions, after its formation, the gravitational atom is destabilized at a large density, which leads to its collapse; this is likely to be accompanied by emission of relativistic bosons (a `Bosenova').

We denote these fields as ultralight if their mass m is smaller than about 1 eV. This implies that their occupation number in galaxies such as the Milky Way is macroscopic. As a result, rather than behaving as a collection of individual particles, their evolution is better approximated by their classical equations of motion (EoM), the Schrödinger-Poisson equations, whose free solutions are scalar waves. See [17] for a recent review of ultralight DM.
The typical parameters assumed in experimental searches for DM, derived from the large-scale properties of the Milky Way halo, are an energy density ρ dm ≃ 0.4 GeV/cm 3 and velocity v dm ≃ 240 km/s (see, for example, [18][19][20][21][22][23]). For ultralight DM (ULDM), the coherence time and other properties related to the ULDM stochasticity are also important for the theoretical interpretation of these searches [24][25][26].
However, overdensities at scales much smaller than the galaxy can modify this picture. In the literature, two possible scenarios have been considered thus far. First, a large fraction of the ULDM could reside in overdense objects supported by their own gravitational interactions [27][28][29][30][31][32][33][34] (see [35][36][37][38] for discussion of non-gravitational self-interactions). 1 However, unless the DM-particle mass is close to the eV scale, these objects have a small encounter rate with the Earth over the typical experimental timescale of a year, and/or a negligible overdensity, and thus unlikely to be detected via terrestrial experiments [44].
In the second scenario, the ULDM is trapped in the gravitational field of a massive astrophysical object, such as a star or a planet, resulting in a halo surrounding the object. The density of such a 'solar halo' or 'Earth halo', which resembles a 'gravitational atom' with the Sun or Earth as its center, 2 could be many orders of magnitudes larger than what is conventionally inferred from the properties of the galactic halo. This can potentially lead to a dramatic change in the way we search for DM [44,45].
However, the question of how such a halo could have formed remained unanswered. The obvious challenge stems from the fact that the typical relative velocity between the DM in our galaxy and the Solar System is of the order of a few hundred km/s which, for instance, is somewhat larger than the escape velocity from the Sun at e.g. an AU distance. Therefore, an efficient mechanism of energy loss is required for the DM particles to be captured by the solar gravitational potential.
In this work we show that there is in fact a class of scalar/axion theories where, in the background of the gravitational potential of astrophysical sources, DM capture becomes efficient. Our basic observation is that there is a regime of masses and quartic self-interactions where processes mediated by such interactions induce the capture of one of two ULDM waves that scatter with each other. The quantity that determines the efficiency of the above capture processes is the ratio between the de Broglie wavelength of an incoming ULDM wave, λ dB , and the gravitational Bohr radius, R ⋆ , characterizing the size of the halo: When ξ foc is larger than unity, the density of the DM waves is increased around and behind the astrophysical object -an effect known as gravitational focusing, leading to an enhanced capture rate (see the next Section for more details, and [46][47][48][49][50][51][52][53] for other consequences of gravitational focusing). In the case Figure 1: The formation of a gravitational atom bound to the Sun, i.e. a solar halo. The dark matter in our galactic halo, traveling with average velocity v dm ≃ 240 km/ s with respect to the Sun (and dispersion σ ≃ v dm / √ 2), experiences 2 → 2 scattering processes mediated by the quartic self-coupling λ that decrease (or increase) the velocity of the particles. When this occurs in the gravitational potential of the Sun, the dark matter can become trapped into a bound configuration, which resembles the ground state (nlm = 100) of the hydrogen atom. Excited states (n > 1) are also populated but are likely to decay into the ground state, see Appendix C.
of attractive self-interactions, as for axions (see e.g. [54,55]), we identify three phases in the evolution of such a halo, with typical timescale of the order of the so-called relaxation time: (i) slow growth from DM capture, (ii) exponential capture phase, and, after the density has grown beyond some critical value (defined below), (iii) collapse leading to a rapid emission of relativistic scalars, known as a Bosenova.
The paper is organized as follows. In the next Section we discuss the basic dynamics giving rise to the ULDM halo and summarize the main results of the paper. After defining the bound states in Section 3, we discuss their formation via DM capture in Section 4, and check and refine these results with numerical simulations in Section 5. Some phenomenological implications of a gravitational atom in the Solar System are discussed in Section 6, and constraints on the relevant values of mass and self-interaction couplings (mostly from structure formation) are discussed in Section 7. We conclude in Section 8.

Basic mechanism
To understand the formation of the ULDM halo, we consider an astrophysical body of mass M , say our Sun, moving in the background of virialized DM. In the rest frame of the body and far from it, the DM has a velocity distribution with mean v dm and variance of order σ 2 ≃ v 2 dm /2. The equation of motion (EoM) of the ULDM field is the Schrödinger equation with an external gravitational potential ∝ 1/r. As a result, the stationary solutions include bound states corresponding to a gravitational atom, with ground state radius R ⋆ ≡ (mα) −1 , where α ≡ GM m is the gravitational coupling of the DM to the body, and is the analogue of the fine-structure constant. Explicitly, the radius is where M ⊙ is the solar mass. Therefore, R ⋆ is fully determined by m and M . These initial conditions can be thought of as all of the ULDM occupying the continuum (unbound states of the atom) at early times. Despite being initially empty, bound states do become populated from Figure 2: A schematic representation of the ULDM waves (black lines) with de Broglie wavelength λ dB , and of the ground state of the atom (spheres), whose radius is R ⋆ ≡ (mα) −1 , where α ≡ GM m. We show the two limits ξ foc ≡ λ dB /R ⋆ ≪ 1 (left) and ξ foc ≫ 1 (right), which correspond to large/small DM velocity respectively. In the latter, the (otherwise plane) waves are distorted and amplified close to and behind the astrophysical body, an effect known as gravitational focusing. Correspondingly, the dark matter density (dashed lines) is enhanced. In the regime ξ foc ≫ 1, dark matter capture overcomes stripping. the DM in the continuum via processes mediated by the quartic self-interactions, which we parameterize as where ϕ is the ULDM field and λ the quartic coupling. Indeed, by treating Eq. (3) as a perturbation in the EoM, in Section 4 we show that the total mass M ⋆ of ULDM bound to the Sun changes in time aṡ where C, Γ 1 , and Γ 2 are positive constants (that can be computed semi-analytically) resulting from the self-interactions, and at leading order proportional to λ 2 . The first term in right hand side of Eq. (4), here referred to as 'capture', contributes positively to the bound mass and can be interpreted as arising from the 2 → 2 process in Figure 1, where two unbound particles scatter into a bound particle and a (more energetic) unbound one, which escapes to infinity. The second, 'stimulated capture', is a consequence of the Bose enhancement of the indistinguishable bosons and represents the same process, but is proportional to M ⋆ itself and so is effective only when the bound state is already populated. The last term is negative and simply represents the depletion of the bound states via the inverse process, 'stripping'. Consequently, the change in M ⋆ can be understood as the net effect of capture and stripping. Substantial DM capture occurs if the rate of stimulated capture exceeds that of stripping, i.e. if Γ ≡ Γ 1 − Γ 2 > 0. Our crucial observation is that this occurs when the dark matter field is significantly gravitationally focused by the external body, namely if is larger than 1 (here we expressed the typical de Broglie wavelength λ dB = 2π/mv dm in terms of the mean DM velocity v dm ). Indeed, if λ dB is smaller than R ⋆ (ξ foc ≪ 1, see Figure 2, left), corresponding to v dm > 2πα, the incoming waves are so fast that the effect of the gravitational potential on them is negligible, and they behave as plane waves everywhere. The kinetic energy mv 2 dm /2 of the incoming particles is much larger than the binding energy of the ground state, −mα 2 /2, and capture is inefficient.
Instead, if λ dB is larger than R ⋆ (ξ foc ≫ 1, see Figure 2, right), the dynamics of the waves close to the Sun is dominated by the Sun itself: Their magnitude increases in the region close and behind it, i.e. they are gravitationally 'focused'. This effect is analogous to the Sommerfeld enhancement, which When λ dB > R ⋆ , i.e. ξ foc > 1, which happens for m ≳ 1.7 · 10 −14 eV, gravitational focusing is relevant (light gray region); the corresponding value of the radius R ⋆ when the two are equal is close to an AU. At large masses, m ≳ 2 · 10 −13 eV (darker region), the radius of the Sun R ⊙ (green line) is larger than the would-be ground state radius R ⋆ : our treatment of the bound states does not apply, and the effect of the gravitational focusing changes.
is the increase of the scattering cross section (of a factor S) of two quantum mechanical particles in the presence of a central potential [56]. Indeed, the DM density at the center increases as ξ foc /(1 − e −ξ foc ), see also Figure 7, which is the same as the Sommerfeld enhancement factor S. The enhancement is ultimately related to the fact that the (initially slow) particles increase their velocity as a result of the gravitational potential, gaining a kinetic energy of the order of their binding energy. In this regime, their kinetic energy is small enough that further (order-one) energy changes, resulting from particles scattering through e.g. their self-interactions, have a chance of trapping them in the gravitational potential well. Additionally, incoming particles are not energetic enough (on average) to ionize the ground state, without getting themselves captured. Thus, we expect that stripping is suppressed compared to stimulated capture in this case. In Section 4.3 we prove that these intuitions are correct by showing that Γ > 0 if ξ foc ≫ 1 and Γ < 0 if ξ foc ≪ 1. By simulating the system numerically, in Section 5 we provide evidence that the transition is likely to happen around ξ foc ≃ 1.
In Figure 3 we compare λ dB and R ⋆ in the Solar System as a function of the mass m. From Eq. (5) and Figure 3, the condition ξ foc ≳ 1 is verified for α ≳ 10 −4 or, equivalently, m larger than approximately 10 −14 eV, for which the corresponding ground state turns out to have a radius of order AU or smaller; see Eq. (2). This coincidence can be easily understood by noticing that ξ foc = 1 corresponds to a gravitational coupling of size α = v dm /2π; see Eq. (5). On the other hand, α can be also interpreted as the effective circular velocity in the bound state at a distance R ⋆ , v b ≡ GM/R ⋆ = α. For R ⋆ ≃ AU, this is v b ≃ 30 km/ s and is indeed about of factor of 2π smaller than v dm ≃ 240 km/ s. Note that our approximation of the 1/r potential and treatment of the bound states break down when R ⋆ is smaller than the solar radius R ⊙ (green line), which occurs for m ≳ 2 · 10 −13 eV (darker region in Figure 3).

Phases of formation
A consequence of Eq. (4) is that the evolution of the system is characterized by two or three phases, depending on whether Γ is positive or negative. The relevant timescale tracks the relaxation time via where in the last equality we defined λ ≡ −m 2 /f 2 a , valid if the ULDM is an axion with decay constant f a . This represents the typical time a particle in a gas with density ρ dm and average square velocity v 2 dm takes to change its velocity by order one via the self-interactions mediated by λ, in the absence of external gravitational potentials [57,58]. This timescale τ rel is therefore the analogue of the gravitational relaxation time [31], with gravitational interactions replaced by the self-interactions. Even extremely weak self-interactions (λ ≪ 1), for which the probability of any single 2 → 2 scattering is small, lead a relatively short timescale, ultimately because the constants C, Γ 1 , Γ 2 in Eq. (4) are proportional to factors of the occupation number of the bosons in the galactic halo (from Bose enhancement). We summarize the evolution below.

(i) Linear growth
Since no bound state is initially present (M ⋆ = 0 at t = 0), the stimulated capture and stripping terms in Eq. (4) are negligible. As a result of direct capture only, the bound mass M ⋆ starts increasing linearly with time, M ⋆ (t) = Ct. This linear growth lasts for a time 1/|Γ|, after which the Γ 1 and Γ 2 terms become important. In Section 4.3 we show that for ξ foc ≳ 1 this critical timescale is in fact similar to the relaxation time, On the other hand, for ξ foc ≲ 1, both processes of capture and stripping are suppressed and the timescale is longer than τ rel and reads 1/|Γ| ≃ 10 τ rel /ξ 4 foc .
(ii) Exponential growth/saturation After the time 1/|Γ|, enough bound mass M ⋆ has been accumulated that the stimulated capture/stripping terms become relevant.
• For ξ foc ≳ 1, as mentioned, stimulated capture dominates over stripping and Γ > 0. As a result of stimulated capture only, the bound mass increases exponentially, i.e. M ⋆ ∝ e Γt , with an exponential timescale 1/Γ ≃ 0.3τ rel similar to the relaxation time. This leads to the formation of a 'dense' gravitational atom. In Section 4.3 and Appendix C we show that also the first few excited states are populated by DM capture at an exponential rate, but then quickly decay (on a timescale fixed approximately by τ rel ) into the ground state. It is thus a good approximation to consider only the ground state.
• On the other hand, for ξ foc ≲ 1, we find that Γ < 0 and the capture and stripping processes rapidly reach equilibrium after t ≃ 1/|Γ|, resulting in a constant bound mass. In this regime, the overdensity at the center of the atom (relative to the DM background) is which is at most a few percent of the local DM density (when ξ foc ≃ 1). The gravitational atom is therefore 'dilute'.

(iii) Collapse and Bosenova
The exponential growth of the dense gravitational atoms taking place for ξ foc ≳ 1 stops (after a few e-foldings) when the DM density in the bound state is so large that the energy in the quartic selfinteractions is of the order of the gravitational potential energy that keeps the particles bound. This occurs at the critical density which can be orders of magnitude larger than the background density ρ dm . For attractive selfinteractions, λ < 0, when the density reaches ρ crit the gravitational atom becomes unstable and collapses. A short time after collapse begins, higher-order self-interaction terms become relevant and, for an axion-like potential, the bound state is expected to undergo a Bosenova explosion, emitting an order-one fraction of its mass into relativistic ULDM particles (as does its self-gravitating counterpart, a boson star [59,60]). The phases (i), (ii), and (iii) then repeat again. On the other hand, for repulsive self-interactions, λ > 0, the density saturates to (at least) ρ crit and there is no collapse or Bosenova.
Summarizing, the role of the two ULDM parameters m and λ is the following: • The mass m sets the radius R ⋆ as in Eq. (2) and the relevance of gravitational focusing through ξ foc , as in Eq. (5), and thus determines the possibility of exponential increase of the bound states, occurring when ξ foc ≳ 1. Larger values of m lead to larger ξ foc , but at the same time to a smaller radius; they also lead to a longer capture timescale (for λ fixed) as in Eq. (6).
In the Solar System, for v dm ≃ 240 km/ s the exponential increase occurs for scalar masses in the range (1 ÷ 2) · 10 −14 eV ≲ m ≲ 2 · 10 −13 eV. The upper limit correspond to the smallest possible atom (with radius equal to the Sun's radius), for which ξ foc ≃ 10, and the lower limit to ξ foc ≃ 1, for which as mentioned R ⋆ is of order AU. However, the lower limit for m decreases for smaller v dm , as we discuss further below.
• On the other hand, λ only sets the capture timescale via the relaxation time in Eq. (6). Intuitively, the self-interactions determine how quickly the particles can change their energy by order one, but only when ξ foc ≳ 1 this energy change leads to their efficient capture.
Finally, λ also determines the critical density ρ crit in Eq. (10): For weaker self-interactions, a larger density is needed for the self-interaction energy to equal the gravitational potential energy.

Dark matter overdensity
To demonstrate the impact of the dense solar halo (ξ foc ≳ 1) on the local DM density in the Solar System, in Figure 4 we show its density profile, which follows ρ(r) ∝ e −2r/R⋆ , for different values of the boson mass between 10 −13 eV and 10 −15 eV. The density is normalized to the background density ρ dm . For each mass we show the profile in the best-case scenario, i.e. at the final stages of the exponential increase, just before the critical density is reached, choosing for each line the value λ such that this is reached at 5 Gyr (roughly the lifetime of the Sun). For smaller |λ|, the maximum density would be reached after a longer time, according to Eq. (6); for larger |λ|, within 5 Gyr multiple Bosenova cycles would have occurred, for each of which the density is smaller than in Figure 4 proportionally to 1/|λ|. The values of λ used, of order 10 −57 ÷ 10 −61 , for an axion correspond to f a within the range 10 7 GeV to 10 8 GeV for the highest and lowest masses, respectively. In the shaded regions the DM velocity parameters are varied from v dm = √ 2σ = 240 km/s to 50 km/s, with the lower edge corresponding to the largest velocity in this range (see text for details; the value 50 km/ s might be reasonable if the DM resides in a dark disk). Dashed lines correspond to the profile for m < 10 −14 eV, for which exponential growth of the bound state occurs only when v dm ≪ 240 km/s (and results are shown for v dm = 50 km/ s). The profile is shown in the best-case scenario, choosing the values of λ such that the critical density in Eq. (10) is reached at 5 Gyr (roughly the lifetime of the Sun). The black and brown dots are the constraints on the maximum DM mass that can be bound to the Sun, arising from the analysis of Solar System ephemerides using planets [61] and asteroids [62] (respectively) as test masses. For each dot, the horizontal axis indicates the radius within which such DM would be confined, and the vertical axis its maximum possible average density.
Additionally, to show the effect of a lower velocity, over the shaded bands in Figure 4 the DM velocity parameters are varied from v dm = √ 2σ = 240 km/s to 50 km/s. Although speculative, this lower value might be plausible if the DM resides in a dark disk, where it could have a smaller average velocity with respect to the Sun, and a smaller dispersion [53,63,64]. The lower/upper edges of the bands correspond respectively to the larger/smaller velocity. 3 For m ≪ 10 −14 eV, the exponential growth occurs only when v dm ≪ 240 km/s; dashed lines show the corresponding profile when the critical density is reached for v dm = 50 km/s. For comparison, in Figure 4 we show with black dots the constraints on the maximum mass that can be bound to the Sun, derived from Solar System ephemerides [44,61]. Dark matter in the orbit of a planet would contribute a force ∝ r and give rise to anomalous perihelion precession; this leads to a point-like constraint on the DM density in the orbits of each planet (see [62]). Although the density in the halo can be orders of magnitude larger than ρ dm , it does not violate these direct constraints. Note that the total bound mass is a small fraction of the solar mass for all the plotted lines.
The dilute gravitational atoms (ξ foc ≲ 1) are also novel targets for experimental searches. Although they only constitute a subcomponent of the DM, this could have a 'coherence' time τ ⋆ larger than that of the DM in the galactic halo, which is τ dm ≡ 2π/(mv 2 dm ). Unfortunately, our analytic treatment in Figure 5: The ratio (ρ/ρ dm )(τ ⋆ /τ dm ), which describes the detection reach of ULDM (linearly coupled to the SM) bound into a gravitational atom relative to that in the galactic halo, as a function of m, assuming v dm = 240 km/sec. The effective coherence timeτ ⋆ = min(τ ⋆ , t exp ) is relevant to an experiment lasting t exp ; τ ⋆ in Eq. (11) is a conservative lower bound on the coherence time of the atom. We show t exp = 1 yr, 5 yr (solid lines). For comparison, we also show the same ratio if the coherence time coincided with that of the DM galactic halo,τ ⋆ = τ dm (dashed line) and also forτ ⋆ = τ ⋆ (long dashed). Thanks to the larger coherence time, even for m ≲ 10 −14 eV, for which the gravitational atom is a subcomponent of the DM, the experimental reach could be better than that of the galactic halo DM.
Section 4 only determines the change in the captured mass over times longer than τ dm , so it does not allow to compute τ ⋆ from first principles. If the bound particles correspond to an exact energy eigenstate of the system, then τ ⋆ is in principle infinite. If this is not the case, then one could estimate a lower bound on τ ⋆ by interpreting the virial velocity v b = α as a velocity dispersion in the field, i.e.
This lower bound on τ ⋆ is larger than τ dm for all ξ foc < 2π. In the Solar System for v dm = 240 km/s, this holds for m ≲ 10 −13 eV, i.e. for essentially all the halos that can form around the Sun. Eq. (11) could be important for experimental searches based on resonance: The sensitivity of these searches to ULDM e.g. linearly coupled to the SM grows with the local DM density ρ and the measurement time t as ϕ √ t ∝ √ ρ t, until t approaches the coherence time (see [65,66]). In Figure 5 we illustrate the ratio (ρ/ρ dm )(τ ⋆ /τ dm ) between the sensitivity reach for a gravitational atom and the DM in the galactic halo; here,τ ⋆ ≡ min(τ ⋆ , t exp ) is an effective coherence time relevant to a search lasting t exp . For m ≳ (1 ÷ 2) · 10 −14 , where exponential growth occurs, there is a large enhancement of the sensitivity, from both a larger density and a larger coherence time. In addition, even for m ≲ 10 −14 eV, where the growth saturates to ρ < ρ dm as in Eq. (9), the detection reach could be better than that of the virialized halo, owing to the larger coherence time.

Implications
The results above imply that, after τ rel , dense gravitational atoms are expected to form around all the massive objects for which the condition ξ foc ≳ 1 is satisfied (including, for instance, solar-mass stars and black holes), which act as 'seeds' for DM capture. For astrophysical systems with M > M ⊙ or smaller DM velocity, the exponential increase occurs also for m smaller than 10 −14 eV, for which τ rel is shorter.
Over a time of the order of τ rel , dynamical relaxation of the DM in the galactic halo via self-interactions is likely to lead to additional effects. Similarly to what happens to fuzzy dark matter (with m ≃ 10 −21 ÷ 10 −22 eV) via the gravitational interactions alone [34,60], after τ rel a solitonic core is expected to form in the galactic center. (This would occur even faster than the formation of the solar halo, as the DM abundance close to the galactic center is expected to be larger, corresponding to a smaller τ rel .) For the masses m of interest, such an object may be absorbed by the supermassive black hole [67,68]. Additionally, self-gravitating boson stars may randomly be produced anywhere in space by relaxation of the scalar field to its ground state via self-interactions [34,57]. In contrast to the gravitational atoms around external gravitational potentials considered here, their accretion rate does not follow Eq. (4) and their density is not expected to increase exponentially (see also [69]). Thus, these bound objects should accrete more slowly.
Sufficiently strong ULDM self-interactions alter the cosmological evolution of DM perturbations. As we discuss in Section 7, attractive/repulsive self-interactions tend to enhance/suppress the perturbation growth and can therefore be constrained from measurements of the linear matter power spectrum [70][71][72]. In any case, as shown in Figure 14, the values of λ required for the formation of the atom to take place within 5 Gyr are consistent with these constraints over a few orders of magnitude. We will also show that these constraints may become considerably stronger if one considers the evolution of the field before matter-radiation equality, an effect not pointed out in previous studies.
Let us now briefly comment about the relic abundance and axions. The simplest version of the misalignment mechanism underproduces the DM for the parameters m ≃ 10 −13 ÷ 10 −14 eV and f a ≲ 10 8 GeV that lead to the formation of the halo around the Sun within the age of the Solar System (equivalently, the expected size of the self-interactions is too small to have τ rel ≃ Gyr). This is because, for a conventional axion potential, f a both sets λ and the axion field range, of order f a . However this is not necessarily the case for scalars with more generic potentials, where in principle the field range is not related to the self-couplings and can extend up to the Planck scale, or in the presence of monodromies [73]. For instance, a scalar field with potential containing the mass term and the self-interaction term in Eq. (3) could lead to the observed relic abundance via misalignment, provided that the initial field value is chosen appropriately. In any case, there is a variety of other production mechanisms leading to a cosmologically nontrivial axion relic density (see for instance [74][75][76]), some of which could lead to the correct DM abundance for the formation to take place over Gyr timescales.
The axion coupling to photons is usually parameterized by the Lagrangian L ⊃ 1 4 g ϕγγ ϕF µνF µν , with g ϕγγ ≡ c ϕγγ α em /(2πf a ), where F µν is the photon field strength, α em is the electromagnetic finestructure constant, and c ϕγγ is a model-dependent constant. For masses m ≲ 10 −13 eV, the strongest bounds currently come from helioscope measurements [77] and several astrophysical observations (see e.g. [78][79][80]) as g ϕγγ ≲ 10 −10 ÷ 10 −12 GeV −1 . 4 Fast enough capture requires f a ≲ 10 8 GeV, so that the unknown constant c ϕγγ needs to be somewhat smaller than O(1). Additionally, as discussed further in Section 6.1, SM couplings can in principle disrupt the formation of the bound state. However, their effect is not amplified by Bose enhancement, and so weak and not expected to affect the formation.
Finally, in Figure 6 we compare, as a function of m, the DM overdensity of the solar halo with that occurring if an order-one fraction of the DM is in boson stars (as a boson star passes through the Solar System). As illustrated, the encounter rate γ of boson stars with the Earth is of experimental relevance (at least one encounter per ∼ few years, with overdensity of at least unity) only for relatively large particle masses m ≳ 10 −9 eV (see for instance [45] for details). On the other hand, a gravitational atom formed around the Sun gives rise to a persistent overdensity, making it a prime target for experimental searches.
The results in this paper open up several signatures and opportunities for discovery of (or constraints on) ultralight dark matter, which we briefly discuss in Section 8 and will explore in future work.

Bound states and waves
Consider a scalar field ϕ with potential V (ϕ) (we will use the terms 'scalar' and 'axion' interchangeably, because in what follows nothing will depend on the parity of the field). We are interested in cases where the amplitude of ϕ is small compared to the typical field range, as expected in our virialized DM halo. Thus, it is appropriate to expand V around its minimum as where the dots include higher orders in the expansion. 5 The sign of λ (negative or positive) determines whether the quartic self-interactions are attractive or repulsive. In the following it will be useful to express λ in terms of a dimensionful coupling g (to be used interchangeably with λ) as If ϕ is an axion with mass arising from weakly-coupled instantons [12,101], its potential is a periodic function dominated by a single harmonic, V (ϕ) = −m 2 f 2 a cos(ϕ/f a ) . In this case, the attractive self-interactions have g = −1/(8f 2 a ) and −λ ≃ 10 −61 (m/10 −14 eV) 2 (10 8 GeV/f a ) 2 . However, the only requirement on V (ϕ) is that it can be approximated as in Eq. (12).
Given its large occupation number, the DM follows the classical equations of motion (EoM) corresponding to the potential of Eq. (12), i.e.
where the metric is g 00 = 1+2Φ, g ij = −(1−2Φ)δ ij , D µ is a covariant derivative, and Φ is the gravitational potential. The (general-)relativistic corrections will be small in the following. Thus, the EoM can be rewritten in a nonrelativistic form, i.e. working in terms of a nonrelativistic field ψ defined by in the limit where ψ is slowly varying,ψ ≪ mψ. In terms of ψ, the EoM reduces to the Gross-Pitaevskii (GP) equation (i∂ t + ∇ 2 /2m − mΦ)ψ = g|ψ| 2 ψ. In the vicinity of a massive object of mass M the gravitational potential is well approximated by that generated by the object itself, i.e. Φ = Φ ex ≡ −GM/r, with r ≡ |x|. 6 As a result, the GP equation simplifies to The gravitational coupling α ≡ GM m measures the strength with which the particles are pulled towards the potential well. Eq. (15) is nonlinear due to the self-interaction term.
In the nonrelativistic limit, the mass density of the field is ρ = m|ψ| 2 , and the number density of particles is n ≡ ρ/m = |ψ| 2 . The (nonrelativistic) energy density, defined as ϵ ≡ T 00 − ρ, with T µν being the stress-energy tensor of ϕ, is given by see Appendix A for the derivation. The terms in the decomposition of Eq. (16) correspond, respectively, to the kinetic, gravitational potential, and self-interaction energies. The latter is negative for attractive self-interactions, as is the gravitational potential energy. If the density ρ is (locally) smaller than a critical density ρ crit ≡ 2|Φ ex |m 2 /|g|, the self-interaction energy is negligible compared to that of gravity (at r ≃ R ⋆ , this reproduces Eq. (10)). In this limit the interaction term in the right-hand side of Eq. (15) can be treated as a perturbation. At the zeroth order, the EoM are the same as those of the hydrogen atom, with α and the number density |ψ| 2 playing the roles of the fine-structure constant and the quantum-mechanical probability density, respectively.
In analogy with the hydrogen atom, Eq. (15) admits quasi-stationary time-periodic solutions of the form ψ = e −iωt Ψ(x), where Ψ solves the eigenvalue equation (−∇ 2 /2m + mΦ ex )Ψ = ωΨ. Their energy is E ≡ d 3 xϵ = N ω, where N ≡ d 3 x|ψ| 2 is the total number of particles, so that ω has the interpretation of energy per particle. Such solutions can be divided into two classes: bound and unbound states of the external potential, depending on whether ω (or E) is negative or positive.

Bound solutions
Bound-state solutions of Eq. (15) are the equivalent of hydrogen atom orbitals. They are discrete states, Ψ ∝ ψ nlm , labelled by the integers n, l, m and have negative (binding) energy ω n ≡ −mα 2 /(2n 2 ). 7 The normalized lowest-energy solution (i.e. the ground state) is This corresponds to a spherically-symmetric density field, which is maximized at the center r = 0. The parameter R ⋆ , given in Eq. (2), is the typical radius. We refer the reader to Appendix A for the form of the excited states, where for example e −r/R⋆ → e −r/(nR⋆) . The bound mass is M ⋆ ≡ d 3 xρ = m d 3 x|Ψ| 2 , so that Ψ = M ⋆ /m ψ 100 if only the ground state is populated. In contrast to the hydrogen atom, the occupation number and therefore M ⋆ can be macroscopic because of the bosonic nature of the particles. In particular, in the limit of g → 0 (and of negligible self-gravity), the mass can grow arbitrarily (until relativistic corrections become relevant). Note that R ⋆ only depends on m and M , and not on M ⋆ , which enters only in the normalization of Ψ. 8 The gravitational atom can be interpreted as a collection of particles bound to the external body, with typical 'virial' velocity v b = GM/(n 2 R ⋆ ) = α/n, set by α (this is consistent with the semi-classical virial relation |ω 1 | ≈ mv 2 b /2).

Unbound solutions
The unbound solutions are traveling waves, which correspond to the Coulomb scattering states. They form a continuum of states, Ψ ∝ ψ k , parametrized by a three-vector k that represents their (asymptotic) momentum, and have positive energy ω k ≡ k 2 /2m. The normalized solution reads [53,103] where Γ[a] is the Gamma function and 1 F 1 [a, b, c] is the confluent hypergeometric function. 9 For later convenience, in Eq. (18) we expressed m and α in terms of R ⋆ = (mα) −1 , although these solutions do not depend on bound-state properties per se. The momentum k of the waves can be also written in terms of their velocity v (with respect to the body) as k = mv. Importantly, m and α appear in Eq. (18) only through As mentioned in the Introduction, this measures the ratio between the de Broglie wavelength 2π/k and R ⋆ . Given than α = v b for the ground state, this is also ∝ v b /v, and therefore controls the relative importance of the virial velocity and that of the wave, modulo a 2π factor. The limits ξ foc (k) ≪ 1 and ξ foc (k) ≫ 1, shown in Figure 2, correspond respectively to large/small R ⋆ compared to 2π/k, or waves that are fast/slow with respect to the virial velocity (v ≪ 2πα and v ≫ 2πα). In these two regimes the scattering states in Eq. (18) have simpler expressions.
• If ξ foc (k) ≪ 1, ψ k in Eq. (18) reduces to a plane wave everywhere, as if no potential were present: In this limit, the density field ∝ |ψ| 2 is homogeneous, with (maximum) value independent of k (see blue line in Figure 7, left).
• If ξ foc (k) ≫ 1, the plane wave is distorted around the body at distances smaller than the wavelength, r ≲ 2π/k. Within this region, an excellent approximation of Eq. (18) is where J 0 is the spherical Bessel function. 10 Instead, for r ≫ 1/k, ψ k → e ik·x and recovers its plane wave form, at least in the region wherek ·x < 1. 11 To better understand the spatial distribution of Eq. (21), in Figure 7 (right) we plot the density normalized to that at the center, where it is maximized, and decays as J 0 [a] → sin a/a for a → ∞. As a result, the maximum density occurs in the line r(1 −k ·x) = 0. This corresponds to a 'tail' of overdensities behind the body along the direction ofk.
The absolute value k of the momentum appears only as an overall multiplier in Eq. (21). Thus, in this limit k only controls the maximum density (as opposed to its full spatial distribution), which increases as |ψ k (0)| 2 = 2π/kR ⋆ = ξ foc (k) at small momenta, as shown by the orange line of Figure 7 (left). Additionally, r only enters as the ratio r/R ⋆ . Given that J 0 [a] changes by order one between a = 0 and a ≃ 1, the amplification of the amplitude occurs around the origin over a region always of size R ⋆ , irrespective of k; see Figure 7 (right). This amplification is the gravitational focusing discussed in Section 2. Importantly for the growth of the bound state in Section 4, the amplified region tracks the size of the gravitational atom for all k in this limit.
The transition between the two regimes is encoded into the hypergeometric function of Eq. (18). A good measure to distinguish between them is the maximum value of |ψ k | 2 , which always occurs at r = 0 also for the full scattering states in Eq. (18) and reads This quantity is shown in Figure 7 (left) together with its limits |ψ k | 2 → 1 and |ψ k (0)| 2 → ξ foc (k) in the two regimes. As mentioned, the turning point between these occurs at around ξ foc (k) = 1. We stress that this corresponds to the wavelength 2π/k being equal to R ⋆ , rather than the velocity v equal to α.
In fact, waves with velocity v = 2πα still get distorted at the order-one level, even though one could naively guess that they are not affected because their speed exceeds the escape velocity at r = R ⋆ (which is √ 2α).

DM in the galactic halo
In the vicinity of the Sun, the DM in our galactic halo is a superposition of scattering states. In the reference frame of the Sun, the field takes the form where ω k = k 2 /2m, ψ k is as in Eq. (18) and a(k) encodes the statistical properties of the DM distribution. Over timescales longer than the 'coherence' time τ dm ≡ 2π/mv 2 dm , a(k) can be effectively treated as a random variable [24]. Its norm describes the DM number density present at each momentum k and, at a given k, it follows a Rayleigh distribution. Its phase should be uniformly distributed in [0, 2π). In particular, if we indicate by ⟨·⟩ the statistical average, which corresponds to averaging over times much longer than τ dm , we expect ⟨a(k)⟩ = 0 and More generally, all such averages vanish unless the same number of a and a * appear (so that the uniformly distributed phases cancel), and in this case they can be reconstructed from Eq. (24) via Wick's theorem [104]. The function f (k) in Eq. (24) is the DM occupation number at momentum k and can be connected to local background DM density ρ dm by requiring that ρ dm = ⟨ρ w ⟩ for |x| ≫ 2π/k (i.e. far from the Sun), which leads to where in the last equality we used the fact that |ψ k (x)| 2 → 1 in the limit |x| ≫ 2π/k. Thus f (k) can be also seen as the average DM number density per phase space volume, with corresponding number of An approximate form of f (k) is a Gaussian with average momentum mv dm ≡ −mv ⊙ and variance (mσ) 2 , where v ⊙ is the Sun's velocity with respect to the galactic rest frame [105,106] and σ is the DM velocity dispersion in the same frame. In other words, We will assume v ⊙ = v dm = 240 km/s (see e.g. [105,106]). In the standard virialized halo model, the density near the Sun is ρ dm ≃ 0.4 GeV/cm 3 [18][19][20][21], and the DM velocity dispersion is σ = v c / √ 2 where v c ≃ 230 km/ s is the circular velocity of the Milky Way at the position of the Sun [106]. Note that although v ⊙ slightly differs from v c (because the Sun has a small relative velocity with respect to the average Milky Way motion), we will assume in reality, σ is slightly smaller than this). In any case, the DM distribution is by no means precisely known, and the above should therefore be treated as a benchmark only; see e.g. [22,23] for further discussion.
In the following, as in Eq. (5), we define the typical focusing parameter in the Solar System as i.e. when ξ foc does not appear with an argument, it is evaluated at k = mv dm .

Bound state formation: analytic approach
In the absence of self-interactions, the dark matter of our galactic halo is described the superposition ψ w of unbound solutions in Eq. (23). If g ̸ = 0, ψ w is however only an approximate solution of the EoM, Eq. (15). We now show that, starting from initial conditions provided by ψ w , the self-interactions lead to capture processes resulting in the growth of the bound-state population: The final result describing this evolution is given in Eq. (41). This is derived in two ways: In Section 4.1 by solving Eq. (15) perturbatively in the self-interaction term, and in Section 4.2 using a much simpler quantum-mechanical scattering calculation.

Perturbation theory
As long as g|ψ| 2 < mΦ ex (i.e. ρ < ρ crit ), the spectrum of bound states is well-approximated by ψ nlm . In this regime it makes sense to calculate the rate of changeṄ nlm (t) in the number of particles bound to the nlm level, by solving Eq. (15) perturbatively. At the end of this subsection we discuss the limitations of this calculation.
Since the bound and scattering states ψ nlm and ψ k constitute a complete basis, it is convenient to expand the field as In this way, the total field can be interpreted as a superposition of bound and wave components: The number of particles bound to the nlm level is the bound part of d 3 x|ψ| 2 and at the leading order in the perturbative expansion is simply N nlm (t) = |c nlm (t) is obtained by solving Eq. (15) at first order in g eff , which reads This is a source equation for ψ (1) : Via their self-interactions, the waves induce a smaller field with both bound and wave components. Diagrammatically, this field can be interpreted as arising from 2 → 2 scattering processes mediated by the quartic self-interactions. Indeed, a generic 'source' term on the right-hand side of Eq. (29) has the form gψ Such a term can be associated to the process i + j → k + ind.
Specifically, substituting ψ (0) = ψ w into Eq. (29), multiplying by ψ * nlm e iωnt and using the othonormality of ψ nlm , we obtain that the bound-state components are generated with coefficients where we used the short-hand notation [dk] ≡ d 3 k/(2π) 3 , and defined the 'matrix element' for the 2 → 2 process k 1 + k 2 → k 3 + nlm as Note that c nlm grows, and the corresponding state gets populated, only for momenta for which the energy difference vanishes, i.e. if energy is conserved in this scattering. The relation ⟨a Eq. (24) allows us to easily compute the average ⟨Ṅ nlm ⟩ over times longer than the DM coherence time τ dm . Wick's theorem leads to where we used the identity d dt | t 0 dt ′ e i∆ωt ′ | 2 = 2πδ(∆ω) valid in the limit t ≫ 1/∆ω. It is useful to represent the results in Eqs. (30) and (33) via the following diagrams: In the first two diagrams, incoming (outgoing) lines represent incoming (outgoing) states in the process k 1 + k 2 → k 3 + nlm and are related to the factors e −iω i t (e iω i t ) in Eq. (30). Solid lines indicate the (three) factors of the coefficients a(k) (and their conjugate for outgoing states), which come from the zerothorder field from the right-hand side of Eq. (29). On the other hand, dotted lines represent the remaining 'sourced' state. The last diagram can be thought as the 'product' of the first two after the statistical average of Eq. (24). In this average, an incoming solid line in the first diagram has been 'contracted' with an outgoing solid line in the second diagram: This provides a double line, associated with the occupation number factor f (k) in Eq. (33). The two ways of performing such a contraction give a factor of 2 for the last diagram, as in Eq. (33).

Stimulated capture and stripping
Eq. (33) shows that the number of bound particles grows (linearly) in time via the the capture one of two scattered DM particles. When enough particles have been accumulated in the bound states, they could stimulate further capture through Bose enhancement, and be simulteneously depleted via scattering with the background DM waves. To understand these effects, we consider a generic later time t = t 0 > 0, when the bound states are populated with occupation number N (0) nlm . At this time the unperturbed solution takes the form representing the bound component. As before, we write ψ = ψ (0) + ψ (1) + ψ (2) + . . . and determine the rateṄ nlm (t 0 ) by solving Eq. (15) perturbatively in a time interval t 0 < t < t 0 + ∆t short enough that the perturbation of the field is smaller than ψ (0) . Note that in regions where the bound-state density dominates, the perturbative parameter is g eff = gρ/(mω) ≃ ρ/ρ crit (and is smaller than unity by assumption), since ω ≃ mα 2 is the typical field's frequency in the region. Using Eq. (28), N nlm (t) reads where in the last expression the dots stand for terms of order g 3 eff . The change in N nlm has, as before, a positive contribution from |c (1) nlm | 2 , of order g 2 eff , which can be read from Eq. (33). However, an additional time-dependent term -the second in the right-hand side of Eq. (36) -appeared from the 'interference' of the induced field with the initial bound component. This term is of order g eff . In any case, c (1) nlm ∝ aaa * from Eq. (30), and if we are interested only in the variation of N nlm over times longer than the DM coherence time τ dm , the contribution proportional to c (1) nlm vanishes given the relation ⟨a * (k 3 )a(k 1 )a(k 2 )⟩ = 0, see Section 3.1. As a result, the full leading contribution to ⟨Ṅ nlm ⟩ arises at second order and reads The second term of Eq. (37) requires calculating the field's perturbation at order g 2 eff , which follows This can be solved similarly to the first-order equation, as explained in Appendix B. We give here only a diagrammatic representation of the result. Eq. (38) is a source equation for ψ (2) , which is generated by two powers of ψ (0) and one power of ψ (1) . Thus, c (2) nlm is represented by a 2 → 2 diagram where one of legs contains the first-order field ψ (1) , and, similarly to Eq. (34), this is incoming/outgoing depending on whether ψ (1) appears without/with the complex conjugate in the right hand side of Eq. (38). The non-vanishing terms are: These are generated from the first and second terms in Eq. (38) respectively. In the first (second) diagram, an outgoing (incoming) wave marked with , which is produced by a 2 → 2 scattering processes involving the bound state nlm, discussed in Appendix B.
As before, solid lines represent factors of a(k) or N nlm . 13 One can also view these diagrams as a combination of two first-order sub-diagrams connected by a 'propagator'. When calculating the average ⟨c (2) nlm ⟩, an outgoing wave is contracted with an incoming one and in the first diagram there are two choices for this. After the contraction, the propagator is 'cut' and provides the square of the amplitude M in Eq. (31), in similar fashion to the optical theorem. A straightforward derivation, see Appendix B, gives where, as before, double lines are associated to f (k) or N (0) nlm . By collecting the contributions in Eqs. (33) and (40) and noting that N where N nlm is evaluated at a generic t = t 0 and the second and third terms come from Eq. (40), where indeed each diagram is proportional to two powers of f and one of N  (33). Thus, at any time until when the critical density is reached, the evolution follows Eq. (41). We dub the three terms in this equation 'capture', 'stimulated capture', and 'stripping', respectively, and summarize the corresponding physical processes in the diagrams in • In the capture process, diagram (a) in Figure 8, one of the incoming particles transfers enough energy to the other that its velocity decreases and it gets trapped into the gravitational potential, therefore becoming a bound state. The capture term is independent of particles already present in the bound states and thus constitutes an unavoidable positive contribution to the bound state occupation (this corresponds to the first term in Eq. (4)). Additionally, it is proportional to the occupation numbers in the initial state, f (k 1 ) and f (k 2 ), but also to that in the final state f (k 3 ), given the three powers of a in Eq. (30).
• Conversely, the stimulated capture and stripping term, diagrams (b) and (c) respectively, arise from the interference term in Eq. (36) and vanish if no bound state is initially present. Thus, they are contributions 'stimulated' by the (same) bound state being already populated, and are indeed proportional to N nlm . Stimulated capture contributes positively to the bound-state occupation, whereas stripping contributes negatively and corresponds to the depletion of a bound state by the scattering of a wave, producing another wave with smaller energy. This term has a multiplier of 2 relative to the others because diagram (c) provides a double contribution, by exchanging k 1 ↔ k 2 , as in Eq. (40).
Observe that, if only the ground state is populated, i.e. N (0) nlm = 0 for nlm ̸ = 100, stimulated capture and stripping would be irrelevant for all levels that are not 100.
Since ⟨Ṅ nlm ⟩ ∝ g 2 , the efficiency of these processes, and the resulting changes to the number of bound particles, do not depend on whether self-interactions are attractive or repulsive.

Expected validity
Eq. (41) describes the instantaneous variation in the number of bound particles, averaged over times longer than the DM coherence time τ dm . For simplicity of presentation, in Eq. (41) we only wrote the contributions toṄ nlm that are linear with N nlm . In reality, the changes in occupation number of different nlm levels are not independent from each other, because of additional (more complicated) terms that appear in the right-hand side of Eq. (41). These represent the exchange of particles from one bound state to another, which can be positive or negative, and are nonlinear with N nlm . They become more important as the number of bound particles grows, leading in principle to an infinite set of coupled equations for all the levels. For instance, these include excitation and relaxation terms, where bound particles transition to higher/lower nlm states, respectively. In Appendix C we discuss the effect of these terms by studying a simplified two-level system consisting of two bound states, nlm = 100 and 200, and show that, if populated, the higher level does not deplete the 100 state and can in fact enhance it through de-excitation, e.g. k 1 + 200 → k 2 + 100. Therefore the captured density described in what follows may be a lower bound to the true result.
Note that the perturbative approach requires the corrections to the energy levels, typically of order gρ/m, to be much smaller than the unperturbed ones, ∼ mα 2 (in Appendix B we give a precise derivation of such corrections). As anticipated, this holds when ρ ≪ m 2 α 2 /g ≃ ρ crit , and as a result this requirement is equivalent to g eff ≲ 1.

Quantum scattering and Bose enhancement
It is interesting to note that the change in the number of bound particles in Eq. (41) can be understood as arising from 2 → 2 scattering processes at the quantum level, upon adding the expected Bose enhancement factors. First, observe that the action of the scalar field in the nonrelativistic limit reduces to where ψ is the propagating degree of freedom, and Φ ex = −GM/r is the external background field. In particular, the EoM in Eq. (15) immediately follows by extremizing the action of Eq. (42). The action S can be quantized by following standard canonical quantization for a non-relativisitic field theory. The main difference is that the quadratic part of S implies that the single-particle states are bound states, with wave functions ψ nlm described by Eq. (70) of Appendix A, or waves ψ k described by Eq. (18), which at large distance behave as plane waves, e ik·x . These single-particle states are obtained by acting with the corresponding creation operators on the vacuum |0⟩ as |nlm⟩ = a † nlm |0⟩ and |k⟩ = a † k |0⟩. (The free-theory field operatorψ is constructed from such creation and annihilation operators, with coefficients given by ψ nlm e −iωnt etc.) The process where T is the Dyson time-ordering operator, M is as in Eq. (31), and we used and similarly for ⟨k 3 nlm|ψ †ψ † . The probability per unit time of the process k 1 + k 2 → k 3 + nlm is obtained by squaring A (stripping off 2πδ(∆ω)), and reads (2π)δ(∆ω)4g 2 |M| 2 .
To calculate the total rate of change in the particle number in the nlm state, one needs to subtract the inverse process, which has (in modulus) the same matrix element, multiply by the occupation number of the particles in the initial state, and integrate over all possible values of the momenta. This leads to conclude thatṄ nlm solves a 'Boltzmann' equation, which effectively describes the rate of change of particle-number per phase-space volume, and equals where a factor of 1/2 has been added in order to avoid double-counting from the exchange of the identical particles 1 and 2. The factors f + 1 and N nlm + 1 have been inserted in the final states in Eq. (44) because of Bose enhancement: The transition rate of indistinguishable particles P a→b from the state a to b is given by (N b + 1)P a→b , where P a→b is the rate for distinguishable particles.
Owing to the Bose-enhancement factors, at the leading order in f and N nlm Eq. (44) coincides with Eq. (41), and reproduces the capture, stimulated capture, and stripping terms derived from the classical EoM. In particular, the presence of a macroscopicṄ nlm can be understood as resulting from the Bose enhancement of identical quantum-mechanical particles.

Dilute vs dense gravitational atoms
We can multiply Eq. (41) by the boson mass to obtain the instantaneous rate of change in the mass M nlm bound to the nlm level, averaged over times longer than the coherence time. This readṡ where the dots represent additional terms (not reported in Eq. (41), and not crucial for the main point of this discussion for the 100 level) that are nonlinear with M nlm and describe the exchange of particles between different levels; see Appendix C for further details. The total bound mass is (45) reproduces and generalizes Eq. (4) of Section 2, with the nlm-dependent coefficients C > 0 and Γ ≡ Γ 1 − Γ 2 given by with M and ∆ω = ω k 1 + ω k 2 − ω k 3 − ω n as in Eqs. (31) and (32).
The matrix element M is, and therefore C and Γ are, independent of the bound mass M ⋆ . Thus, the solution of Eq. (45) is straightforward: As mentioned in Section 2, at early times the bound mass increases linearly, M nlm (t) = Ct (dashed line in Figure 9), with C representing the mass captured per unit time. At t ≃ 1/|Γ| the bound-state-stimulated processes become relevant. If Γ < 0 the bound mass saturates to the constant value M eq ⋆ ≡ C/|Γ| (blue line in Figure 9), while if Γ > 0 it increases exponentially as M nlm ∝ e Γt (purple line), starting from the value C/Γ. 14 C and Γ depend on the gravitational coupling α (e.g. through ψ k and ψ nlm ) and on the DM velocity and dispersion, v dm and σ (through f (k)). We now focus on the ground state and calculate these two quantities in the limits ξ foc ≡ 2πα/v dm ≪ 1 and ξ foc ≫ 1, where the integrals in Eq. (46) are tractable. In the latter limit 1/Γ is positive and will track the relaxation time introduced in Section 2 (reproduced here for simplicity), Eq. (47) can be easily derived by calculating the inverse of the rate of the process k 1 + k 2 → k 3 + k 4 , similarly to Section 4.2. The result is that 1/τ rel ≃ n g Σv dm ×f b , where n g = ρ dm /m is the number density of the gas, Σ ≃ m 2 g 2 is the cross section of the 2 → 2 scattering mediated by the self-interactions, and f b ≃ n g /(mv dm ) 3 is the occupation number enhancement factor. Some of the nlm ̸ = 100 states are also populated, but, as discussed below and in Appendix C, they are either subdominant or likely to decay to the ground state (when M nlm is large enough) over a time comparable to τ rel , via the terms in the dots in Eq. (45), and lead to a faster increase of the 100 level. Figure 9: The time evolution of the mass M nlm bound to the nlm level (neglecting exchange between energy levels), which follows Eq. (45). The bound mass initially grows as M nlm = Ct (dashed line), as a result of the capture process; this behaviour is valid at t ≪ 1/|Γ|. If Γ < 0, which occurs in the limit ξ foc ≪ 1, the mass saturates at the time t ≃ 1/|Γ| to the asymptotic value C/|Γ| (blue line). If Γ > 0, which occurs in the opposite limit ξ foc ≫ 1 due to to gravitational focusing and the dominance of the stimulated-capture term, the linear growth is instead superseded by the exponential growth M nlm ∝ e Γt (purple line). When M nlm is large enough, excited states are likely to decay into the 100 level via terms in the dots of Eq. (45), see Appendix C.
Dilute gravitational atoms: ξ foc ≪ 1 We consider the standard DM halo with v dm = √ 2σ (see Section 3.1). Although the integrals in Eq. (46) extend over all values of k 1 , k 2 and k 3 , the occupation number f (k) in Eq. (26) is a peaked at mv dm with dispersion mσ. Thus, the dominant contribution to C and Γ comes from momenta of order mv dm . For C, this is ensured by the factors of f alone, while for Γ this is also a consequence of the δ function, which sets all the k i to be of the order mv dm in the limit 2πα/v dm ≪ 1. 15 In other words, in this limit the processes in Figure 8 primarily involve waves with momenta of order mv dm , because the binding energy of the bound states, −mα 2 /2n 2 , is negligible with respect to the kinetic energy mv 2 dm /2 of the DM waves. Consequently, as discussed in Section 3, it is a good approximation to assume that all scattering states ψ k -that enter in Eq. (46) through M -are plane waves e ik·x , and the matrix element M in Eq. (31) can be computed analytically. This matrix element represents the spatial superposition of three plane waves with a bound state, and for the 100 state it reads This depends on the momenta only through the combination ∆kR ⋆ = ∆v/α, where ∆k ≡ |k 1 +k 2 −k 3 | = m∆v is the modulus of the momentum 'transferred' from the waves to the bound state. Over most of the integration region in Eq. (46), ∆v is of order v dm . Thus, ∆kR ⋆ ≃ v dm /α ≫ 1 (since ξ foc = 2πα/v dm ≪ 1) and the matrix element in Eq. (48) is suppressed in this limit. This suppression can be understood by the rapidly oscillating phase ∆k ≫ R −1 ⋆ of the plane waves, which averages out in Eq. (48) when integrated over the effective domain fixed by the ground-state radius R ⋆ . The matrix elements of higher-n states are even more suppressed, with asymptotic behaviour at ∆kR ⋆ → ∞ given by M k 1 +k 2 →k 3 +n00 → 8 πR 3 ⋆ /n 3 /(∆kR ⋆ ) 4 . (For these states the cancellation from the oscillating phase is stronger because of the larger region of support for the spatial integration.) 15 More explicitly, energy conservation requires k3 = k 2 As a result of this suppression, in the limit 2πα/v dm → 0 all the processes of capture, stimulated capture and stripping are inefficient and both the C and Γ are suppressed. Intuitively, the particles in the DM waves are so fast (relative to those in the bound state) that it is improbable that they lose enough energy via 2 → 2 scatterings to acquire a velocity of order of the ground-state virial velocity v b = α, and thereby become trapped. This is even more convincing for higher-n states, whose virial velocity v b = α/n is smaller. Therefore, in the following we will consider only the 100 level; see Appendix C for the effect of the excited states.
The coefficients C and Γ can be calculated by rewriting Eq. (46) in terms of the momenta normalized to the inverse ground-state radius, p i ≡ k i /(mα), and switching to spherical coordinates p i = p ir +θ iθ +ϕ iφ .
Choosing v dm to point e.g. towardsẑ, Γ simplifies to In the previous equation ∆p ≡ |p 1 + p 2 − p 3 | and it is understood that p 3 = p 2 1 + p 2 2 + 1, which arises from the integration over p 3 due to the δ(∆ω) factor. 16 Additionally, the integral over ϕ 3 has been performed explicitly, 17 so it is understood that ϕ 3 = 0 in the integrand of Eq. (49).
In the limit 2πα/v dm ≪ 1, the terms dependent on p i in the exponential factors of Eq. (49) are negligible, because the exponentials are small unless p i ≃ v dm /α ≫ 1, but for such large values of p i the matrix element suppresses the integral (as mentioned before, for typical momenta in the DM distribution the matrix element is small). Thus, the two addends in Eq. (49) coincide except for a factor of −2, ultimately due to the multiplicity of the diagram in Figure 8(c). As a result, in this limit Γ 2 ≃ 2Γ 1 , and stripping dominates over stimulated capture. At the leading order in α/v dm , the integral in Eq. (49) is a constant and can be estimated as −A(2π) 2 2 3 ; a numerical evaluation gives A ≃ 0.21. The corresponding saturation time is (The benchmark value of m is smaller than 10 −14 eV, which is more appropriate in this regime.) Due to the ξ 4 foc ≪ 1 factor at the denominator, it takes a time parametrically longer than τ rel to reach the steady-state regime (roughly a factor of 10 4 longer, at the benchmark above). In Figure 10 (left) we show the full calculation of Γ in Eq. (49), where we can see the (small) corrections to Eq. (50).
Similarly, C reads where the integral can be estimated as B(2π) 2 2 3 with B ≃ 0.08. The mass at equilibrium is M eq ⋆ = −C/Γ ≃ 0.4 · 8π 3/2 ρ dm /(mv dm ) 3 with central density ρ eq (0) = M eq ⋆ /πR 3 ⋆ ≃ 0.4ρ dm ξ 3 foc /π 5/2 , as in Eq. (9). The bound state only constitutes a small fraction of the DM density, which can be also understood 16 We used the identity δ(p 2 1 + p 2 2 − p 2 3 + 1) = (2p3) −1 s=±1 δ p3 + s p 2 1 + p 2 2 + 1 . 17 We took advantage of the freedom in choosing e.g. the x-z plane to coincide with the p3-z plane (the y axis is fixed by this choice, and all remaining integrations are therefore nontrivial).  47). If ξ foc ≪ 1, Γ is negative and suppressed, so the bound mass saturates over a time 1/|Γ| ∝ (α/v dm ) −4 that is parametrically longer than τ rel ; see Eq. (50). If ξ foc ≫ 1, 1/Γ is positive and coincides with τ rel : the bound mass increases exponentially with timescale ≃ τ rel . Right: The central density of the atom in units of the background density at the time 1/|Γ|, when it reaches its equilibrium value for ξ foc ≪ 1 or starts increasing exponentially for ξ foc ≫ 1. In the first case, the equilibrium density is ρ eq (0)/ρ dm ∝ (α/v dm ) 3 and is a subpercentage of the DM, see Eq. (9). In the second case, it is of order one, except for very large α/v dm , where it is suppressed because the direct capture process is kinematically disallowed. These small values will in any case exponentially increase in time with rate Γ. The behavior of the system in transition region ξ foc ≃ 1 (dotted lines) is uncertain, and studied with numerical simulations in Section 5. We show with dotted lines naive interpolations of the two regimes, which likely overestimate 1/|Γ|. from the fact that the bound mass M eq ⋆ is similar that of the background DM contained in a de Broglie wavelength-sized volume λ 3 dB , but R ⋆ ≫ λ dB in this limit. In Figure 10 (right) we show the overdensity with the full calculation of the integrals, which matches well with the estimate in Eq. (9).
Since the ratio between the capture and stripping terms is independent of g, the overdensity at equilibrium is only a function of α/v dm . On the other hand, as noticed in Section 2, the timescale to reach equilibrium |Γ| −1 is strongly dependent on g; see Eq. (50). As we show in Appendix C, higher levels are also populated at much later times. The mass bound to the higher levels is the same (when they saturate), which means that the density at the center and at r ≃ R ⋆ will be dominated by the 100 level (because the excited states have radius larger than R ⋆ ) so that the estimate in Eq. (9) is a good approximation of the full density.
Dense gravitational atoms: ξ foc ≫ 1 In the limit ξ foc = 2πα/v dm ≫ 1, the processes of capture and stripping of Figure 8, i.e. C and Γ 2 in Eq. (46), are kinematically suppressed, but stimulated-capture (Γ 1 ) is not. Indeed, given the density factors in Eq. (46), C obtains its dominant contribution from momenta of order k 1 ≃ k 2 ≃ k 3 ≃ mv dm ≪ 2πmα, and larger momenta are exponentially suppressed. However δ(∆ω) sets k 2 1 + k 2 2 = k 2 3 − α 2 m 2 /n 2 , forcing k 3 to be of order mα/n, therefore it is in the UV tail of f (k 3 ) (except for very large n, in which case, as described later, the matrix element is suppressed). In other words, as evident from Figure 8(a), the particle '3' must be emitted in the final state with (large) energy k 2 3 /2m ≃ mα 2 /2n 2 , to compensate for the energy of order −mα 2 /2n 2 that has been 'lost' to the bound state. As a result, the factor f (k 3 ) provides a suppression in C: direct capture is not enhanced by the Bose enhancement.
Similarly, in the process of stripping, energy conservation forces the incoming particle, '3' in Figure 8(c), to have (large) momentum of order mα/n, in order to strip a bound state while still providing an outgoing state with positive energy. As before, the factor f (k 3 ) suppresses Γ 2 in Eq. (46), i.e. particles in the DM halo are typically not energetic enough to strip out a bound state. On the other hand, for stimulated capture energy conservation requires the outgoing momentum, k 3 in Figure 8(b), to be of order mα/n. Since no density factor enters for k 3 , this term is not suppressed. This leads us to conclude that Γ > 0.
From the discussion above, the relevant momenta for all the integrals are either mα/n or mv dm ≪ mα. As a result, the velocities v at play are α/n or v dm and all satisfy 2πα/v ≫ 1. Thus, as explained in Section 3 and Figure 7 (left), we can approximate all scattering states in Eq. (46) with their Bessel limit of Eq. (21) (including those with 'large' velocity v ≃ α). Note that as 2πα/v dm approaches 1, for instance 2πα/v dm = 3 ÷ 4, C and Γ 2 are not exponentially suppressed anymore; however, δ(∆ω) still implies that k 1 , k 2 , k 3 satisfy 2π/v i ≳ 3 ÷ 4, and the Bessel approximation still holds for the integrands.
The Bessel form in Eq. (21) dramatically simplifies the matrix element in Eq. (31). For instance, for the 100 level this becomes where ζ ≡ x/R ⋆ is the coordinate in units of R ⋆ , and we omitted a phase that drops out in the absolute value. 18 This quantifies the spatial superposition of three scattering states ψ k (of the form in Figure 7 (right), oriented in three generic directions) with the ground state. In this superposition, the dependence on α has completely dropped out. This occurs because the typical region where ψ k is of order-one follows the ground-state radius R ⋆ (as α/v dm changes), as in Figure 7 (right); as a result, there is always overlap of order-one between the two. Additionally, as α/v dm gets larger, R ⋆ and therefore the volume of integration in M gets smaller, but at the same time the amplitude of ψ k increases. These two effects compensate to make M independent of α/v dm . Thus, contrary to the previous case, M is not suppressed. Note that the magnitudes of the momenta k i only appear as an overall factor. The remaining integral in Eq. (52) is, by rotational invariance, a calculable (order-one) dimensionless function F [k i ·k j ] of the three scalar products ofk i , clearly maximized when all momenta are parallel. 19 For the 100 level, similarly to Eq. (49), Γ = Γ 1 − Γ 2 becomes 8g 2 ρ 2 dm π 2 m 3 v 2 dm dp 1 dp 2 d cos θ 1 d cos where we redefined p i = k i /(mv dm ) and it is understood that p 3 = p 2 1 + p 2 2 + α 2 /v 2 dm and ϕ 3 = 0. The integral in Γ 1 is independent of α and v dm , while Γ 2 depends on α 2 /v 2 dm via p 3 . Thus, where a numerical calculation gives c ≃ 0.34, and s(a) arises from the stripping term: as anticipated, it is exponentially suppressed for a → ∞, and of order-one otherwise (see Eq. (113) in Appendix C). The result of Eq. (54) is that, when ξ foc ≫ 1, the timescale 1/Γ is of order τ rel regardless of the value of α. The external body acts as a 'catalyzer' for the formation of the gravitational atom, but the rate of accretion does not depend on its mass. We show this characteristic time, including the correction from s(a), in Figure 10 (left). This correction is of order one for 2πα/v dm ≃ 3 ÷ 4 (for smaller values the Bessel approximation is not valid anymore), but still leads to a positive Γ. The exponential increase starts from the mass M in ⋆ ≃ C/Γ accumulated by DM capture. However, given that capture is kinematically disallowed for ξ foc ≫ 1, this mass and the corresponding bound density are also suppressed in this limit. More precisely where h(a) ≃ 0.12 s(a)/(1−s(a)). We show the density ρ in (0)/ρ dm in Figure 10 (right). The function h(a) is exponentially suppressed for a → ∞. Consequently, if ξ foc ≲ 10, when the exponential growth starts the gravitational atom has a density comparable to the local DM density, as also shown in Figure 10. However, if ξ foc ≳ 10 the suppression of C is more significant and the initial density is much smaller. 20 Despite the suppression of M in ⋆ from direct capture, an overdense bound state will still form. In fact, the 100 state is irreducibly populated, possibly via the decay of higher levels (see below) or by quantum fluctuations. These will get exponentially enhanced and reach a macroscopic density -comparable to the local DM density -at the latest by t ≃ O(10)τ rel .
Let us finally briefly comment on the higher levels, focusing on n00. If the typical radius of the excited state (of order n 2 R ⋆ ) is smaller than λ dB = 2π/mv dm , which occurs for (small) n ≲ ξ 1/2 foc , the Bessel approximation of ψ k is accurate to evaluate M (see also footnote 18). In this case the matrix element for the n00 states is similar to Eq. (52), but with an additional overall factor n −3/2 , and e −ζ → e −ζ/n P n (ζ), where P n is a polynomial of ζ; see Eq. (70) in Appendix A. Depending on the form of P n , M could be of the same order as for the ground state.
In fact, from a direct calculation, the matrix element M for e.g. nlm = 200 is of the same order as that of 100 level in the limit ξ foc ≫ 1, so this level is populated at a similar rate. By studying a system with the 100 and 200 states only, we show in Appendix C that for ξ foc ≳ 1 the higher level in any case quickly decays (in a timescale comparable or shorter than τ rel ) into the ground state, via the nonlinear terms omitted from Eq. (45). Note that direct capture is less exponentially suppressed for these higher n00 states (because α/v dm in Eq. (55) is substituted with α/(nv dm ) with n > 1); therefore, the ground state will be also populated by relaxation of higher levels. Consequently, the result of this Section should be thought as a lower bound to the captured mass.
On the other hand, for (large) n ≳ ξ 1/2 foc , the effective integration region in Eq. (31), set by the radius of the excited state, is larger than λ dB and one should use the full scattering states in Eq. (18). These resemble plane waves for r ≳ λ dB and, similarly to Eq. (48), this behavior suppresses the integral in Eq. (31). Consequently, for states with n ≫ 1, M is expected to be small compared to the first few levels, and these are not efficiently populated to start with.
Intermediate regime: ξ foc ≃ O (1) Due to the complex hypergeometric function in Eq. (18), it is not feasible to calculate C and Γ analytically for intermediate values of ξ foc ≃ O(1); in Figure 10 we show a naive interpolation of this region by the dotted lines. Nevertheless, we expect that Γ changes sign at some critical value ξ foc = O(1), such that for larger ξ foc the density increases exponentially, and saturates otherwise. In fact, by simulating the system, in Section 5 we provide evidence that this is the case, and we estimate the critical value to be ξ foc ≃ 1. The numerical simulations also suggest that |Γ| −1 remains of order τ rel even for ξ foc ≃ 1, and as a result the interpolation in Figure 10 likely overestimates |Γ| −1 in the transition region.

The fate of the gravitational atom and instability
As anticipated in Section 2, the exponential growth for ξ foc ≳ 1 stops once the density grows beyond In the second equality of Eq. (56) we fixed r = R ⋆ . Although the sign of g does not play a role in the formation of the gravitational atom, when ρ ≃ ρ crit this does in fact affect the bound state evolution and stability. Note that in this regime the self-interactions can no longer be treated perturbatively.
• Once the density reaches ρ crit , attractive self-interactions (g < 0) cannot be compensated by the gradient pressure and destabilize the bound state, leading to its collapse in a manner analogous that of a self-gravitating boson star, widely studied in [59,60,[107][108][109]. During this collapse, the density increases until the point where higher order terms in the potential of Eq. (12) become relevant in the dynamics, and a full relativistic treatment is needed. In the related case of self-gravitating boson stars for an axion-like potential, in the final moments of the collapse the typical field value is of order f a and the self-interactions lead a rapid emission of relativistic scalar particles (through e.g. 3 → 1 processes [110,111]), known as a Bosenova explosion [59,60]. 21 The analogous process for a collapsing gravitational atom is expected to release an order-one fraction of the bound-state mass into relativistic scalar radiation, which would present novel opportunities for detection. The collapse of a solar halo is worthy of dedicated study, though we discuss some implications in Section 6.
• If instead the self-interactions are repulsive (g > 0), as soon as the density reaches the critical value, the balance of forces supporting the gravitational atom changes, such that the external gravitational potential is balanced by the self-interactions rather than the gradient pressure. The critical mass (corresponding the density ρ crit ) at which this happens is By solving mΦ = gρ/m at r = R (g>0) ⋆ , using ρ(R ⋆ ) ≃ M ⋆ /(πR 3 ⋆ ), we find that this new bound state has radius The radius R (g>0) ⋆ is larger than the Bohr radius (mα) −1 and grows as √ M ⋆ . 22 At this point, our analytic analysis of the capture rate no longer directly holds; numerical simulations (see next Section and Appendix D) suggest that the density quickly saturates after this time. A detailed analysis is worthy of further study.
As mentioned in Section 2, for the values of f a and m in Eq. (56) that predict τ rel of the order of the age of the Solar System, the critical density is orders of magnitude larger than ρ dm . 23

Comparison with numerical simulations
To check the analytic expectations of Section 4, we evolve the EoM in Eq. (15) on a discrete lattice. The simulation requires one to resolve the UV scale of the α/r divergent potential, which is the radius R s of the Sun, 24 so the EoM are solved with α/r → α(3 − (r/R s ) 2 )/2R s for r < R s . The spacing between grid points should be small enough to resolve both R s and the radius R ⋆ , and the box needs to contain enough de Broglie wavelengths to resemble the effectively-infinite volume of our galaxy.
We set the initial conditions as a random realization of the wave superposition ψ w in Eq. (23) with f (k) fixed to its Gaussian form in Eq. (26) and dispersion σ = v dm / √ 2. As shown in Appendix D, by rescaling coordinates and time as x → mv dm x and t → mv 2 dm t, the EoM and the initial conditions depend only on the dimensionless combinations ξ foc = 2πα/v dm ,g ≡ gρ dm /(m 2 v 2 dm ), andR s ≡ mv dm R s ; in particular, the strength of the self-interactions g always enters together with ρ dm . Althoughg ≃ O(10 −8 ) for the values of m and g for which dense atoms form around the Sun in 5 Gyr, in the following we will fix |g| = 0.006 andR s = 0.3. These values do not affect the interpretation of the simulation results, but allow (i) the relaxation time τ rel to be short enough (and therefore the atom to form within the available simulation time), and (ii) the Sun to be easily resolved by the grid. More details on the simulations can be found in Appendix D. Figure 11 (left) shows the evolution of the overdensity ρ(0)/ρ dm at center of the Sun r = 0 for a fixed initial condition and different values of ξ foc = 2πα/v dm . All the lines have g < 0 except for ξ foc = 2.5, 1.5 for which g > 0 (see below for discussion). At early times, the density fluctuates over the coherence time of the DM waves, 2π/mv 2 dm ; therefore, we also plot the average over times longer than this period (solid thick lines). Although the fluctuations depend on the stochastic phase of the waves in the initial conditions, their average should be independent of this phase. To compare with our expectations from Section 4, we show time in units of the relaxation time τ rel in Eq. (47). We also indicate with a solid disk the prediction of the exponential timescale 1/Γ in Eq. (47) based on the Bessel approximation, valid when ξ foc ≳ 3 ÷ 4 (see Section 4 and Figure 7), and with an empty disk (placed at the end of each simulated range) the value of the critical density ρ crit /ρ dm = ξ 2 foc /(2π 2g ) in Eq. (56), multiplied by an order-one factor equal to 0.7.
For values of ξ foc definitely larger than 1, the density increases at the predicted timescale. In fact, the growth might be faster than exponential, because also the first few excited levels also experience exponential growth, and then quickly decay to the 100 level, as discussed in Section 4.3 and Appendix C. 22 With an abuse of notation, we wrote g = 1/(8f 2 a ) in Eq. (57) for repulsive self-interactions. Note that Eq. (58) is analogous to the Thomas-Fermi radius [116] of repulsively-interacting boson stars, and the mass in Eq. (57) is analogous to the critical mass of a self-gravitating boson star [36][37][38]; both are easily obtained in the limit M → M⋆. In the repulsivelyinteracting case, there is an instability stimulated by general-relativistic effects, though the critical mass is too large to be relevant here; it is also easily obtained using the M → M⋆ limit of the self-gravitating case [35]. 23 Note that we can safely ignore the instability in the regime ξ foc ≪ 1, because the saturated density is much smaller than the critical density, as can be understood by comparing Eqs. (9) and (56) for the appropriate values of m and fa. 24 Formally we distinguish the physical radius of the Sun, R⊙ ≃ 7 · 10 5 km, from the variable Rs used in the simulation. Note that the density continues to oscillate also during the exponential increase. These oscillations are not captured by our analytic derivation, which only determinesṄ ⋆ averaged over times much larger than 2π/mv 2 dm ; see Eq. (41). 25 For g < 0, once the density is large enough, the halo quickly collapses and the nonrelativistic energy d 3 x ϵ (see Eq. (16)), is not conserved anymore in the simulation regardless of the time and space-steps, and the simulation is stopped; to describe the collapse and subsequent evolution of the bound state, the full relativistic EoM are needed. In particular, the density when this occurs in the simulation beautifully matches our predicted critical density, as shown by the empty disks. On the other hand, for g > 0, energy conservation continues to hold, as the exponential increase stops and the density saturates shortly after the critical density is reached; see also Figure 21 in Appendix D.
As expected, as ξ foc decreases, the exponential growth timescale increases, but the latter remains of the order of τ rel even for ξ foc ≃ 1, where the Bessel approximation of Eq. (21) has fully broken down. For values of ξ foc close to 1, given the large values ofg that are feasible for the simulation, the critical density is so small (just a few times the background density) that energy non-conservation occurs soon after the increase. Thus, for these values of ξ foc it is more feasible to test the exponential increase for g > 0. Such simulations, shown in Figure 11 (left) for ξ foc = 2.5, 1.5, provide evidence that the exponential increase (i.e. Γ > 0) holds for ξ foc close to 1, or even as small as ξ foc ≃ 1 as suggested by Figure 21 in Appendix D. The exponential timescale continues to be of order the relaxation time, as expected, although we do not For the sake of illustration, we assume that the contraction and Bosenova explosion -expected to occur once the critical density is reached -are immediate, and that the process then repeats from the beginning. The shaded region represents the direct constraint on the density from Solar System ephemerides [61]. attempt to extract it precisely in this work.
In Figure 11 (right) we also show the 'spherically averaged' density profile, obtained by averaging ρ(x) over a spherical volume shell of radius r centered at r = 0 for ξ foc = 1.5, 4, 8 during the exponential increase. Also the profile oscillates in time, so we show its time-average over times larger than 2π/mv 2 dm . The profile matches that of the ground state with the expected radius given in Eq. (2) (dashed lines). This holds throughout the exponential increase, as shown in Figure 20 of Appendix D, where the profile is shown at two subsequent times. These results nicely confirm our analytic analysis, in particular that for ξ foc ≳ 1 the mass of the atom increases exponentially over a timescale comparable to the relaxation time.

Solar gravitational atoms
As anticipated in Section 2, the exponential growth of the gravitational atom around the Sun occurs for (1 ÷ 2) · 10 −14 eV ≲ m ≲ 2 · 10 −13 eV. The upper value corresponds to a radius R ⋆ smaller than the radius of the Sun, for which the point-like approximation of the Sun is not applicable; see the intersection between the blue and green lines in Figure 3. In the following we will focus on the mass range above.
As an illustrative example, we sketch in Figure 12 the expected time evolution of the overdensity at the position of the Earth, ρ(AU)/ρ dm , resulting from the formation of the gravitational atom bound to the Sun, for m = 10 −14 eV (left) and m = 2 · 10 −14 eV (right), and λ = −m 2 /f 2 a < 0. In both cases we show the result for two values of f a . Figure 12 assumes that the exponential time scale is of order τ rel also for m as low as m ≃ 10 −14 eV (corresponding to ξ foc slightly smaller than 1) and that the collapse and Bosenova explosion are immediate.
As can be seen in Figure 12, the self-interaction strength determines the capture timescale via Eq. (47) and the maximum density via Eq. (56). For m ≃ 10 −14 eV, self-interactions with f a ≃ 10 7 ÷ 10 8 GeV are strong enough for the exponential time-scale to be of order the age of the Solar System, while predicting a density ρ(AU) > ρ dm at the position of Earth. Additionally, larger m requires a smaller f a to achieve the same density. As mentioned in Section 2, the DM overdensity in our Solar System is constrained by direct observation using Solar System ephemerides, at the level of ρ/ρ dm ≲ 10 5 (10 7 ) at a distance of r ≃ 1 (0.3) AU from the Sun [61]. This is shown by the gray regions in Figure 12. 26 Note that for a fixed value of m, there is an 'optimal' value of f a , which we denote asf a , for which the density after t = 5 Gyr is maximized (this is obtained by requiring ρ(t = 5 Gyr) = ρ crit ). For f a ≳f a , the growth timescale increases as in Eq. (47), while for f a ≲f a , the critical density in Eq. (56) decreases; in both cases, the final density is smaller than that for f a =f a . The result for the maximum overdensity as a function of m is shown in Figure 13, where we also report in the upper axis the corresponding optimal valuef a . We illustrate the density at a distance of r = {1, 0.39, 0.1} AU (blue, yellow and green lines); the latter two may be relevant for proposed DM searches aboard space missions [117]. The dashing between m = (1 ÷ 2) · 10 −14 eV indicates the uncertainty in the point where exponential growth becomes relevant, near ξ foc ≃ 1. We use the DM parameters of the standard halo model v dm = σ/ √ 2 = 240 km/s (left panel) and, as an example, also v dm = σ/ √ 2 = 50 km/s (right panel), which could be plausible is the DM resides in a low-dispersion dark disk (see e.g. Refs. [53,63,64]).

Solar halos with λ > 0
For λ > 0, as the density reaches the critical value in Eq. (56), the (repulsive) self-interaction balance the (attractive) gravitational potential of the Sun, and the radius is given in Eq. (58). While no Bosenova occurs in this case, our prediction of a large, stable overdensity is striking and worthy of further study.

Effect of Standard Model couplings
We have so far ignored possible direct couplings of ϕ to the SM. These, however, determine the DM direct detection prospects. In this Section we provide evidence that such couplings do not significantly alter the formation of the gravitational atom in the Solar System. We will consider the case where ϕ is a pseudo-scalar field, so that these couplings take the form where ψ is a SM fermion field and F µν is the photon field strength tensor, withF µν = ϵ µναβ F αβ /2. The coupling constants g ϕγγ and g ψ in Eq. (59) are model-dependent and not necessarily related to the self-interactions, although in many plausible axion theories they are fixed in terms of f a . The interaction with the photon in Eq. (59) induces the decay ϕ → γγ, which could deplete the gravitational atom. 27 The corresponding rate in vacuum, is far too slow to be relevant on Solar-System time scales (the benchmark value g ϕγγ ≃ 10 −11 GeV −1 in Eq. (60) is close to the upper limit for this coupling, see e.g. [78][79][80]). On the other hand, 'stimulated' decay of the scalars can be induced in the background of a radiation field, e.g. by the large number of photons around the Sun. This enhances the decay rate by a factor of the occupation number of (soft) photons with frequency ω ≃ m/2; see e.g. [118]. Since the temperature of the Sun and therefore the typical photon energy are ≫ eV, this process is not enhanced in our case.
In the background of a large scalar density ρ, also parametric resonance enhances the decay ϕ → γγ. This effect has been analyzed for both homogeneous and inhomogeneous field configurations and applied to boson stars [119][120][121], and is negligible if where R ≃ R ⋆ = 1/(mα) is the typical size of the region where the density is ρ. This condition is satisfied for the values of m and ρ for the solar halo, as shown for instance in Figure 4. 28 SM couplings could also disrupt the bound state by inducing 2 → 2 scattering processes that eject bound axions out of the gravitational atom. One example is X + 100 → X + k, with X either the photon or the electron. Because the typical energy of photons and electrons from the Sun (of order keV) is larger than the binding energy of the atom, the amplitude for this process is suppressed compared to that of the capture processes in Figure 8. Additionally, unless the emitted axion remains nonrelativistic, the final state will see negligible Bose enhancement, resulting in further suppression of the scattering rate from SM particles.
A second possibility is that the scalars in the bound state get absorbed as in a process of the form e + 100 → e + γ (where e is the electron field), studied for example in [123]. We merely note that absorption is the inverse of a Compton-like production process, e + γ → a + e, which was studied in e.g. [124,125]; its effect was found to be negligible for m ≪ keV due to a very large absorption timescale τ abs ∝ (e m/T − 1) −1 , where T is the temperature in the Sun. In our case, there will be further suppression due to the fact that the overlap region between the gravitational atom and the Sun is a small fraction of the total volume of the bound state (as long as R ⋆ ≫ R ⊙ , or equivalently m ≲ 10 −13 eV). We therefore expect axion absorption to be too slow to significantly impact the formation of the gravitational atom.
The discussion above suggests that, for the allowed values of the SM couplings in Eq. (59), the formation of the solar halo is unlikely to be affected by solar dynamics of e.g. photons, electrons, or other SM particles.

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 δ ≡ (ρ − ρ)/ρ follows [126][127][128][129][130][131] where ρ ∝ a −3 is the average DM energy density, k is the comoving momentum, H ≡ȧ/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, ± refers to attractive/repulsive self-interactions). The two comoving momentum scales k λ (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 and k J ≡ a 16πGρm 2 1 4 ≃ 10 Mpc −1 a a eq where a eq is the scale factor at the time of matter-radiation equality. Note that Eq. (62) holds only when the mode δ 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 δ k ∝ 1 + (3/2)a/a eq ) for k larger than the critical momenta k λ 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 δ k (until δ k ≃ 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 δ 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-α data), i.e. approximately k ≃ 1 Mpc −1 (see e.g. [132][133][134]). Given that k λ ∝ a and k J ∝ a 1/4 , see Eqs. (63)(64), the corresponding bound on m and λ 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 eq λ ≡ k λ | a=aeq ≳ 1 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 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.
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 λ. 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 δ k once the modes enter the horizon, so that δ k increases logarithmically from the time when the mode enters the horizon until matter-radiation equality. In particular, the mode at k = 1 Mpc −1 enters the horizon at redshift z i ≃ 5 · 10 5 , corresponding to the temperature T i ≃ 100 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 where y ≡ a/a eq , δ ′ ≡ dδ/dy and we neglected the quantum pressure term, which is an excellent approximation for m ≳ 10 −21 eV. The solution of Eq. (65) deep into radiation domination starting from a = a i (when the perturbation enters the horizon) with initial conditions δ k (a i ) = 1 andδ k (a i ) = 0 is Eq. (66) shows that, soon after the time when a = a i , perturbations with momentum k ≳ k eq λ a i /a eq increase substantially during radiation domination, receiving an enhancement factor much larger than one, in particular of order (a eq /a i ) 1/4 exp [ 6a eq /a i k/k eq λ ] by the time of matter-radiation equality. 29 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 eq λ a i /a eq ≳ 1 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/ a i /a eq ≃ 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 ≳ k eq λ a i /a 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 |Γ| −1 is 5 Gyr and 150 Myr respectively: For m ≲ 10 −14 eV, Γ 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 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 ÷ 10 8 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/ρ < 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 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].

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 29 The expression in Eq. (66) is a solution only in the limit k ≳ k eq λ a/aeq, 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 ≳ k eq λ ai/aeq. 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 ξ foc ≡ λ dB /R ⋆ = 2πα/v dm between their de Broglie wavelength, λ dB , and the gravitational Bohr radius, R ⋆ in Eq. (2). In particular, the bound state either saturates to an (underdense) steady-state configuration (when ξ foc ≲ 1, i.e. m ≲ 10 −14 eV for the Sun) or it exponentially grows in density (when ξ foc ≳ 1). We have confirmed this behavior through numerical simulations. On the other hand, the size of the self-interaction coupling λ only dictates the formation time scale, which in the regime ξ foc ≳ 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 ρ dm and v dm ), and the magnitude of λ. Note that the basic mechanism depends neither on the sign of λ nor the parity of the field.
The most striking implications occur in the exponential-growth regime, when ξ foc ≳ 1. For instance, for axion DM with mass m ≃ few · 10 −14 eV and decay constants f a ≃ 10 7 ÷ 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 AU from the Sun; see Figures 4 and 6. The density at r ≪ AU can be even larger for smaller masses, m ≃ 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 (λ > 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 (λ < 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 ⋆ is larger than the radius R of the astrophysical body that captures the ULDM. For large enough axion masses (m ≳ 2 · 10 −13 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 ⋆ ). Here we merely note that when ξ foc ≳ 1, the Bohr radius R ⋆ of the atom bound to a few different astrophysical systems is given by: where for reference R J = 7 · 10 4 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 ≲ 10 −9 eV [44]. However, for these masses and M = M ⊕ (the mass of the Earth), ξ foc ≲ 0.1 and therefore the exponential growth does not occur. Using Eq. (9) one finds that ρ eq ≪ ρ 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π/v dm ) than those for which DM capture starts to be effective, ξ foc ≃ 1. 30 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 ∼ O( 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 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 ν ≡ v dm /v b (ν is analogous to 2π/ξ foc ): Boson stars with ν ≫ 1 decrease in mass as a result of processes analogous to stripping, see Figure 8(c). In the other limit, ν ≪ 1, their mass increases, but at a rate that decreases in time, as M 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. 30 The analysis of solar-mass and supermassive spinning black holes led to constraints in the mass range m ≃ 10 −13 ÷ 10 −11 eV [138,140] and m ≃ 10 −21 ÷ 10 −17 eV [141], respectively.

A Energy density and bound states
In this Appendix we derive the energy density ϵ in Eq. (16) and provide the expression of excited states ψ nlm .
Energy density. As mentioned in Section 4.2, in the nonrelativistic limit the action of the scalar field reduces to S ≡ dtd 3 xL in Eq. (42). If the self-gravitational potential Φ se (neglected in the main text) is relevant, L contains the additional terms −(8πG) −1 |∇Φ se | 2 − mΦ se |ψ| 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 (with Φ = Φ se + Φ ex ) is conserved on the EoM, when integrated over the whole space. Eq. (68) can be simplified using |∇Φ se | 2 = ∇ · (Φ se ∇Φ se ) − Φ se ∇ 2 Φ 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 is globally conserved.
Bound states. The quasi-stationary solutions ψ = e −iωt Ψ(x) of the EoM in Eq. (15) satisfy the eigenvalue equation (−∇ 2 /2m − α/r)Ψ = ωΨ, whose most general solutions with negative energy (ω < 0) are the well-known hydrogen-like bound states with eigenvalues ω = ω n = −mα 2 /(2n 2 ). In Eq. (70), Y m l (θ, φ) are the spherical harmonics and L b a [x] the generalized Laguerre polynomials of degree a and order b. The bound states ψ nlm and the scattering states ψ k in Eq. (18) are orthonormalized, and are all orthogonal to one another, as

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 ⟨Ṅ nlm ⟩ in Eq. (41). This is obtained by solving the equation of motion with initial condition ψ(t = 0) = (ψ w + ψ b ) t=0 in Eqs. (23) and (35), perturbatively in the (local) effective self-interaction parameter g eff = gρ/(mω), where ρ and ω are the typical values of the local density and frequency of the field. As mentioned in Section 4.1, in the regime ρ < ρ 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 λ = g = 0, i.e.
where ψ nlm is given in Eq. (70) and the coefficients c nlm and c k depend on λ (or g) and vary in time more slowly than the typical field oscillation periods, given by 1/ω nlm and 1/ω k for bound and scattering states respectively. We can expand c nlm and c k in a series in g eff as where the superscripts indicate the order in g eff of each term. In generalω nlm andω k in Eq. (73) are different from the original eigenfrequencies ω n and ω k , which indeed receive time-independent corrections (or 'shifts') when λ ̸ = 0; see for instance [138]. These corrections may be also expanded as where ω n = −mα 2 /2n 2 and ω k = k 2 /2m are the original λ = 0 frequencies, and |ω (1) nlm | ≪ |ω n |, etc.. As a result, we can organize the full solution ψ in Eq. (73) in the perturbative series where, for instance, 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 and implies that c k (t) are in fact time-independent. Imposing that ψ equals the initial condition ψ w + ψ b at t = 0, we find ψ (0) = ψ w + ψ b , and c nlm Re[c (2) nlm (t)] + . . . . (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 ⟨Ṅ nlm ⟩.

First-order solution
The first order EoM is which is a source equation for ψ (1) . We solve Eq. (80) with initial conditions containing a generic configuration of scattering states a(k), but only particles in the ground state, i.e. N (0) nlm = 0 unless nlm = 100, or equivalently With these initial conditions, ⟨Ṅ 100 ⟩ will describe the evolution of the bound-state occupation number both before the ground state is populated (upon taking N Using the orthonormality relations of Eq. (71) we can single out the contribution of the individual boundstate levels and waves: Substituting the expression of ψ (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(x)e −iωst ⊂ g|ψ (0) | 2 ψ (0) that oscillates with frequency ω s contributes to a significant change in c nlm or c 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 k have a chance of growing with t instead of oscillating. In particular (see also [138]): 100 , as all other diagrams do not conserve energy or do not change the particle number in the 100 state.
1. In Eq. (83), terms with ω s < 0 and ω s = ω 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 ω (1) nlm . If instead ω s ̸ = ω n for all n, the term is a driven off-resonance oscillator, and represents a (suppressed) off-resonant production of particles.
2. In Eq. (84), terms with ω s > 0 correspond to a resonant production of scattering states or a correction to their frequency ω (1) k . As we will see, only resonant terms will effectively contribute to c nlm and c 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 k correspond to the factors of ψ 100 or ψ 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 ψ * 100 (or in general ψ * 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: • Diagram (4) in Figure 15 corresponds to g N (0) 100 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(ω 1 −ω 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 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 * (k 2 )a(k 1 ) as ⟨a * (k 2 )a(k 1 )⟩ + [a * (k 2 )a(k 1 ) − ⟨a * (k 2 )a(k 1 )⟩]. Using Eq. (24), the first part becomes 2g N (0) 100 and thus represents a shift of the energy of the ground state due to the dark matter waves, ω (For instance, in the plane wave limit of ψ k , the frequency correction is ω (1) 100 ⊃ 2gn dm , with n dm = ρ dm /m being the number density of the waves. Thus the perturbative condition |ω (1) 100 | < |ω 1 | is equivalent ρ dm < ρ crit /8.) The second part is proportional to a * (k 2 )a(k 1 ) − ⟨a * (k 2 )a(k 1 )⟩, 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 ⟨Ṅ 100 ⟩ at this order. Although it actually affectsċ 100 , it can be verified that its contribution to N 100 (t) at second order in g eff is canceled between ⟨|c (1) 100 | 2 ⟩ and 2 N (0) 100 Re[⟨c (2) 100 ⟩] 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 → k 3 + 100, and is the only diagram in Figure 15 that leads to an effective change in c where This reproduces Eq. (30) in the main text, and, correspondingly, diagram (5) reproduces c (1) nlm in Eq. (34). 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 ̸ = 100 and c (1) k , whose relevant diagrams are shown in Figure 16 (left and right, respectively). The result for c (1) nlm is 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 |ω n | < |ω 1 |. The result for c 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 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ñlm, i.e.
in which case it is sufficient to substitute 100 →ñlm 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

Occupation number and second-order solution
As mentioned in Section 4.1, the contribution at order g eff to ⟨N 100 ⟩ (the second term in Eq. (79)) vanishes, since ⟨c (1) 100 ⟩ ∝ ⟨a(k 1 )a(k 2 )a * (k 3 )⟩ = 0, and the first nontrivial correction to ⟨N 100 ⟩ comes at order g 2 eff and reads ⟨|c (1) 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 ⟨a * (k 1 )a(k 2 )⟩ in Eq. (24), along with the identity lim with ∆ω ≡ ω k 1 + ω k 2 − ω k 3 − ω 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: Similarly to Eq. (80), this is a source equation for ψ (2) , generated by a combination of ψ (0) and ψ (1) . By plugging the explicit expression of ψ (2) in terms of c (2) nlm and c (2) k and using again the orthonormality relations, we obtain (as in Eq. (85)) iċ (2) where the dots stand for frequency corrections analogous to those in Eq. (85). It is via ψ (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 ⟨Ṅ 100 ⟩.
We can substitute into Eq. (92) the expressions of ψ (0) in Eq. (81) and the solution of ψ (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 ψ (0) , i.e. [dk]a(k)e −iω k t ψ k ⊂ ψ (0) , gives rise to processes that modify the particle number.
• Additionally, only the unbound component of ψ (1) , i.e. [dk]c 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 ψ (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 k 1 , k 2 , 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 ⟨Ṅ 100 ⟩, see below.
Using the identity which follows from Eq. (24) and Wick's theorem, Eq. (93) simplifies to The time integration can be performed explicitly using t 0 dt ′ t ′ 0 dt ′′ e iω(t ′ −t ′′ ) + c.c = 4ω −2 sin 2 (ωt) → 2πδ(ω)t, and yields ⟨c (2) 100 + c 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 (2) 100 ; this term gives Plugging the last two equations into Eq. (89) we obtain the full result: 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 [dk] is substituted with nlm (this is because the diagrams (2, left) and (2, right) in Figure 16 are identical upon substituting k ↔ 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 (0) 100 ). 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 ξ foc = 2πα/v dm ≫ 1. Consequently, N 100 increases. As a result, the excitation terms are not expected to affect the discussion in the main text.

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 M ≡ M k 1 +k 2 →k 3 +nlm in Eq. (31) is suppressed for n ≫ 1, in both limits of ξ foc ≫ 1 and ξ foc ≪ 1, which inhibits the population of very high-n states.
However, we also argued that for ξ foc = 2πα/v dm ≫ 1 the magnitude of M for the first few levels (those satisfying n ≲ ξ 1/2 foc ) could be comparable to that of nlm = 100 (this can be seen, for instance for nlm = 200, by using the explicit expression of ψ 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 ξ foc ≪ 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 ξ foc ≫ 1, the excited state is rapidly populated but transitions to the ground state over a timescale similar to the exponential increase time, τ 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: 200 = 0 at t = 0 in our Solar System, but, as in Section 4.1, keeping them generic will allow us to study the rates ⟨Ṅ 100 ⟩ and ⟨Ṅ 200 ⟩ at an arbitrary time after such levels get populated from the DM capture.
As in Appendix B, the occupation number ⟨N 100 ⟩ is given in Eq. (89), and requires the calculation of c The ⟨|c (1) 100 | 2 ⟩ part of ⟨N 100 ⟩ is simply given by the square of the amplitude of each term in Eq. (99) (the cross terms vanish upon average) and reads ⟨|c (1) 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 → 100 + k 2 and 200 + 200 → 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  As in Appendix B, one could proceed to calculate the second order coefficient c (2) 100 , and thus the second part of ⟨N 100 ⟩. The only difference is that in ψ (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 (0) 200 . Rather than going through this lengthy derivation, we can obtain (heuristically) the form of the terms in ⟨Ṅ 100 ⟩ that depend on N (0) 200 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 ↔ 100 + k 2 and 200 + 200 ↔ k + 100 (which include direct and inverse processes) is consistent with Bose enhancement. In other words, the contribution to ⟨Ṅ 100 ⟩ from the former should be proportional to while for the latter to where the approximations are valid for N nlm ≫ 1 (for simplicity of notation, here and in the following we omit the superscript '(0)'). Then, the evolution of the occupation number ⟨Ṅ 100 ⟩ 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 iṡ One can similarly obtain the evolution of ⟨N 200 ⟩, which readṡ and is analogous to Eq. (103), provided that N 100 and N 200 are exchanged, appropriate minus signs added and a factor of 2 in the last line of Eq. (104). Note that: • The first two lines of Eq. (103) 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. (97). 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. (103-104) 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 (0) 300 ̸ = 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. (103)(104)) describing the simultaneous evolution of the occupation number of all levels. However, in the following we will prove that: • In the limit 2πα/v dm ≪ 1, the two-level system in Eqs. (103)(104) 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πα/v dm ≫ 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 lowerlevel 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: ξ foc = 2πα/v dm ≪ 1 We first consider the limit 2πα/v dm ≪ 1. As mentioned in Section 4.3, M k 1 +k 2 →k 3 +n00 ∝ n −3/2 ; thus, the capture and stripping rates for the 200 state are approximately 1/8 times those for 100. Energy conservation, ω k 1 + ω 1 − ω 2 − ω k 2 = 0, forces k 1 ≃ k 2 in the 2πα/v dm ≪ 1 limit. As a result, since the sum of the excitation and de-excitation terms is proportional to f (k 2 ) − f (k 1 ), these cancel between each other. Consequently, Eqs. (103) and (104) As in Section 4.3, the coefficients a 1 and a 2 are with A ≃ 0.21 and B ≃ 0.08. The other coefficients may be evaluated in a similar way: where f 4 (a) = 2 √ 2a/e for a ≪ 1. 31 These equations are greatly simplified when rewritten in terms of t ′ = t/τ rel and of the normalized occupation number Note thatÑ n00 represents (modulo a πn 3 factor) the contribution ρ n00 from the state n00 to the density of the bound state at the center r = 0, normalized to ρ dm (i.e. ρ(0)/ρ dm ). In particular, for the singlelevel system in Eq. (41) For ξ foc = 2πα/v dm ≪ 1 the equilibrium solution (i.e.Ṅ nlm = 0) of this system occurs at N 100 = N 200 ≡ 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Ñ 100 −Ñ 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Ñ 2 ∝ (α/v dm ) 7 and (v dm /α) 2Ñ 3 ∝ (α/v 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Ñ ≪ 1. Figure 19 (left) shows the numerical evolution of the system in Eqs. (109)(110) starting from N 100 = N 200 = 0 for ξ 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.
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  Figure 19: The time evolution of the normalized occupation numbersÑ 100 (blue) andÑ 200 (orange), by numerically solving the two-level system, Eqs. (109)(110), with ξ foc = 2πα/v dm = 0.1 (left), and Eqs. (111)(112) with ξ foc = 10 (right).Ñ n00 /(πn 3 ) represents the contribution from the state n00 to the central density, normalized to ρ 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 ≃ τ rel .
powers of N 100 and N 200 . Ignoring all the terms that are exponentially suppressed except for direct capture, in the 2πα/v dm ≫ 1 limit Eqs. (103)(104) are well approximated by: where the coefficients can be calculated similarly to Section 4.3, 32 and read where to a good approximation s(x) ≃ e −x+ √ 4+2x−2 , and D ≃ 0.035g 2 α 2 mρ dm v dm , R ≃ 1.3 · 10 −4 g 2 α 4 m 5 .
The solution of Eqs. (111)(112) with initial conditions N 100 = N 200 = 0 can be understood qualitatively. First, notice that Γ 100 and Γ 200 are of the same order (because the matrix elements M k 1 +k 2 →k 3 +n00 with n = 1, 2 are similar), so we expect also N 200 to increase exponentially, with similar timescale τ rel . The capture coefficients C 100 and C 200 are comparable, but since s(x) is exponentially decreasing function of α/(nv dm ) with n = 1, 2, the latter is less suppressed at large 2πα/v dm . Ultimately, this is because the exponential suppression from f (k 3 ), with k 3 = k 2 1 + k 2 2 + α 2 m 2 /n 2 , is smaller for n = 2 than n = 1. Thus, as mentioned in Section 4.3, for 2πα/v dm ≫ 1 the mass captured until the time 1/Γ 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 ≃ τ 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 . 33 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 ∝ N 200 ∝ e t/τ rel , we can estimate that the latter start dominating over the former, inducing the decay of the 200 level, at times t c,D ≃ τ rel log D and t c,R ≃ 2τ 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 ≃ log(15v 2 dm /f 2 α 2 ) and log R ≃ 2 log(27v 2 dm /f 2 α 2 ); here f 2 is the density of the 200 state in units of the DM density at time t ≃ τ rel . We conclude that, soon after the exponential increase starts, 200 decays into 100. Figure 19 shows the time evolution ofÑ 100 andÑ 200 for 2πα/v 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 ≳ 2τ 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.

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 {α/v dm ,g ≡ gρ dm /(m 2 v 2 dm ),R s ≡ mv dm R s }. This can be easily seen by rewriting (i∂ t +∇ 2 /2m−mΦ)ψ = g|ψ| 2 ψ in terms of the rescaled variablesx ≡ mv dm x,t ≡ mv 2 dm t andψ ≡ ψ/ ρ dm /m (note that space and time are normalized to the typical wavelength and period of the DM waves, and mψ 2 to their average density; the relaxation time is simply mv 2 dm τ rel = 1/g 2 ). In this way, the EoM become i∂t +∇ whereΦ = − α/v dm r forr >R s (i.e. away from the Sun) andΦ = − α/v dm 2Rs (3 −r 2 /R 2 s ) forr <R s (inside the Sun). These redefinitions are useful because the initial conditions only depend on α/v dm , and read where we took σ = v dm / √ 2 and k dm parallel toẑ, and it is understood that R ⋆ → v dm α in ψk(x). In practice, a realization of these initial conditions is (at t = 0) Evolution and systematics. We solve Eq. (115) with initial conditions in Eq. (117) in a cubic box of length L with N 3 x = 64 3 ÷ 128 3 points. 34 The coefficient A in Eq. (117) is in this case A = 8(πL) 3/2 , whereL ≡ mv 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,58,146], to which we refer for all the details. (Note that the two additional parameters {N x ,L} need to be specified in a simulation.) The EoM imply the conservation of the total mass, m d 3 x|ψ| 2 , (which is automatic with the pseudospectral method we use) and of the nonrelativistic energy E ≡ d 3 xϵ, where ϵ is the energy density in Eq. (16). The time-independence of E is therefore a nontrivial check of the numerical evolution. We compute E throughout all the simulations, check its conservation and rely on simulations where this is conserved better than δE/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.
(a) The lattice spacing ∆ ≡ L/N x should be small enough to resolve the wavelength of the DM waves and the Sun radius. We checked that this requires ∆mv dm =L/N x ≲ 1 and ∆/R s =L/(N xRs ) ≲ 1 respectively, i.e. there should be at least approximately one grid point per inverse DM momentum and per Sun radius.
(b) The length L should be large enough for the box to be in the infinite volume limit, which requires L/(2π/mv dm ) =L/2π ≫ 1. In other words, the number of wavelengths per box side should be large;L/2π ≳ 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.
(c) Finally, R s should be much smaller than R ⋆ for the point-like approximation of the Sun to hold, i.e. R s /R ⋆ = (α/v dm )R s ≪ 1.
In the simulations of the main text ( Figure 11 in Section 5), we used |g| = 0.006 (g > 0 for the two lowest values of 2πα/v dm andg < 0 otherwise) andR s = 0.8 (except for the two highest values of 2πα/v dm , for whichR s = 0.4). We choseL = 40 and performed simulations with both N 3 x = 64 3 and N 3 x = 128 3 to check the lattice spacing requirement (a). Finally, the function γ(k) in the initial conditions in Eq. (117) has been fixed to be the same (random) function ofk for all simulations in Figure 11. 35 Note that the conditions (a) and (b) are both satisfied with this choice of parameters. As for (c), although R s < R ⋆ for all simulations, for the largest values of 2πα/v dm the ratio R s /R ⋆ is not too much smaller than 1 (the worst case is R s /R ⋆ ≃ 0.8, for 2πα/v 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. 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 ξ foc = 2πα/v dm = 4, 8. This is obtained by taking the spherical average of the field |ψ| 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 34 In practice, to compute the initial conditions we used the plane wave form ψk(x) = e ik·x , which is correct except close tor = 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. 35 We use time-step ∆tmv 2 dm ≃ (∆mv dm ) 2 /4, which is small enough for the energy E to be conserved as described above. Note that a much larger time-step is sufficient for 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πα/v dm . 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 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 ρ(0)/ρ dm at the position of the Sun as a function of time (left) and the density profiles (right) in three simulations with 2πα/v dm = {6.5, 2.5, 1} andg = 0.006 · {1, 1, 0.5}. 36 First, the exponential increase of the density is evident, given the logarithmic scale in the y-axis. Soon after t ≃ τ rel , the overdensity tends to saturate to the expected critical value ρ crit /ρ dm = (2πα/v dm ) 2 /(2π 2g ) 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. Γ > 0, is likely to occur for values of 2πα/v 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 |g| andR 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 ofg andR s , we checked that the time (relative to τ rel ) when the density increases exponentially is independent of these parameters, as long as R s ≲ R ⋆ , which confirms our analytic analysis of Section 4.3. An example of this is the red curve of Figure 21, for whichg 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 .