Dark Photon Stars: Formation and Role as Dark Matter Substructure

Any new vector boson with non-zero mass (a `dark photon' or `Proca boson') that is present during inflation is automatically produced at this time from vacuum fluctuations and can comprise all or a substantial fraction of the observed dark matter density, as shown by Graham, Mardon, and Rajendran. We demonstrate, utilising both analytic and numerical studies, that such a scenario implies an extremely rich dark matter substructure arising purely from the interplay of gravitational interactions and quantum effects. Due to a remarkable parametric coincidence between the size of the primordial density perturbations and the scale at which quantum pressure is relevant, a substantial fraction of the dark matter inevitably collapses into gravitationally bound solitons, which are fully quantum coherent objects. The central densities of these `dark photon star', or `Proca star', solitons are typically a factor $10^6$ larger than the local background dark matter density, and they have characteristic masses of $10^{-16} M_\odot (10^{-5}{\rm eV}/m)^{3/2}$, where $m$ is the mass of the vector. During and post soliton production a comparable fraction of the energy density is initially stored in, and subsequently radiated from, long-lived quasi-normal modes. Furthermore, the solitons are surrounded by characteristic `fuzzy' dark matter halos in which quantum wave-like properties are also enhanced relative to the usual virialized dark matter expectations. Lower density compact halos, with masses a factor of $\sim 10^5$ greater than the solitons, form at much larger scales. We argue that, at minimum, the solitons are likely to survive to the present day without being tidally disrupted. This rich substructure, which we anticipate also arises from other dark photon dark matter production mechanisms, opens up a wide range of new direct and indirect detection possibilities, as we discuss in a companion paper.


Introduction and Summary
Despite the overwhelming evidence for the existence of dark matter (DM), from scales spanning from astrophysical to cosmological, its specific nature remains unknown. Candidates range from primordial black holes or massive super-Planckian composite states, through sub-Planckian particles with mass 1 eV, to particles with mass 1 eV which are best described as semi-classical 'wave' dark matter in galaxies such as the Milky Way. So far, in the particle case, we have no information about the DM spin, and only extremely limited information about its mass and possible non-gravitational interactions with the Standard Model (SM). Among the numerous candidates, a minimal possibility is a new bosonic particle of spin-0 or spin-1. The presence of such particles is expected from string theory compactifications [1][2][3][4][5][6].
Irrespective of their couplings to the SM, elementary spin-0 particles that exist as states at high scales are automatically produced in the early Universe via the so-called misalignment mechanism [7][8][9], and, if stable, form a component or all of the dark matter. There are a variety of other production mechanisms that could lead to a cosmologically interesting relic density of spin-0 particles, though misalignment production is attractive because of its minimality.
For vector bosons, the abundance from misalignment is generically suppressed [6]. However, if a massive vector state is present during inflation, its longitudinal component is automatically produced by inflationary fluctuations [10]. The resulting relic abundance Ω A is where m is the mass of the vector, H I is the Hubble scale during inflation and Ω DM is the observed DM abundance [10]. This expression assumes that the vector mass does not change during evolution of the Universe from the inflationary epoch until today, and that there are no charged or Higgs-like states close in mass or lighter than the vector boson. 1 Given the current upper bound on H I from non-observation of gravitational waves and the requirement that m H I , such production can provide the observed DM density for 10 −5 eV m 10 8 GeV.
The vector is produced because the longitudinal component acts as a scalar during inflation by the Goldstone equivalence theorem and -in the same way as any other scalar present at this time -obtains an approximately scale-invariant spectrum of energy density perturbations. 2 Crucially, after inflation, the perturbations in the vector at large scales redshift faster than they would for a scalar. Consequently, the primordial power spectrum reaches a form that, in addition to the standard adiabatic perturbations, has an isocurvature component peaked at cosmologically tiny scales 10 11 km (10 −5 eV/m) 1/2 , corresponding to the size H −1 of the Hubble horizon when H = m. The suppression of the perturbations at large scales makes the spectrum automatically consistent with bounds on isocurvature perturbations, since the primordial power spectrum is essentially unconstrained at those small distances. Production by inflationary fluctuations is therefore both robust and, at least in the absence of further interactions, unavoidable. 3 The purpose of this paper is to investigate the model-independent dynamical process of structure formation from these small-scale primordial inhomogeneities. We will show that they lead to a strikingly rich structure of gravitationally bound objects, depicted in Figure 1, which are normally absent in con-1 Such a situation is possible in, e.g. the string theory context. 2 The transverse components are not produced since they are equivalent to a massless vector. 3 Very massive scalars are automatically produced from inflationary fluctuations by a different process [11]. Figure 1: A schematic representation of the dark photon dark matter distribution today. At subgalactic scales, roughly parsec and smaller, the dark matter is bound in compact halos that arise from the collapse of small-scale primordial inhomogeneities (centre panel). The left panel shows structure at smaller scales: solitons, i.e. dark photon stars, their surrounding fuzzy halos, and dark matter filaments connecting them.
ventional cold dark matter structure formation. 4 As we will see in Section 3, the typical length scale of the inhomogeneities is so tiny that the quantum pressure of the bosons is relevant in the dynamics. Indeed, although the small-scale inhomogeneities are already non-linear before matter-radiation equality (MRE), quantum pressure prevents their collapse until after MRE. At that point, the overdensities collapse into gravitationally bound objects fully supported by quantum pressure, with typical mass of order 10 −16 M (10 −5 eV/m) 3/2 . The number of bound particles is N 10 55 (10 −5 eV/m) 5/2 , and for all vector masses m few × 10 7 GeV the occupation numbers of the quantum states are so large that this bound object is well-described as a quantum soliton of the semi-classical vector field.
We find that these dark photon star or Proca star solitons make up an order one fraction of the DM, a remarkable possibility given that the solitons and both macroscopic and intrinsically quantum. Moreover they attract surrounding DM, which therefore gets bound around the solitons in a 'fuzzy' halo, see Figure 1 (left). Additionally, the solitons are produced with excited quasi-normal modes, which are long lived and decay via the emission of dark photon (spherical) waves.
The evolution during inflation also leaves the vector field with overdensities at larger length scales, with a k 3 power spectrum. Overdensities on increasingly large scale become non-linear after MRE, and when they collapse they induce a 'primordial' structure formation, producing small halos, which we dub compact halos. These have mass of order 10 −11 M (10 −5 eV/m) 3/2 and will contain some of the solitons, see Figure 1 (centre). The compact halos will become subhalos of the larger galactic halos that are produced when the standard adiabatic almost-scale invariant fluctuations sourced by the inflaton become non-linear.
We will see that the solitons, fuzzy halos and the densest compact halos are likely to remain undisrupted during the formation and evolution of the Milky Way, and therefore are likely to persist to the present day. The dark matter substructure in Figure 1 is an unavoidable property of vector boson DM produced by inflation fluctuations. It could therefore lead to smoking-gun signatures for dark photon dark matter, singling out its particle nature and production mechanism. This is possible even in the most challenging scenario in which the dark photon has no non-gravitational interactions with the SM. The dark matter substructure can be investigated with gravitational-only probes. Depending on the mass of the substructures, the possibilities include, at minimum, pulsar timing arrays [13][14][15][16], microlensing, [17][18][19][20][21][22][23][24] photometric microlensing [25][26][27][28], and extra-galactic strong gravitational lensing [27]. In particular, we will see that the size and mass of the solitons and compact halos are in one-to-one correspondence with the dark photon mass. As a result, experimental evidence of such DM substructure will allow us to infer the mass of the dark photon. This would lead to a prediction of the Hubble scale during inflation, which -if confirmed e.g. by observation of tensor modes in the cosmic microwave background -would be compelling evidence for this type of dark matter.
In the presence of direct interactions of the massive vector with the SM, e.g. kinetic mixing with hypercharge and thus the SM Z-boson and photon [29][30][31][32][33], the substructure leads to a plethora of signatures in both direct and indirect detection experiments. Solitons (and the fuzzy halos around them) have typical energy density of the order of the Universe's density at MRE, many orders of magnitude larger than the local DM density in the vicinity of the Earth today. Similarly, compact halos have average density of a few orders of magnitude larger than the local density. We will see that solitons can encounter the Earth and other astrophysical objects frequently. These will lead to significant changes to direct detection prospects and new astrophysical and cosmological signals. In a companion paper we discuss the resulting detection and observational signals in detail. 5 Additionally, the primordial perturbations of a vector produced from inflationary fluctuations have surprisingly similar qualitative features to those in other new physics scenarios. For example, axions in the post-inflationary scenario have qualitatively similar primordial perturbations, although the dynamics that produces them is totally different. There are also other model-dependent mechanisms that could lead to a relic abundance of vector dark matter [36][37][38][39][40][41][42] (in the same way as there are other production mechanisms besides misalignment for scalars). These might also lead to similar initial conditions, e.g. if the vector is produced via parametric resonance or topological defects. Therefore, although in this paper we focus on a particularly minimal and predictive theory, we expect that our approach and results may be at least qualitatively applicable to a wide range of scenarios.
The paper is structured as follows: In Section 2 we describe the production of a vector boson during inflation and calculate the power spectrum of primordial density inhomogeneities during radiation domination, following [10] (see also [43][44][45] for related work). A reader familiar with [10] or solely interested in the later development of small-scale structure and solitons may safely skim rapidly through this Section, the primary results being eqs. (9) and (10) along with the precise numerically-derived spectrum shown in Figure 2 (left panel). In Section 3 we first discuss the post-inflation dynamics of the fluctuations and the importance of the 'quantum pressure' term arising from the Heisenberg uncertainty principle. We then highlight aspects of the physics of the vector solitons that are particularly important for phenomenology, before discussing in detail the non-perturbative evolution of inhomogeneities around the time of MRE and the formation of solitons from their collapse. In Section 4 we describe the 'primordial' structure formation at larger scales and the compact halos. In Section 5 we discuss the survival of the vector dark matter substructure until today and the collision rate of solitons with the Earth. In Section 6 we summarise our results, describe improvements and extensions, and propose future directions, some of which will be covered in our companion paper which focuses on the potential observational and experimental implications of the solitons and their fuzzy halos. Appendices provide details on the initial conditions from inflation, analytic analysis of the evolution of overdensities, our approach to numerically solving the Schrödinger-Poisson equations, analytic treatment of the soliton and compact halo mass distributions, and finally an extensive discussion of the survival of solitons, fuzzy halos and compact halos to the present day.

Initial Conditions from Inflation
Consider a vector boson A µ (we will also use the name 'dark photon' interchangeably) with field strength F µν , described by the action with metric ds 2 ≡ g µν dx µ dx ν . We remain agnostic about the dynamics giving rise to the mass as long as m remains constant during and after the inflationary epoch, and that no other light fields significantly coupling to A µ are present. All that matters in the following is that the action in eq. (2) describes the vector field during inflation and in the subsequent evolution of the Universe. 6 As we will see, efficient inflationary production requires m H I . 7 For the purpose of determining the relic abundance it is sufficient to consider the homogeneous background FRW metric ds 2 = −dt 2 + a 2 (t)d x 2 .
The three propagating degrees of freedom of the vector field are A ≡ A i , while the component A 0 does not have a kinetic term in eq. (2) and corresponds to an auxiliary field. We can eliminate A 0 by writing eq. (2) in terms of the Fourier modesÃ µ = d 3 x exp(−i k · x)A µ , where k is the comoving momentum (we will drop the tilde in the following). The equations of motion of A 0 become algebraic, can be solved explicitly, and their solution plugged back into eq. (2) to get rid of A 0 . This leads to S = S T + S L with where we have decomposed A in terms of the longitudinal and transverse modes A L and A T , defined by k · A = kA L and k · A T = 0. The actions S T and S L describe transverse and longitudinal modes, which are decoupled from each other given the FRW form of the metric. Eqs. (3) and (4) hold both during inflation -when the Hubble parameter H ≡ȧ/a = H I is approximately constant -and subsequently in the early Universe (although we will see that the transverse and longitudinal modes become coupled together around MRE).
We consider the evolution of a generic mode, starting from when it is subhorizon during inflation (k/a > H I ) in the vacuum state to when it becomes nonrelativistic (k/a < m) and subhorizon, in the radiation dominated era. First note that all the modes that start subhorizon during inflation (k/a > H I ) are relativistic during inflation, since k/a > H I m, and remain relativistic at horizon exit (i.e. when a = k/H I ). Using conformal time dη = dt/a, the action for the transverse modes reads S T = (2π) −3 d 3 kdη 1 2 |∂ η A T | 2 − (k 2 + a 2 m 2 )| A T | 2 . Since m is negligible for relativistic modes, S T is 6 Swampland conditions might require m 0.3 eV for such an effective theory to be embedded in a UV completion [46], however there are constructions that claim to evade these limits [47]. Moreover, we will see that the case m eV leads to interesting phenomenology and experimental signals. 7 If the mass m is produced by a dark Higgs mechanism, our analysis applies if the mass of the dark Higgs is much larger than the inflationary Hubble scale, so it can be neglected during inflation. time-translation invariant (in fact, conformally invariant). The vacuum state of the transverse modes therefore does not change, and they are not produced. We therefore set A T = 0 in the reminder of this Section. 8 On the other hand, in the relativistic limit the action for the longitudinal modes reduces to that of a free real scalar ϕ ≡ (m/k)A L , i.e. S L = a 3 d 3 xdt 1 2 [(∂ t ϕ) 2 − |∇ϕ| 2 /a 2 ], where we Fourier transformed back to coordinate space. It is well known that the vacuum of this theory is time-dependent [48,49]. Modes that are in the vacuum at early times get populated -after they exit the horizon -with a Gaussian power spectrum P ϕ = (H I /2π) 2 , where we defined the power spectrum P X of a generic field X as This implies that P A L = (kH I /2πm) 2 at horizon exit, a = k/H I . The expectation value of the energy density ρ ≡ T 00 reads and the energy density spectrum at horizon exit is scale invariant, ∂ρ/∂ log k H 4 I /(2π) 2 . This is the well known result for a scalar field, and indeed the longitudinal component of the vector in the relativistic limit reproduces a massless scalar by the Goldstone Equivalence theorem.
After inflation the modes evolve classically, following the equations of motion of eq. (4), We can calculate the evolution of each mode with initial condition A L = A L,0 andȦ L 0 at a = k/H I , until it reenters the horizon and becomes nonrelativistic. Since eq. (7) is linear, A L will still have a Gaussian distribution during such evolution, with a power spectrum at a generic time given by where we used eq. (5). 9 Unfortunately, the solution for A L cannot be evaluated analytically for a generic k. However, it is possible to understand the behaviour of A L (t, k)/A L,0 analytically. Here we summarise the main results, leaving complete derivations to Appendix A. All the modes of interest will eventually become subhorizon and nonrelativistic, and it is convenient to classify them into two classes: (1) Those that become nonrelativistic while still superhorizon ('low frequency') and (2) Those that reenter the horizon while relativistic ('high frequency'), and become nonrelativistic only afterwards. We define the mode that becomes nonrelativistic exactly when it reenters the horizon k /a ≡ m = H(a ), so the low (high) frequency modes satisfy k < k (k > k ) respectively.
While relativistic and superhorizon ρ ∝ a −2 for all modes. Modes with k < k are suppressed because they become nonrelativistic while superhorizon and subsequently their energy density decreases as ρ m 2 A 2 L /a 2 ∝ a −2 until they enter the horizon. This is the crucial difference compared to a scalar field, for which ρ is frozen for nonrelativistic superhorizon modes. The difference is due to the form of the mass term, which controls the energy density in such modes: for a scalar, 1 2 m 2 ϕ 2 ∝ const, while 1 2 m 2 g ij A i A j ∝ a −2 for a vector. Meanwhile, the modes with k > k are suppressed because they enter the horizon while still relativistic. They have ρ ∝ a −4 after they enter the horizon but before they become nonrelativistic and ρ ∝ a −3 subsequently. The result of this is that the spectrum is peaked at momentum k . Modes with k k are suppressed since they stayed in the superhorizon nonrelativistic regime the longest, and those with k k are suppressed because they underwent subhorizon relativistic redshift before becoming nonrelativistic, with larger k suppressed more since they were subhorizon and relativistic for longer. The power spectrum of A L is, to a very good approximation, given by The exact form of P A L (t, k), plotted in Figure 2 (left), can be extracted by solving eq. (7) numerically. The least suppressed mode, k , corresponds to subgalactic scales today: defining λ ≡ 2π/k where a 0 is the FRW scale factor today. Moreover, the misalignment mechanism is related to the energy density in the zero mode, which gets a huge suppression, and is therefore ineffective. The energy density of the vector behaves as matter at late times, and forms a component of the DM abundance. Given the peaked form of P A L , the DM abundance is approximately given by the energy density in modes of momentum k , redshifted from horizon exit at a = k /H I to a = a and then to today when the scale factor is a ≡ a 0 . Therefore, ρ(a 0 ) . A full calculation leads to the relic abundance given in eq. (1). The Hubble scale during inflation is bounded by the non-observation of tensor modes in the cosmic background, assuming single field slow roll inflation. The latest data combination from Planck and BICEP2 [50] bounds H I < 6 · 10 13 GeV. 10 As a result, the vector can account for the full dark matter abundance if m 10 −5 eV.
The structure in the power spectrum of A L leads to small-scale overdensities in the energy density field ρ = (Ȧ 2 L + m 2 A 2 L )/2. During radiation domination and once all the relevant modes have become non-relativistic, the properties of typical fluctuations remain constant. In Figure 3 we plot a section of the vector's energy density over a line at this stage, normalised to its average energy density ρ ≡ ρ . There are obvious O(1) fluctuations in the energy density (the peaks have ρ/ρ 2 ÷ 4), with spatial size of order λ .
It is convenient to introduce the overdensity field The distribution of inhomogeneties is encoded in the power spectrum of overdensities, P δ , defined from δ(x) as in eq. (5). Once all the relevant modes are non-relativistic, P δ is constant during radiation domination. Using the fact that ∂ t A L and A L are independent Gaussian fields, one finds Figure 2: Left: The power spectrum of the longitudinal vector field component as is automatically produced by inflationary fluctuations, during radiation domination and after the modes have become subhorizon and nonrelativistic. The spectrum is peaked at k , corresponding to the momentum equal to the Hubble parameter H when H = m. Right: The power spectrum P δ of the (non-Gaussian) overdensity field. Importantly for later soliton formation, at matter radiation equality the quantum Jeans momentum k J (ρ) associated with the mean DM density (defined in Section 3.1) coincides with k up to an order one factor, regardless of the dark photon mass.
which is a useful analytic approximation that captures the asymptotic limits k/k → {0, ∞} exactly (see Appendix A and [10] for full expressions). In Figure 2 (right) we plot P δ , which as expected is peaked at k/k 1, and decreases as (k/k ) 3 and (k /k) at small and large k, respectively. Since δ(x) does not have a Gaussian distribution (indeed, it is asymmetric around δ = 0) it is not fully described by P δ . 11 The power spectrum however still provides useful information about the variance of the field and the magnitude of the overdensities. Note that at length scales much larger than λ = 2π/k , δ(x) is Gaussian and can be fully reconstructed from P δ alone.
The fluctuations in the vector's energy density are isocurvature perturbations, since they are induced only in the vector during inflation. As mentioned in the Introduction, these fluctuations are only allowed to be O(1) because perturbations at much larger scales, which are strongly constrained by observations, are automatically suppressed thanks to the k 3 behaviour. This is in contrast to a scalar field, for which the power spectrum from inflationary fluctuations is flat at k < k , and order one fluctuations are completely excluded unless the scalar is a tiny fraction of the total DM. The smallest scales at which the power spectrum has been observed are roughly k obs = 7Mpc −1 from Lyman-alpha [51], so k obs /k 10 −11 ( eV/m) 1/2 , and the observed modes are far off the left of the plot in Figure 2 (right) for all relevant dark photon masses.
Inflation also sources perturbations of the inflaton, which are metric perturbations with a change of gauge, i.e. ds 2 = −(1 − 2Φ)dt 2 + (1 + 2Φ)a 2 d x 2 . The gravitational potential Φ has small differences in different patches after inflation. These lead to the same relative perturbations in all form of energy, i.e. adiabatic perturbations, including in the vector overdensity field δ(x). Its power spectrum P δ therefore automatically acquires also the almost-scale-invariant contribution, as shown by the purple line in Figure 2 (right), as is necessary to be consistent with observations.
We have considered a vector with only the action of eq. (2), in which case the relic density discussed Figure 3: A one-dimensional slice of the vector field energy density ρ, during radiation domination, relative to its mean density ρ. The typical fluctuations are on scales λ = 2π/k , eq. (10), corresponding to the P δ peak in Figure 2. There are also small overdensities at scales larger than λ (not distinguishable in the plot), which correspond to the (k/k ) 3 part of P δ , and considerable fluctuations on smaller scales because of the slow (k /k) fall-off of the power spectrum for k k .
in this Section is an unavoidable contribution to the dark matter abundance. The situation is more complicated if the vector has self-interactions, couples to other particles or has a non-minimal coupling to gravity [6]. Such possibilities do not necessarily affect the production during inflation, and indeed it would be surprising if fluctuations of order H I were prevented. However, the subsequent evolution will be affected in some theories, which could alter the relic abundance. For instance this is the case if there are 'dark electrons' [52]. 12 A key direction for future work is to systematically study the impact of different possible interactions.

Collapse of Inhomogeneities: Vector Solitons
During radiation domination the vector field evolves freely, i.e. as in the absence of a gravitational potential. 13 The gravitational potential becomes relevant at around MRE. If the vector constitutes a sizeable fraction of the dark matter abundance, the large overdensities corresponding to the peak of P δ at k , start evolving nonlinearly under the effect of the gravitational interactions at around this time (more precisely, at a a eq /δ where a eq is the scale factor at MRE, see Appendix C.1 for further details).
In what follows we will assume that the vector makes up all the DM. Since they are already O(1), as soon as the gravitational potential becomes important the overdensities cannot be studied using the standard perturbative treatment of the density field. They are expected to clump into bound objects, and -as we will see next -this happens in the regime where the fuzzy dark matter properties of the boson, i.e. quantum pressure, are important. 12 If an interaction leads to an operator of the form F 4 /Λ 4 in the effective theory at H , when the most important modes re-enter the horizon, the subsequent evolution is unaffected provided Λ (mHI ) 1/2 (m/ eV) 3/2 MeV where the last equality holds if ΩA = ΩDM. We leave the effects on the evolution of the superhorizon modes for future work. 13 This is the reason we neglected Φ in the previous Section.

Dynamics of fluctuations and quantum pressure
The modes relevant for the DM abundance are nonrelativistic at around MRE. It is therefore convenient to rewrite the equations of motion (EoM) of the Lagrangian in eq. (2) and the component 00 of the Einstein equations in their nonrelativistic form. To do so, we work in terms of ψ i defined by in the limit where ψ is slowly evolving soψ i mψ i andψ i m 2 ψ i . In terms of ψ i the EoM become the Schrödinger-Poisson (SP) system [53] Apart from involving the three propagating components of A, these equations have exactly the same form as those for a nonrelativistic scalar field. (The fourth vector field component, A 0 , is non-dynamical and may be recovered from the Lorenz-Proca constraint ∇ µ A µ = 0 which follows as a consistency relation from the EoM.) At around MRE when Φ is non-negligible, the EoM are nonlinear and couple together all the components of the vector through the gravitational potential. In particular the longitudinal and transverse modes are coupled and no longer evolve independently. Consequently, although initially vanishing, the transverse mode is sourced by the longitudinal mode via Φ. 14 For a full treatment, eqs. (14) and (15) need to be solved numerically with initial conditions in eq. (9). However we can qualitatively understand the evolution by doing the Magdelung transformation, i.e. writing eqs. (14) and (15) in terms of density and velocity fields, ρ i and v i , defined by ψ i = a 3 ρ i e iθ i with v i ≡ (am) −1 ∇θ i . The imaginary and real parts of the Schrödinger equation become the continuity and Euler equations of a three-component perfect fluid with local density ρ i = |ψ i | 2 /a 3 and velocity v i , i.e.
where ρ = i ρ i and ρ ≡ ρ is the average energy density of the vector field, and we defined where for clarity we have restored the factors of in this one expression. Their appearance is due to the fact that the mass m of a particle is independent of while, due to the that appears in the Planck-Einstein-de Broglie relation, the mass parameter that appears in the field action, eq. (2), is really m/ .
From now on we return to = c = 1 units. From eq. (17), the fluid is subject to the gravitational potential Φ and to the 'quantum pressure' potential Φ Q . The gravitational potential tends to increase the overdensities, while the quantum pressure tends to make overdensities fluctuate. This can be seen from the fact that the quantum pressure has the opposite sign as the gravitational potential, or more rigorously in the perturbative treatment of small overdensities, which we summarise in Appendix B.
The importance of Φ Q with respect to Φ can be understood by comparing the last two terms of eq. (17), i.e. if ∇Φ ∇Φ Q the quantum pressure dominates. Taking the divergence of this relation and using eq. (18) we get −8πGρm 2 a 4 ∇ 2 (∇ 2 √ ρ/ √ ρ). Going to Fourier space, for comoving momenta k much larger than the comoving 'quantum' Jeans momentum associated with physical density ρ Φ Q is much more important than Φ, and in the opposite limit it is irrelevant. This implies that overdensities at comoving length scales smaller than λ J ≡ 2π/k J are prevented from collapsing, while those at length scales much larger than λ J evolve as in the absence of quantum pressure. Note that this is not the conventional Jeans scale, which is proportional to the inverse of the sound speed and in this context is infinity.

The remarkable coincidence
Crucially, if the vector boson makes up an O(1) fraction of the dark matter, at the time of MRE the quantum Jeans scale k J (ρ) (corresponding to the average density) is parametrically close to the typical scale of spatial fluctuations k , independently of m. In particular, the ratio between k J (ρ) and k at MRE is where Ω M is the matter energy density normalised to the critical density, and ρ tot eq is the total matter energy density at MRE. In eq. (21) we defined g R ≡ (g eq, /g , ) 1/4 (g ,s /g eq,s ) 1/3 , where g s , g denote the effective number of relativistic degrees of freedom for entropy and energy, and as usual eq and denote quantities at MRE and when H = m. This factor accounts for the change in the number of degrees of freedom between H and MRE, which affects the value of k . The cancellation of m, G and ρ eq in eq. (21) occurs because, up to numerical factors, a /a eq T eq /T Gρ tot eq /m 2 1/4 (making the excellent approximation that only radiation contributes to the energy density at T ). 15 For dark photon masses of interest, m 10 −5 eV, T 200 GeV, so we set g ,s and g , to their Standard Model high temperature values, which gives g R 1.27. 16 Therefore, if the vector comprises the full dark matter abundance, Eq. (22) means that the quantum pressure term (i.e. the wave-like properties of the boson) cannot be neglected and will affect the evolution and collapse of the O(1) overdensities in Figure 3, regardless of the value of m, and even of G. In Figure 2 we show the value k eq J of k J (ρ) at MRE. In fact, the numerical value is close to the peak of P δ , which is at about 2k . 15 More explicitly, a /aeq = g −1 R (4πGρ tot eq /3m 2 ) 1/4 . We also note for future use that a /aeq can be rewritten as a /aeq = g −1 R (Heq/( √ 2m)) 1/2 , using H 2 eq = 8πGρ tot eq /3. 16 This is a reasonable assumption provided there are not 10 2 new degrees of freedom close to the TeV scale.

Dynamics around matter radiation equality
Eq. (22) leads us to conclude that at a = a eq : • For modes with k k eq J k the dynamics is as in the absence of quantum pressure. In particular, over cosmological scales -much larger in length than λ -the quantum pressure is irrelevant and the field behaves as conventional DM. This applies to the adiabatic modes and those in the k 3 IR tail of P δ , see Figure 2. Since the fluctuations in the k 3 tail are small at MRE, at least initially the field at these distances is in the perturbative regime and its evolution can be calculated analytically. Deep in radiation domination such isocurvature fluctuations (i.e. the k 3 modes) are frozen, and once in radiation domination they grow linearly. To a good approximation δ ∝ 1 + 3 2 a aeq and P δ ∝ (1 + 3 2 a aeq ) 2 (we give further details in Appendix B).
• The dynamics of modes with k k eq J is dominated by the quantum pressure. As mentioned, some of these modes are already nonlinear, but are prevented from collapsing into bound objects, and instead just oscillate.
• Modes k k eq J are on the boundary at which quantum pressure is relevant at MRE. Owing to the coincidence between k k eq J and k , these are at the peak of P δ and there are order one fluctuations on these scales. Such modes oscillate because of quantum pressure before MRE, and -as we will see next -collapse around MRE.
The fact that the modes with k k are linear allows us to (at least initially) treat the evolution of the modes with k k independently. As time increases the comoving quantum Jeans momentum increases as k J (a) = k eq J (a/a eq ) 1/4 , see eq. (20), i.e. the dashed line in Figure 2 moves to the right (here and in the following, whenever k J takes a scale factor a as an argument we mean k J (ρ(a))). As a result, the modes at larger comoving momentum, previously prevented from collapsing, will be able to clump and are expected to form compact bound objects. In other words, collapse of the overdensities with k k is prevented until the comoving quantum Jeans length is smaller than the comoving size of the overdensity, when they start being dominated by the gravitational potential instead of the quantum pressure. In particular, given the coincidence in eq. (22), the modes at k k (where P δ is peaked) collapse into bound objects at around a a eq , when the two independent factors preventing their collapse (quantum pressure and radiation domination) both pass. 17 Once formed the objects rapidly decouple from the Hubble flow and will be described by a stationary solution of eqs. (14) and (15) in the absence of expansion. 18 Since the collapse happens at a scale where the quantum pressure is still relevant in the dynamics, we expect that for such solutions the last term in eq. (17) is of the same order as the others, in particular of ∇Φ. In fact, as we will see shortly, these bound objects are solitons. For these, the gravitational potential term is fully balanced by the quantum pressure term (Φ Q = −Φ), instead of being balanced by the velocity terms as in a conventional halo. Moreover, we expect that the mass of the objects is parametrically set by the dark matter mass inside a region of the size of the collapsing perturbation, which is given by the quantum Jeans scale at the time of collapse, so a soliton formed at time a has mass and c M a dimensionless time-independent O(1) coefficient. As time increases, M J (a) decreases, and solitons with smaller and smaller masses are produced, from the collapse of the smaller scale fluctua-tions. When integrated over a > a eq , this leads to a relic abundance of solitons with a nontrivial mass distribution. As we will discuss in detail in Section 4, as time progresses also the modes with k k become nonperturbative and collapse to form halos, which unlike solitons are supported by the velocity term in eq. (17). We dub these compact halos. These are expected to include some of the solitons previously produced. Inside the compact halos the modes with k k can no longer be treated separately, as they are no longer decoupled. However, outside the compact halos the field at k k is still perturbative and solitons are still expected to form following eq. (23). Once the majority of the DM is bound in compact halos at around a/a eq 30, z 100 (see Figure 18 in Section 4), the production of solitons is expected to decrease, eventually approaching zero.
Note that, crucially, it is thanks to the coincidence in eq. (22) that (dense) solitons form. If P δ were peaked at a momentum k much smaller than k eq J , the overdensities would have collapsed into halos, as quantum pressure would have been negligible. If instead k eq J k , the collapse would have been prevented until late times (if at all) and resulted in much less dense solitons (since their density is set by the DM density at the time of formation, as discussed below).

Vector Solitons (Quantum Dark Photon Stars)
A vector soliton is the self-gravitating stationary solution of eqs. (14) and (15) in (asymptotically) flat spacetime a = 1 that minimises the total energy E at fixed, finite, value of the vector-boson particle number This particle number is conserved as the fundamental action for the massive vector boson, eq. (2), contains no number-changing self-interactions. (The implied interactions with gravitons, h µν , following from the action are ineffective as they do not lead to number-changing processes such as A i A j → h, 2h, . . . for any of the solutions we consider.) 19 In any case, the validity of a semi-classical analysis in terms of fields requires N 1. Since we will always be in a regime where the gravitational binding energy is small compared to the total bare rest-mass M ≡ mN it is convenient for later purposes to use M as an effective conserved mass parameter of the soliton. In the limit in which we work, and in terms of M , the total soliton energy is Although not mathematically proven, it has been conjectured that ground state vector solitons have the form ψ i ∝ u i where u i is a generic complex vector with unit norm (u * i u i = 1), and have spherically 19 Note that a free theory possesses an infinite number of conserved quantities as the particle and antiparticle numbers for every single wavenumber and spin direction are individually conserved. In the presence of classical gravitational interactions a only a finite number of these quantities, including N , remain conserved due to non-linear mode mixing. In the presence of additional interactions with, e.g., Standard Model fields, the effective conservation of particle number over the cosmological timescales of interest must be reconsidered. We emphasise that contrary to statements in the literature the existence of a dark photon, or equivalently Proca star does not require Aµ to be promoted to a complex field. This conservation law is made explicit in the non-relativistic reduction, eqs. (13), (14) and (15), as the fields ψi are complex (despite the fact that Ai are real) and the effective action for ψi associated to these SP equations has a global U(1) invariance leading to a conserved particle number. Thus vector solitons are an example of a non-topological Q-ball soliton [54] stabilised against decay due to their being the lowest energy state at fixed N , in particular lower than the energy E = N m of a collection of N far-separated vector bosons. symmetric energy density [53,55]. 20 In this case, the solution space of ground state solitons as a function of their mass can be constructed by a rescaling of the basic ansatz where r ≡ | x| is the distance from the soliton's centre and γ a fixed real constant. Here χ 1 (x) and Φ 1 (x) are monotonic functions that satisfy χ 1 (x) → (1, 0) and Φ 1 (x) → (const, 0) for x → (0, ∞). Numerical solution of the equations arising from the use of this ansatz shows that for the ground state soliton, γ −0.65, while an excellent analytic approximation for . Given the form of χ 1 , the soliton's energy is localised at its centre. All other ground state soliton solutions of eq. (14) and (15) have the form in eq. (26) with the rescalings Given eqs. (25) and (26) and A i ∝ ψ i e −imt , the quantity −γα 2 can be interpreted as the binding energy per unit mass. Our non-relativistic, weak-gravitational field, approximation thus requires α 1 for consistency. 21 The corresponding rescaled soliton solution has mass and half-mass radius respectively. All the soliton solutions that we consider in this work have α 1 and are stable. It will be useful to us that the combination M R is independent of the soliton mass and is given by The corresponding central density of the solitons is The quantum Jeans length associated to this density is λ J (ρ s ) 2.3R, and is therefore of order of the size of the object. When studying the properties of solitons here and subsequently λ J (ρ s ) = 2π/(16πGρ s m 2 ) 1/4 refers to the physical rather than comoving quantum Jeans length. As discussed, this means that the quantum pressure term is relevant for this solution. In fact, since v = 0 from eq. (26), solitons are configurations fully supported by quantum pressure, in the sense that Φ Q = −Φ (this can be used as an alternative equivalent definition of solitons). In Figure 5 (left) we plot the soliton density profile in terms of λ J (ρ s ).
Solitons can be thought of as bound configurations of bosons with gravitational energy balanced by their intrinsic kinetic energy. This comes from their 'irreducible' velocity, of order 1/(mR), that stems from the uncertainty principle associated to the Schrödinger eq. (14) and the fact that the particles are localized in a finite region of space. Indeed, the de Broglie wavelength of the particles is approximately equal to the size of the soliton. Note that given eq. (28), this intrinsic velocity is the order of the virial velocity (GM/R) 1/2 that the particles would have in a gravitationally bound configuration, consistent with the particle interpretation.
Given our expectation that the masses of cosmologically produced solitons are determined by the quantum Jeans length at the time of creation, eq. (23), the density of the solitons is parametrically set by the average DM energy density when they are created, independently of m. Indeed, plugging the expected value of the solitons masses in eq. (23) into eq. (29), the density of the solitons produced at time a is ρ s = 4.5 · 10 4 c M ρ(a).
It is noteworthy that ground state vector solitons with a fixed mass form an infinite set parameterized by u i , i.e. a generic direction in coordinate space, which indicates the direction of the vector field. This moduli space of ground state solutions is due to the U (3) symmetry of eqs. (14) and (15). For a fixed direction, the soliton profile resembles the well known scalar soliton, see e.g. [57].
We also observe that a soliton is only an exact solution in the idealised case where the total mass d 3 x i |ψ i | 2 is finite. When the mass is infinite, for instance in the presence of a constant density background ρ, solitons of central density ρ s ρ are still good approximations of stationary solutions near their cores. However, away from their centres the solitons' profiles will be deformed, and, since they gravitationally attract the background matter, they will be surrounded by a halo.
Importantly, as argued in detail in Ref. [55], depending on the choice of u i , and in the absence of other interactions, the ground state solitons can either carry zero total angular momentum J, or possess macroscopically large J without energy cost or alteration of the radial profile χ 1 (r), at least at leading order in an expansion in G. Since all known physically-realised astrophysical objects carry angular momentum, it is interesting to discuss in more detail the physics of these J = 0 vector solitons.
The total angular momentum J is composed of two parts, an orbital angular momentum contribution, L, and an intrinsic spin density which integrates to S. In terms of ψ i , and in the non-relativistic weakgravity limit in which we are working, they are given by On the ground state soliton solutions, eq. (26), the first of these expressions gives L = 0. On the other hand, the intrinsic angular momentum of these solitons is S = i( u × u * )N , where once again we have temporarily restored the factors of for clarity. Since | u × u * | ≤ 1, S can acquire any value in [0, N ] and can even vanish depending on u i . This can be made explicit by employing an orthonormal basis, n , for the λ = ±1, 0 spin quantization states along a fixed axisn. Choosing, without loss of generality,n =ẑ these polarisation states are Then one may expand ψ i and thus A i in this basis; specifically for the L = 0 solutions where in principle the functions χ (λ) ( x) could differ. However in the non-relativistic weak-gravity limit of interest to us, χ (λ) ( x) = a (λ) χ 1 (mr) with a (λ) complex coefficients satisfying λ |a (λ) | 2 = 1, and with χ 1 (mr) satisfying exactly the same equation as before. 22 For instance, for |a (±1) | = 1 the resulting total spin angular momentum parallel to theẑ axis is maximum and reads J z = ±N . Linearly or partially polarised ground state solitons are possible if non-trivial complex linear combinations of the angular momentum eigenstates result from the non-linear dynamics of structure and vector soliton formation (in the absence of additional interactions the decoherence time can be parametrically long). In fact, we will see in the next Section that, in the minimal theory we study, the solitons are generically formed with a nontrivial polarisation. In particular, they have intrinsic angular momentum uniformly distributed in the range S ∈ [0, N ], and, as expected, negligible orbital angular momentum L. This is not surprising given that, as mentioned, there is no energy difference between solutions with different intrinsic angular momentum. We note that, if present, interactions of the dark photons with the Standard Model, themselves or other new light states could lead to solitons with a particular angular momentum being energetically favoured.

Comparison with numerical simulations
To confirm our analytic expectations, measure the unfixed numerical coefficient c M in eq. (23), and determine the soliton mass distribution, we study the dynamics of the system numerically. It is straightforward to show that the equations of motion (14) and (15), the initial condition in eq. (9) and a(t) only depend on the number of relativistic degrees of freedom at T and Ω A /Ω DM (see Appendix C.1). The evolution is therefore independent of the values of m and G, which is why the ratio between k and k eq J in eq. (22) does not depend on m explicitly. We continue to fix the number of relativistic degrees of freedom at T to the SM high temperature value and we assume that the dark photon makes up all the DM.
We solve eqs. (14) and (15) numerically on a discrete lattice in a periodic box of constant comoving size (3.75λ ) 3 , starting deep in radiation domination at a/a eq 0.01 (this choice is not important as long as a/a eq 1), from a realisation of the initial conditions with the Gaussian power spectrum in eq. (9). The final simulation time is limited by loss of resolution of the soliton cores and by the growth of density perturbations on the scale of the box, which once non-perturbative lead to finite volume systematic uncertainties. With our available numerical resources we can run to a/a eq 7. Our main results are obtained averaging over approximately 100 individual simulations. Further details of the evolution and tests of the systematic uncertainties are given in Appendix C.
To verify that the broad features of the dynamics of Section 3.1 do indeed occur, in Figure 4 (left) we show the time evolution of the maximum density ρ max in a single, typical, simulation run. As expected, prior to MRE ρ max /ρ is, on average, constant. Due to quantum pressure, density fluctuations on small scales, k k J (a), oscillate even during radiation domination, which leads to small oscillations in ρ max /ρ during this time. This would not happen in the limit k J /k → ∞, in which case ρ( x) would be almost completely frozen for a nonrelativistic field. 23 In the absence of quantum pressure, the largest overdensities in the initial conditions δ 3 would collapse at a/a eq 1/δ 1/3. As expected, collapse is actually delayed to a later time, a a eq , when k J (a) has exceeded k . At a > a eq , the maximum density increases fast, until it reaches an approximately constant value, while the mean dark matter density continues to decrease. This indicates that a bound object, in our case a soliton, has formed and decoupled from the Hubble flow. Additionally, the maximum density has clear oscillations, which are due to the soliton forming with excited quasi-normal modes. We Figure 4: Left: The evolution of the maximum value ρ max of the energy density field of the vector through MRE, in a simulation with volume (3.75λ ) 3 . Time is parameterised by the scale factor a relative to that at MRE, a eq . The mean vector energy density ρ is also plotted. At early times ρ max follows ρ, with small fluctuations due to the oscillation of modes with k k J , driven by quantum pressure. The collapse of overdensities with δ 1, which in the absence of quantum pressure would occur at a/a eq 1/δ, is hindered until after MRE. Once k J /k ∝ a 1/4 has grown sufficiently, overdensities collapse. After the collapse, the maximum density is at a point inside a soliton. The soliton is produced with excited quasinormal modes, so the maximum density subsequently oscillates. Right: A slice of the energy density at a/a eq = 7, in the same simulation as is plotted in the left panel. The slice passes through the point that has the largest density at this time, which is at the centre of a soliton. The soliton (red region in inset) is surrounded by a spherical 'fuzzy' halo (yellow/green region) and there are cosmic filaments connecting it to other solitons. Spherical waves can be seen around the soliton, which are due to the emission of energy from quasinormal modes. A video showing the evolution can be found at [58]. study the growth of density perturbations and the evolution of the density power spectrum in more detail in Appendix D.
In Figure 4 (right) we plot the density field ρ through the slice of the same simulation that contains the point with the largest density, at a/a eq = 7. There is a central soliton (red region). The soliton is surrounded by a spherical fuzzy halo (yellow/green region) extending far from its core, the maximum density of which is about two orders of magnitude smaller than the soliton core density. Finally, the early stages of a cosmic web connecting different solitons have formed (see also Figure 1 left, where we show a 3D version of the same energy density). Spherical waves can be seen beyond the halo. These are due to energy released by the decay of the soliton's quasi-normal modes.
To understand the nature of the collapsed objects, in Figure 5 (left) we plot the spherically averaged density profile around the centre of the objects at a/a eq = 5, averaged over all the objects in our full set of simulations. To enable the profiles of objects with different mass to be combined, for each object the density profile is normalised to its central density ρ s and the distance from its centre to the quantum Jeans length λ J (ρ s ) corresponding to its central density ρ s . As it is clear from Section 3.1, in terms of Figure 5: Left: The spherically averaged density profile of collapsed objects as a function of the distance from their centre (blue line). The density is plotted relative to that at the object's centre ρ s , and the distance is normalised to the (physical) quantum Jeans length λ J (ρ s ) = 2π/k J (ρ s ). We average over all objects, with barely visible statistical uncertainties. The profile in the inner region closely matches that of a soliton (dashed black line, see eq. (26)). Right: The time evolution of the masses M of the solitons, normalised to the parametric expectation , where a i is the time when the soliton forms and λ J (a i ) = 2π/k J (ρ(a i )). The solitons are binned based on a i , with statistical error bars. Once produced the soliton masses are approximately constant. Moreover, on these variables the soliton density profile is χ 2 1 (x/λ J (ρ s )) and is independent of the soliton mass. Evidently the collapsed objects have a profile that is remarkably close to the soliton, out to a distance λ J (ρ s )/2 1.5R. This confirms our expectation that the objects formed around MRE are supported by quantum pressure. Spatial angular momentum initially carried by the dark matter that later becomes a soliton could be lost during formation or transferred to the fuzzy halo that surrounds the solitons. 80% of the total mass in the pure soliton solution is within λ J (ρ s )/2, so to a good approximation we can identify the mass in the soliton-like part of the collapsed objects as being the same as the total mass of the vacuum soliton. That gravitational interactions result in the dark photon no longer being purely longitudinal and indeed the vector field in the soliton is not longitudinal.
As a further check of the nature of the solitons, we have evaluated the variance of ψ i ( x)/| ψ| 2 over the soliton cores (defined as the region in which the spherically averaged density exceeds ρ s /2). This would represent the spatial variation of the (would-be) constant unit vector u i in eq. (26). As expected, the variance is tiny ( 0.02) for almost all of the objects that form, confirming that the solutions are indeed close to the pure soliton solution. We also calculated the intrinsic angular momentum of the solitons, defined in eq. (30), in units of N . We find that all the components of S/N in the soliton core have a flat probability distribution between −1 and 1.
We note that the initial perturbative linear growth of the k k modes matches the analytic expectation (we analyse this in Appendix D). However, due to the limited simulation time available only few such modes have become nonperturbative by the final time. The compact halos resulting from their collapse are still almost non-existent and contain only a very negligible fraction of DM (see also Figure 9 in the next Section), so we do not attempt to study them in these small scale simulations.
Finally, as mentioned, initially the solitons are produced with quasinormal modes. Although the quasinormal modes are expected to have disappeared by today, they are long-lived with effective Qfactor > 10 3 [59] and could have interesting, possibly observable, consequences, which would be worth exploring in the future. In particular, a substantial fraction of the DM is in spherical waves emitted due to the quasinormal modes (since quasinormal modes initially store an order one fraction of the soliton energy density). Their wavelength is of the physical size of the solitons when they are first emitted.

The soliton mass distribution
To and As we will see in more detail shortly, the solitons produced through the evolution have masses approximately in the range M (0.05 ÷ 0.5)M J (a eq ) = (2.6 · 10 −24 ÷ 2.6 · 10 −23 )M ( eV/m) 3/2 . Ultimately, the solitons have masses inversely proportional to m 3/2 , because the total mass initially contained in a region of volume k −3 has this dependence.
Classically the soliton ground state is stable. Quantum mechanically there is the possibility of decay to gravitons, but this process is exponentially suppressed [60]. Therefore, unless destroyed by e.g. tidal disruption (studied in Section 5.1), they constitute an irreducible component of the DM abundance. One of the most interesting quantities is the soliton mass distribution, which counts the fraction of dark matter in solitons per unit log mass (in plots we use the base 10 logarithm so that f s can easily be estimated). In eq. (36), f s (a, M ) is the fraction of dark matter in solitons with mass less than M , and can be related to the number density of solitons n(a, M ) with mass less than M at time a. 25 The soliton mass distribution in eq. (36) can be evaluated in numerical simulations, at different times. The result is shown in Figure 6. As expected, the distribution is time-dependent, and has a break at low masses that tracks M = M (a) ≡ c M M J (a) (indicated with ticks above the lower axis), because only a  Figure we show the density of the solitons ρ s at their core, which is related to their mass via eq. (29). ρ s is independent of m and parametrically coincides with the average DM density at the time of formation. In particular, it corresponds to the density of the Universe at MRE for the most massive solitons, and is many orders of magnitude larger than the local DM energy density ρ local today.
Unfortunately the limited time range of simulations means we cannot capture the solitons that are produced at a/a eq 7 (evidently, the soliton mass distribution is still evolving at small M in Figure 6). Given that no additional heavy solitons will be subsequently produced, Figure 6 gives a lower bound on f s , and additional production will only strengthen e.g. direct detection signals. We will also see in Section 5.1 that the densest solitons, for which we do have a reliable prediction, are most likely to survive to the present day in the Milky Way.
Interestingly, despite the complicated dynamics, the shape of the soliton mass distribution can be understood via a simple analytic argument based on the initial power spectrum P δ of Figure 2. This gives theoretical control of the soliton mass distribution, and also allows us to estimate the extrapolation of the numerical results in Figure 6 to smaller soliton masses. Estimating the mass distribution requires two inputs: 1) the mass and 2) the number density of the solitons produced, as a function of time a. To fix 1) we crudely approximate that the solitons produced at every time have a unique mass M = M (a) given by eq. (23), with c M 0.45. For 2), we assume that as soon as the comoving quantum Jeans scale drops below the size of an order one fluctuation, this will collapse into a soliton. Therefore, the number of solitons with masses M ÷ M + dM that are produced is expected to be proportional to the 'frequency' with which there are corresponding fluctuations in the initial conditions that are larger than some critical value δ c of order one. To estimate this frequency, we make the crude assumption that the dark photon density field is Gaussian, with power spectrum of Figure 2 (right), although in reality this is not the case.
In Appendix E we derive the resulting soliton mass function.
Despite involving rough approximations, in Figure 6 we see that our analytic argument reproduces the data at large masses, where this has already converged to its late-time value, remarkably well. 26 The analytic prediction does not account for the decrease in the soliton production rate due to DM becoming bound in compact halos at larger scales, which becomes important around a/a eq 30, corresponding to z 100. We therefore indicate on the plot the solitons that are produced before this time, for which the prediction applies.
We note that the argument we used is very similar to that usually employed to estimate the abundance of primordial black holes from small scale curvature perturbations that could be produced during inflation (see, e.g. [61]), with the quantum Jeans scale functioning as an effective 'horizon', in the sense that both prevent the collapse of the order one fluctuations, until they cross the size of that perturbation.
Although we have focused on solitons produced from initially large density fluctuations, we note that additional solitons might form later in compact halos. This could happen directly when the compact halo forms, although in Section 4 we will see that only a small fraction of the mass in the halos will end up in a soliton this way. 27 Alternatively it could occur later, on longer timescales, by gravitational relaxation, see [62][63][64]. Additionally, the solitons already present might increase their mass by accreting the background DM via gravitational relaxation, or solitons might merge together if they become bound into compact halos. This lead to an uncertainty on their mass. We do not try to study these potentially important issues in our present paper.

Fuzzy halos around solitons
Although the bound objects seen in simulations resemble the pure soliton solution of eq. (26) at small distances, as anticipated in Section 3.1 their profile deviates at distances larger than the soliton half-mass radius. Indeed, soliton cores are surrounded by a halo. We dub this a 'fuzzy' halo, because the quantum Jeans length associated with its typical density is only marginally smaller than the size of the halo and the wave-like properties of the vector are still relevant. Indeed, the fuzzy halo is a (time-dependent) solution of eq. (17) where the gravitational potential is locally balanced by a combination of both the quantum pressure and velocity term, in particular v i = 0. Such a halo can be seen surrounding the soliton in Figure 4. Fluctuations on distances of order the de Broglie wavelength can be seen.
In Figure 7 we plot the density profile averaged over relatively heavy solitons, with mass M 26 In this we have set the unfixed parameter in our analysis δc = 0.22 to reproduce the soliton production rate measured in simulations. 27 During collapse, the density in the compact halo increases, usually by a factor 200. If the fluctuation corresponds to spatial scales only slightly larger than k −1 this could be enough for quantum pressure to become relevant and most of the mass might end up in a soliton with a relatively large mass might. Figure 7: The (spherically averaged) density profile around solitons relative to the density at their core ρ s , as a function of the distance from the centre in units of the quantum Jeans scale at the core λ J (ρ s ). Results are shown averaging over all heavy solitons with M > 0.3M J (a eq ), which form relatively early. The density profile matches that of the pure soliton solution for r/λ J (ρ s ) 1/2, and the 'fuzzy' halo outside is well fit by an NFW profile, see eq. (37). The profiles are cut off at a radius where the dark matter density drops below 20 times the mean dark matter density. 0.3M J (a eq ). As in Figure 5, we combine different objects' profiles by plotting their density relative their core density ρ s and measuring the radius in units of λ J (ρ s ). The profile at a given time is cut off at the radius where the density drops below 20 times the mean DM density.
The fuzzy halo is evident as a deviation from the soliton profile for r λ J (ρ s )/2. The scaling property of the soliton solution appears to also apply to their halos, and the profiles from different objects take a universal form (this suggests the existence of a soliton-fuzzy halo relation). At later times the halo extends further from the soliton core, with the region closer to the soliton remaining with a fixed density. The inner part of the halo might form when the overdensity collapses. However, at a/a eq 7 the fuzzy halo extends almost down to the mean DM density at a/a eq 3, so the external parts of the halo are most likely related to the accretion, as the background DM is attracted to the soliton.
The fuzzy halo turn out to be well described by an NFW profile where ρ 0 and r 0 are dimensionless parameters that are universal for all the halos. Fitting the density profile at a/a eq = 7 in the interval 1 ≤ r/λ J (ρ s ) ≤ 7 (where the lower limit is due to the profile transitioning to soliton form) we obtain that ρ 0 0.042 and r 0 1.56 accurately reproduce the data. 28 We do not know how far out the fuzzy halos will grow at times beyond the reach of simulations. In Section 5.1 we will see that the outer parts of the fuzzy halos are destroyed in the late Universe and the parts that are most likely to survive are mostly already formed at a/a eq 7, the final simulated time.
Eq. (37), together with the soliton mass distribution and an input about how far the fuzzy halos extend, allows the distribution of fuzzy halos to be calculated. The fuzzy halos surrounding solitons contain much more DM mass than the solitons themselves. As an indication, the part of a fuzzy halo within 8λ J (ρ s ) (which is a typical distance out to which a fuzzy halo is likely to survive disruption in the late Universe) has mass 6M where M is the mass of the central soliton. The corresponding mass and size in physical units can be read off from eqs. (35) and (49).
Finally, we note that Figure 4 shows that there are overdense filaments connecting solitons analogous to the standard cosmic web, which forms much later at much larger scales. However, these are much less dense than the fuzzy halos and are probably destroyed in the subsequent evolution.

Compact Halos and Primordial Structure Formation
In this Section we focus on scales larger than λ = 2π/k and reconstruct the evolution of the modes in the k 3 part of the spectrum in Figure 2 (right). As we will see, as they become nonperturbative, they induce the formation of a chain of heavier and heavier compact halos: a 'primordial' structure formation. This happens before (and is normally not present in) canonical structure formation, because of the additional small-scale inhomogeneities.
These dynamics are similar to the formation of compact halos in the case of post-inflationary axions, which also has a density power spectrum with a k 3 dependence in the IR. 29 The process of compact halo formation in this case has been studied extensively, both analytically [19-21, 25, 28, 68-70] and numerically [19,71,72]. Additionally [15,28] studied compact halos in vector dark matter produced by inflationary fluctuations using the Press-Schechter approach, similar to our analytic analysis. 30

Dynamics of the IR modes
As shown in Section 3.1, after MRE, modes with k < k (both the k 3 and adiabatic modes) are not affected by quantum pressure and grow linearly. Once a k 3 mode becomes of order one, the linear approximation to eqs. (14) and (15) is no longer applicable. Qualitatively, we expect that at this point the DM density contained in that perturbation collapses into a gravitationally bound object -a 'compact' halo -incorporating already bound objects inside the region, including solitons. Since perturbations on larger and larger scales become nonlinear at later and later times (given that P δ (k) ∝ k 3 ), heavier and heavier compact halos will progressively form. When the adiabatic modes become nonperturbative, at around z 15, they trigger canonical structure formation. We expect that after this time few new compact halos form, and most of those already present eventually become small subhalos of much larger galactic halos, as in Figure 1. 31 Standard structure formation occurs as in cold dark matter cosmology, with the only difference that some of the DM is bound in compact halos.
The small-scale substructure of the compact halos (e.g. any solitons they contain) is clearly affected by quantum pressure and the same may be true very close to the centre of the halos. However, given that they are generated from scales larger than the quantum Jeans scale, if we focus on their properties at scales sufficiently larger than λ (effectively smoothing out small scales), the effect of quantum pressure is expected to be mostly irrelevant. Moreover, at large enough scales the (effective) initial velocity of 29 With a peak of roughly P δ 1 at a scale approximately set by k [65][66][67], although there are large uncertainties related to the decay of the string-domain wall system. 30 In some theories fluctuations in the inflaton can lead to similar dynamics [73,74]. 31 Adiabatic modes on small spatial scales collapse slightly before those on larger scales because adiabatic fluctuations grow logarithmically during radiation domination once they have re-entered the horizon. the field is negligible (because the IR modes do not oscillate, being unaffected by quantum pressure, and the field is highly nonrelativistic). This means that at these scales the equations of motion in eqs. (16), (17), (18) reduce to those of a single component perfect fluid (i.e. with Φ Q = 0) subject only to gravitational interactions, with density ρ = ρ 1 + ρ 2 + ρ 3 , and v v i ( 0 at a < a eq ), i.e.
Despite still being nonlinear, the dynamics of the system at these large scales is simpler and can be analysed by combining standard analytic and numerical approaches. On the analytical side, the so-called Press-Schechter (PS) method handles the evolution of eqs. (38), (39), (40) by determining the number of halos present at every time, based on the power spectrum in the initial conditions. Although this is a model rather than a first principles calculation, it has been shown to capture the main qualitative features, and provide a reasonable quantitative prediction of halo mass functions in many settings [75]. The same equations can be also investigated via N-body simulations. By evolving a system of discrete particles interacting only gravitationally, these reproduce the dynamics of a perfect fluid. Owing to the discretisation into particles, N-body simulations do not lose resolution of collapsed objects in the way that direct SP simulations do, and are therefore better suited to studying the successive chain of structure formation.

Formation of compact halos
In the following we treat the k 3 modes separately from the adiabatic modes, which will be accurate until the adiabatic fluctuations start becoming non-linear, and determine the abundance of compact halos. A halo is a set of gravitationally bound matter. In the remainder of this Section we will not consider the substructure of compact halos in terms of subhalos or solitons. As we will discuss in Section 5.1 we expect that many of the solitons survive intact inside compact halos, although they could also be destroyed, merge with each other or increase in mass by accretion. Compact halos will subsequently be bound inside adiabatic halos. We can predict the distribution of compact halo masses using the Press-Schechter approach. At distances larger than λ , while the linear approximation is valid the overdensity field grows as δ(t, x) = δ(t t eq , x)D[a] with D[a] = 1 + 3 2 a aeq . We expect that at every time regions of space where the overdensity has exceeded an order one critical value δ c have collapsed into a halo. To assign a mass to these regions, we consider the field smoothed over a distance R, and the mass contained in each of them is M = (4π/3)R 3 ρ. Since these regions are expected to collapse into halos of mass up to M , the Press-Schechter anzats [76] is that the fraction of DM in halos with mass > M equals the probability that the field δ s smoothed over R = (3M/4πρ) 1/3 is larger than δ c . This probability is fully determined by the power spectrum P δ since δ is initially Gaussian at scales larger than λ , and so only by the variance σ 2 s ≡ δ 2 s (t, x) , with P δ at t t eq as in eq. (12). In Appendix G, we show that the resulting fraction of DM bound in compact halos f h then satisfies where ν(M ) ≡ δ c /σ s (M ) and we introduced an extra factor of 2 to address the well known cloud-in-cloud problem, as originally done in [76]. In the last equality we defined M ≡ (4π/3)ρa 3 (2π/k ) 3 6.9M eq (see eqs. (22) and (23)) and approximated the spectrum with a single power law k 3 , considering only modes with k < k since modes with k k are affected by quantum pressure, and form solitons. As evident from eq. (41), the halo mass distribution is peaked at with an exponential cutoff at higher masses and a power law suppression (∝ M 1/2 ) at lower masses. Therefore, as anticipated, as time increases the most frequent halos are increasingly heavy. Also as expected, the compact halos are much heavier than the solitons, see eq. (35), and in fact will contain some of the solitons. Eqs. (41) and (42) are reliable only at large enough masses, which come from the largest modes that are least affected by quantum pressure.
To check the validity of this analysis, in Figure 8 we compare the halo mass distribution in eq. (41) using δ c = 1.7 (as predicted by the so-called spherical collapse model [77]) with results of N-body simulations, at different times, i.e. redshifts. The simulation starts at a/a eq = 0.01 1 with initial conditions given by the full (non-Gaussian) density field in eq. (11) -converted to a particle configuration -and vanishing initial velocity v 0. These can be run until the most IR modes in the finite box start to become non-linear, around z = 30, when finite volume systematic errors become significant. 32 From the plot, it is evident that the PS method reproduces the peak and the high mass cutoff of the halo mass distribution at all times. Although the height of the peak is slightly different between the two, the overall agreement is impressive.
We can also use N-body simulations to estimate how large a mass a compact halo must have in order for it to be unaffected by the dynamics of the most UV modes, near the k peak, which are influenced by quantum pressure. To do so, in Figure 22 (left) of Appendix E we compare the halo mass distribution from N-body simulation with initial conditions given by a Gaussian field with the power spectrum of Figure 12, but with a UV cutoff at k > 0.5k (for this momentum range, P δ 0.05, so the initial fluctuations on such scales are indeed very close to Gaussian, and for them the quantum pressure is irrelevant). The two distributions coincide for M 50M eq J . This is only a rough test because, even with initial conditions with a UV cut, once they become non-perturbative the IR modes will source higher k modes, which will still, incorrectly, evolve as in the absence of quantum pressure. Nevertheless, the difference between the two data points in Figure 22 (left) gives a feeling of the impact of high k modes and suggests that eq. (41) can indeed be trusted for large enough M .
The evolution of compact halos according to eq. (41) continues up until z 10÷15, when the adiabatic modes also become nonperturbative. 33 Unfortunately our computational resources do not allow us to simulate the evolution of the adiabatic part and the k 3 modes at the same time, therefore we cannot directly investigate how compact halos are incorporated into the much larger adiabatic halos. We estimate the distribution of compact halos today (if not disrupted) by evaluating df h /d log M at z 10 ÷ 15, when most of the DM is bound in adiabatic halos. To improve accuracy, instead of using the PS analysis we fit the N-body simulation data with a Sheth-Tormen (ST) inspired form [79], which provides a better fit than eq. (41) (see Appendix G for details) to the data points in Figure 8. We indicate the mass distribution reconstructed at z = 10, 15 from the ST fit with gray lines in Figure 8. The peak occurs at around Finally, in Figure 9 we plot the fraction f i (a) ≡ dM df i (a,M ) dM of dark matter collapsed in soliton i = s (from the analytic estimate in eq. (67)), in compact halos i = h and in adiabatic halos i = a, as a function of time. The fraction of DM bound in compact halos is obtained from the ST analysis in Appendix F (with data points from simulations also plotted) excluding objects of mass M < 50M eq J which will be strongly affected by quantum pressure and not correctly described by our analysis. The fraction in adiabatic halos is estimated from a PS analysis, extrapolating the observed power spectrum to smaller scales. 34

Compact halo profiles
We can also study the density profiles of the compact halos. To do so we use an argument normally employed to reconstruct the density of adiabatic halos produced during canonical structure formation (e.g. [80]), which has been used in the context of axion miniclusters by [25,69,70]. This is based on three simple assumptions. 33 After they reenter the horizon, the adiabatic perturbations grow logarithmically in radiation domination because of the gravitational potential generated by the photon perturbation, see e.g. [78] (this does not happen for the isocurvature ones). The logarithmic increase is the largest for the highest adiabatic modes (because they have been subhorizon the longest), and this changes the slope of the adiabatic spectrum in Figure 2, making it larger than 0. Standard structure formation therefore proceeds similarly to the 'primordial' one induced when the k 3 modes become non-perturbative, but with a power spectrum that is much flatter and is initially much smaller. 34 In particular, we put a UV cut the power spectrum at the point where the isocurvature spectrum dominates, which depends on m. In Figure 9 the width of the adiabatic collapse fraction line corresponds to the, barely noticeable, effect of varying m from 10 −5 eV to 1 eV.  Figure 9: The fraction of dark matter collapsed in solitons (blue), compact halos (green) and adiabatic halos (purple) as a function of redshift. The fraction of mass in solitons is calculated from the analytic analysis of Section 3.6 as in Figure 6. The fraction of DM bound in compact halos is evaluated from N-body simulations (data points) and from the analytic Sheth-Tormen analysis (line), with parameters fit to the data. The fraction collapsed in adiabatic perturbations is calculated using the Press-Schechter method from an extrapolation of observed power spectrum. This plot ignores the fact that the dynamics of objects on larger distance scales might halt the formation of structure on smaller scales: Soliton formation is prevented after z 100 when most of the DM is bound into compact halos (indicated by the predicted soliton fraction being plotted dashed for z < 100); similarly, compact halos are expected to stop forming at z 15.
(1) We assume that all the compact halos that form at a given time have a particular mass, determined parametrically by peak value PS distribution in eq. (41), i.e. halos with mass M are created at a c (M ) = (8π 3/2 /3 7/4 )(δ c /ν c ) M /M a eq with ν c an order one parameter. Halos are assumed to remain subsequently undisturbed. Clearly this is a rough approximation, since at a given time halos with some range of masses will form, and also because a halo of a certain mass could form from mergers of lighter halos.
(2) Since halos grow slowly from the low-density fluctuations (as opposed to from the collapse of already large scale fluctuations), it is reasonable to expect that they all have a fixed overdensity with respect to the background density at their moment of formation. We will therefore assume that the mean density of the halos, which we call ρ M , is ∆ times larger than the average DM density at the time of formation, i.e. ρ M = ∆ ρ(a c (M )), with ∆ > 1 a universal time-independent number for all the halos.
As we will see from N-body simulations, the halo profiles are well described by the NFW form [81,82] with scale radius r 0 and density ρ 0 . As is well known, the total mass in such a halo 4π ∞ 0 dkk 2 ρ(r) is logarithmically divergent. The halo extends up to an edge defined by the virial radius R ∆ , related to the scale radius r 0 via the so-called 'concentration' parameter c ∆ > 1 as R ∆ = c ∆ r 0 (c ∆ parameterises how large the halo is with respect to r 0 ). The average density is therefore ρ M = 4π R ∆ 0 dkk 2 ρ(r)/(4πR 3 ∆ /3). Note that the halo is completely specified by the three parameters {ρ 0 , r 0 , c ∆ }. where all halos within 1% of the labeled mass are averaged. Also plotted is the analytic expectation of the halo profiles given in eqs. (45) and (46). The profiles are plotted in terms of the scale parameter r 0 of the NFW profile (which depends on the halo mass), left, and the physical distance from the center, right. Note that lighter halos, which are produced earlier, have larger density parameter ρ 0 as can be seen in the left panel, despite having smaller density at a fixed physical distance from their centres. In the right panel, we also plot the density ρ J (r), defined by λ J (ρ J ) = r. Quantum effects are relevant for densities ρ ρ J , since λ J ∝ ρ −1/4 .
(3) Given that the formation of halos is self-similar during the evolution, we assume that all the halos are formed with a universal concentration parameter, irrespectively of the time when they are created.
The assumptions above allow the parameters ρ 0 and r 0 of the NFW profiles to be calculated in terms of ν c , ∆ and c ∆ . We report the full analytic formulas in Appendix G. The results for ν c 0.67, c ∆ 30 ÷ 50 and ∆ 300 ÷ 500, which we will see fits the numerical results well, are 35 The lighter halos are denser, since they are created at earlier times when the background density is larger. Note that ρ 0 represents the density of the halo at r r 0 . Since the halo extends a factor c ∆ > 1 beyond r 0 the average density is smaller than suggested by eq. (45) and, for c ∆ 1, is approximately 36 The energy density profile of the compact halos can be calculated in the same N-body simulations discussed in Section 4.2. In Figure 10 we show the (spherically averaged) profile of halos of different masses at z = 70, together with the predicted analytic form in eqs. (44), (45) and (46). Albeit not from first principles, we see that the analytic analysis is in remarkable agreement with the data, which indeed 35 For canonical structure formation, the parameters are νc O(1), c∆ 4, ∆ 200 [83]. 36 Note that the average density of compact halos with mass M > 5000MJ is smaller than the local dark matter density.
follow an NFW profile. Although in Figure 10 we only show the profiles of the halos present at a fixed redshift, we checked that the analytic expectation in eqs. (45) and (46) matches simulation results at other accessible values of the redshift, for approximately the same values of ν c , c ∆ , ∆. Given the rough nature of our analytic analysis (e.g neglecting that some compact halos will merge, be bound inside bigger halos, accrete mass) we do not attempt a global fit, and the values given above are sufficiently accurate for our purposes.
We also note that, similarly to the halo mass distribution, the halo profile should be trusted only for compact halos with mass M 50M eq J , which come from modes k k so are unaffected by quantum pressure (at least during their linear growth). Additionally, the density profiles of compact halos in N-body simulations are only reliable in regions where the quantum pressure is negligible. 37 In Figure 10 (right) we plot the curve ρ J (r) ∝ r −4 defined by λ J (ρ J ) = r. Since λ J ∝ ρ −1/4 , ρ J (r) is the minimum average density a region of size r can have before quantum pressure becomes relevant. 38 Therefore, regions of the halos for which ρ(r) > ρ J (r) are expected to be mostly unaffected by quantum pressure. We see that the majority of the compact halo mass is self-consistently in a region where quantum pressure is negligible, especially for relatively heavy halos, although this fails at the centre of the halos. For instance, for M/M eq J 200 less than 5% of the mass of the NFW profile is inside r c .

Late Time Evolution
Possible observational and experimental signals of the solitons, their fuzzy halos, and the compact halos (which we collectively refer to as dark matter clumps) rely on them surviving undisrupted to the present day. The resulting discovery potential will depend on how often collisions between a clump and e.g. an observer on Earth occurs.

Survival of the clumps
There are various processes that could destroy a clump, the most important of which occur when it is bound inside a larger structure, e.g. a larger clump or the Milky Way itself. As previously discussed, the cosmological history of the clumps is extremely complicated. The majority of solitons (and also the relatively low mass compact halos) will be bound in a series of compact halos of increasing mass and decreasing average density. Subsequently the dark matter clumps will fall into adiabatic halos. For direct detection we are interested in clumps that are eventually bound in the Milky Way, with orbits that cross the neighbourhood of the Earth. Such orbits could have a range of forms (mostly in the disk, on a rosette orbit, etc. [84]), which will affect the survival probability.
Given this complexity we do not attempt to directly track the survival probability of a clump from when it formed to the present day in full detail. Instead we estimate whether typical clumps survive by calculating the rate of the processes that are most likely to lead to their destruction. Most of our analysis will be done treating the clumps as being made up of classical particles, although this will not be accurate for the solitons. Despite these approximations, our analysis will be sufficient to indicate whether or not solitons and compact halos are likely to survive in various environments. On the other hand, further dedicated study would certainly be worthwhile, especially of the possible destruction during hierarchical structure formation, which (although challenging to study) might be important. 37 In Figure 10 we plot the data only down to radius r soft below which the N-body results do not reproduce the equations of motion of even cold dark matter, for numerical implementation reasons. Further details are given in Appendix F. 38 If in a region of size r the density is ρ < ρJ (r), due to quantum pressure the configuration will relax (into a soliton) with radius larger than r (and density even smaller than ρ).
Among the mechanisms that can lead to destruction are tidal forces from a central potential, dynamical friction and tidal shocking by a central potential. For an object in the Milky Way there is also destruction by encounters with stars and by tidal forces from crossing the galactic disk. Destruction by tidal shocks is also possible when a soliton or compact halos falls into a larger compact halo (or a relatively small halo from a small scale adiabatic fluctuation) during hierarchical structure formation. We analyse each of the destruction processes in detail Appendix H. Here we summarise the main conclusions: First, the only property of a clump that the rate of destruction depends on is its mean density ρ(r) = r 0 ρ(r)4πr 2 dr/ 4/3πr 3 (for all the important processes). Consequently the disruption is independent of the dark photon mass. Regions of clumps within which the mean dark matter density ρ(r) is larger than about 0.05 eV 4 are likely to survive to the present day if they end up on a typical orbit that passes through the Earth.
The solitons have central densities ρ s (0.1 ÷ 100) eV 4 (see Figure 6). From Figure 7, the density profile switches from soliton-like to the fuzzy halo NFW form at radius of r edge 2λ J (ρ s )/3. At this radius ρ(r edge ) 0.2ρ s (with local density ρ(r edge ) 0.05ρ s ). Consequently, the majority of the solitons, and especially the relatively heavy ones, which form at the times accessible in the simulations in Section 3, are likely to survive.
Given the profiles of the fuzzy halos around solitons, ρ 0.05 eV 4 corresponds to regions of these with local densities ρ(r) 10 −3 eV 4 . Therefore, the outer part of the fuzzy halos around solitons is likely to be destroyed. This is to be expected given that the halos extend out to densities barely greater than the dark matter density in the neighbourhood of the Earth. However, the central part of the fuzzy halos are likely to survive. The density in the part of the fuzzy halos that survives corresponds to an enhancement over the local dark matter density (which for definiteness we fix to 0.5 GeV/ cm 3 ) of ∼ 10 3 .
The compact halos are far less dense than the solitons and much more likely to be destroyed. Compact halos with masses in the range 10 2 ÷ 10 4 M eq J (see Figure 17) have typical mean densities ρ (10 −6 ÷ 10 −3 ) eV 4 and the local densities at their edges are 10 −7 ÷ 10 −4 eV 4 (see Figure 23 right), which at the lower end is even smaller than the local dark matter density. Consequently, the majority of such halos, which contain 70% of the dark matter at z 10, are likely to be destroyed. The relatively dense core of the compact halo might survive, but these will contain a smaller fraction of the dark matter. Additionally, less massive (more dense) compact halos that are subsequently bound in the most massive compact halos might survive.

Collision rate
We now analyse the rate at which a point-like observer collides with a compact object. Given the uncertainty about whether the compact halos survive to the present day, we focus on collisions with the solitons and the parts of their surrounding fuzzy halos with local density ρ 0.01 eV 4 (corresponding to a mean density ρ 0.05 eV 4 ), which are likely to persist to today. With an eye to future analysis of direct detection signals, we present our results assuming a local dark matter density ρ local = 0.5 GeV/ cm 3 (this is subject to significant uncertainties, see e.g. [85]), corresponding to the Earth's local environment, but they can easily be rescaled to other densities.
For convenience, we note that a soliton of mass M has central density, given by eq. (29), that corresponds to and a size set by λ J (ρ s ), which, from eq. (20), is λ J (ρ s ) = 4.6 · 10 3 km eV m where M eq J is given in terms of solar masses, as a function of m, in eq. (34). We can easily obtain an analytic estimate of the collision rate. Approximating that the solitons all have a single mass M , their local number density is where, as in Section 3, f s is the fraction of DM in solitons. 39 The number of collisions per unit time between a point in space (e.g. a dark matter detector) and solitons is therefore where R is the maximum radius that the soliton profile extends to, which in the second line we have set to λ J (ρ s ) given the results of Section 3.7.
Apart from the explicit dependence, all the factors in the second line of eq. (51) are independent of the dark photon mass m. Consequently the interaction rate is larger when the dark photon mass is larger. This has a straightforward interpretation: the fraction of DM in solitons, f s , and the densities of the solitons are independent of m; therefore the fraction of time that a point spends inside a soliton is independent of m. However, the physical size of the solitons (with fixed M/M eq J ) is inversely proportional to m 1/2 (see eq. (49)), so each encounter with a soliton lasts less time when m is larger. In particular, a collision with impact parameter b with a soliton lasts for roughly Consequently the collision lasts roughly a minute for m eV. From Figure 6 we see that during the collision the dark matter density is typically enhanced by a factor ∼ 10 4 ÷ 10 7 compared to the local density. Given that the halos around solitons might survive out to approximately 10λ J (ρ s ), the interaction rate with these is expected to be substantially larger and the interaction time significantly longer (although the typical enhancement over the local density will be smaller).
We can make more accurate predictions using the simulation and analytical results for the soliton mass function from Section 3.6. We continue to consider the rate at which a single point collides with a clump. Given the size of the solitons eq. (49) this will be appropriate for direct detection signals, and also e.g. neutron stars provided m 10 5 eV. It is straightforward to repeat our calculations to determine the rate of collision between an object with size comparable to λ J (ρ s ) and clumps.
In Figure 11 we plot the rate Γ at which collisions that result in a dark matter density enhancement of at least ρ/ρ local occur. In other words, a point is expected to experience a collision that results in a density enhancement (at its peak) of at least ρ/ρ local roughly once per Γ −1 time. In this plot we assume ρ local = 0.5 GeV/ cm 3 , and a mean relative velocity of 10 −3 between the solitons/fuzzy halos and the clump. The scaling with Γ ∝ m 1/2 in the estimate of eq. (51) is exact also for the full analysis, so we factor this out on the vertical axis.  Figure 11: The expected rate Γ at which a point-like observer collides with a dark photon clump in such a way that the (maximum) resulting dark matter density enhancement seen is ≥ ρ/ρ. Results are shown considering collisions with the central regions of solitons, r/λ J (ρ s ) < 0.6, (blue) and due to colliding with either a soliton or its surrounding halo (purple). As discussed in the main text, Γ is proportional to m 1/2 . Results are plotted using the soliton mass function at z = 500 obtained directly from simulation data (solid) and from the analytic extrapolation of the soliton mass function to z = 100 when the collapse of the k 3 fluctuations is expected to prevent further formation of solitons (dashed). The fuzzy halos around the solitons are assumed to have the NFW form seen in simulations. The local dark matter density is fixed to 0.5 GeV/cm 3 and the corresponding physical densities that the observer passes through are shown on the upper axis. Regions of the halos surrounding solitons with densities 10 −2 eV are likely to be disrupted in the local environment. The values of ρ/ρ local where this will affect Γ are plotted in grey. The majority of the solitons contributing to the rate are expected to survive to the present day. Figure 11 shows Γ considering only collisions with the solitons (given Figure 7, we define this as the region r < λ J (ρ s )) and allowing collisions with the halos surrounding the solitons. Additionally we plot results obtained from the soliton mass function at the final simulation time and from the analytic extrapolation, see Figure 6. The rate of collision with the halos surrounding the solitons are obtained by assuming these take the form of the NFW halo plotted in Figure 7. 40 We also indicate the values of ρ/ρ local such that some of the halos that would contribute to the rate are likely to have been destroyed. For a dark photon mass of m eV collisions with solitons leading to enhancements in the dark matter density of a factor of 10 3 ÷ 10 4 occur on reasonable experimental timescales, and collisions with the surrounding halos are even more frequent.
As well as the maximum dark photon density experienced, the time that an enhancement lasts for Figure 12: The dark matter density as a function of time encountered by a point-like observer colliding with a soliton, and its fuzzy halo. The soliton's central density is fixed to ρ s = 10 eV 4 . We show results for impact parameter b varying from 0 (passing through the centre of the soliton) to those that only pass through the fuzzy halo, measured in units of λ J (ρ s ), the quantum Jeans distance at the soliton core (the soliton-like part of the clump extends to roughly λ J (ρ s )). We also indicate the parts of the fuzzy halo that are likely to have been disrupted prior to the collision, in the Earth's local environment. and the density profile experienced is also potentially important for detection and observation. This can easily be extracted from the density profile of the solitons/fuzzy halos (and likewise the compact halos if they are assumed to survive). In Figure 12 we plot an example of the density seen as a function of time during a collision with a soliton/fuzzy halo of central density 10 eV 4 (corresponding to soliton mass of 0.16M eq J , in the main part of the soliton mass function of Figure 6), for different impact parameters that both lead to collisions with the soliton-like part of the object and that are just with the fuzzy halo. We also indicate the part of the signal corresponding to regions of the fuzzy halo that are likely to have been disrupted prior to the collision. Such disruption could lead to more diffuse streams of dark matter around the object [87], resulting in interesting features in the density as a function of time when an observer passes through the edge of a dark photon clump, although we do not investigate this further in our current work.
An additional potentially important feature is that the dark photon field is coherent over the entire soliton core (since the field's de Broglie wavelength inside the core coincides with the size of the soliton). This is far larger than for unbound dark matter, for which it is 2 · 10 −4 m ( eV/m) 10 −3 /v . Even in the fuzzy halo the typical velocity (from the NFW profile) is tiny, and the de Broglie wavelength is typically less than a factor of 10 smaller than in the soliton core, typically in the range 10 2 km( eV/m) 1/2 (see eq. (49)), still many orders of magnitude larger than it would otherwise be. Finally, we note that we have assumed that the soliton retains its form throughout the collision. In reality, the tidal forces experienced by a clump during a collision with the Earth are likely to lead to its eventual destruction. However, the clump's profile remains mostly unchanged until long after passing the Earth and this effect will not alter the density profile seen by a detector [27]. 41

Discussion and Future Directions
There are a plethora of potential observational and experimental implications of the solitons and their fuzzy halos, which we will study in detail in a companion paper. For m eV, a typical observer in the local environment collides with solitons regularly, and during each collision the dark matter density is boosted by a factor of up to 10 6 and is far more coherent than it is outside solitons. For smaller m, collisions are rarer, but the clumps contain more mass and are larger in size. Although less relevant for direct detection, there is a range of possible striking astrophysical and indirect signals in this case. If the vector has self-interactions, the solitons might be unstable and could explode, leading to further potential signals, similar to the analysis for axion stars in [89]. Moreover, compact halos have a macroscopic mass (of the order of the mass of a planet for m 10 −5 eV), and -being abundant in our Galaxy -can in principle be observed and constrained via gravitational measurements.
Notice that the classical description of the field is valid when the number of particles N within a de Broglie wavelength λ dB ≡ 2π/(mv) (of the order of the inverse typical momentum) is much greater than 1: i.e. N = ρλ 3 dB /m 1. As is well known, for this to be valid today requires m 10 eV if the field is made out of waves. However, inside the solitons the occupation number is much larger, and thereforeas we now show -the validity of our results extends to much larger m.
The crucial condition for the solitons to form is that the vector bosons are classical at around MRE when the overdensities start collapsing, as described in Section 3. Since the soliton energy density and occupation number remain constant, solitons never exit the classical regime afterwards. Classicality at MRE requires This condition is parametrically the same as requiring that the occupation number in the core of the soliton is large (ρ s R 3 /m M/m 1), since the density and size of the solitons match the typical density of the field and the size of the fluctuations at MRE. Eq. (53) implies that our results will apply for all the masses in the range m 10 8 GeV. The upper limit happens to coincide with the largest m that does not exceed the Hubble parameter during inflation, requiring that the dark photon makes up the full relic abundance, see eq. (1).
For m 10 eV in most of the Universe one needs to consider the heavy gauge bosons as quantum objects (similarly to WIMPS), however the solitons are perfectly captured by the classical description. 42 In this regime, the soliton radius is tiny, but the number density of solitons and their encounter rate with the Earth (and other astrophysical objects) is huge. There are potentially dramatic signatures, e.g. at direct detection experiments that aim to detect the collision of single particles.

Possible improvements and extensions
There are a number of ways that our analysis could be refined or extended. One direction is to better test, and if necessary improve, our analytic prediction of the soliton mass function from the density power spectrum. 43 This would be particularly valuable when studying similar dynamics in a theory where the initial power spectrum changes as the theory's parameters vary, or in the context of a theory where the focusing of ultralight DM, studied in [88] if the field is made of waves. 42 We expect the solitons to be stable, even if surrounded by a non-classical background (indeed they might accrete via gravitational relaxation), although a detailed analysis would be useful. 43 Given the local (quadratic) non-Gaussianities in the initial energy density field and the similarity with primordial black holes discussed in Section 3.5, it could be useful to follow approaches already developed in this case [90]. initial power spectrum is uncertain. It would also be interesting to study the fuzzy halos that surround the solitons in detail. For example, simulations of ultra-light axions seem to show a relationship between the masses of solitons that form at the centre of galaxies in such models and their halos [91], and it would be nice to understand if our fuzzy halos follow a similar relationship. Additionally, it would be useful to study the transition between collapse to objects supported by quantum pressure and collapse to objects well approximated as being composed of cold DM more systematically. Although we have analysed the compact halos assuming that quantum pressure is negligible, which is self-consistent over most of their profile especially for the most massive objects, solitons are likely to still form at their cores, and these could be relevant for observation or detection signals. It would also be interesting to study the effect of the vector having self-interactions, or interactions with the SM or new fields. For example, this could lead to the solitons decaying on cosmological timescales, potentially leading to observational signals [92][93][94][95], and it could also affect the rate at which solitons form or gain mass by dynamical relaxation [96].
A key direction for future work is to determine the probability that the solitons, fuzzy halos and compact halos survive to the present-day more accurately. In our analysis we had to make multiple approximations and focus only on particular processes. It would certainly be valuable to better understand the probability of destruction or possible accretion during hierarchical structure formation and also to determine whether wave-like effects alter the probability of the solitons or fuzzy halos surviving. One could also carry out a more detailed study of the probability of destruction by collisions with stars, e.g. along the lines of [69,97,98]. The solitons, as mentioned in Section 3.5, could also merge in the late Universe or when they are bounded into halos, which could change their mass distribution [99].
On the numerical side, there are several directions in which our approach could be developed, which could help address some of the issues above. Our simulations of the SP system could only reach a/a eq = 7 and it was impossible for us to use them when analysing the compact halos in Section 4. An adaptive mesh approach, e.g. as employed in [100,101] could enable a much greater range in a/a eq . One could also employ a hybrid SP -N-body approach, evolving the SP equations only in regions where quantum pressure is relevant, or modify the equations of motion in N-body simulations to attempt to reproduce the effects of quantum pressure. With such work, one might be able to see whether the solitons survive if they are bound in a compact halo.
It would be valuable to understand whether vector solitons could form with other production mechanisms, e.g. by gravitational condensation analogous to the way that QCD axion stars are thought to form in the post-inflationary scenario [62,63]. We have focused on a very minimal theory, but it would also be interesting to study the impact of additional interactions. As mentioned at the end of Section 2, these could lead to changes at all stages, from the production of inflationary fluctuations through to the dynamics around H , MRE and the present day (see also [52]). Additionally, there could be new types of compact object that form when the vector has interactions [102].

A Details of the Initial Conditions from Inflation
In this Appendix we discuss some more details on the production of the vector field during inflation, which was first analysed in [10] (see also [45]). In particular, we give the derivation of eq. (8) and report the useful analytic approximation eq. (12).
To derive eq. (8), we estimate the solution A L (t, k) of eq. (7) in the limits k k and k k . While relativistic and superhorizon, for all modes eq. (7) is approximated by (∂ 2 t + H∂ t )A L = 0, with obvious solution A L A L,0 (i.e. the modes are frozen) and ρ m 2 A 2 L /a 2 ∝ a −2 . However: (1) Even after they become nonrelativistic, the modes with k < k while superhorizon still follow (∂ 2 t + H∂ t )A L = 0. Therefore A L A L,0 , and their energy density decreases as ρ m 2 A 2 L /a 2 ∝ a −2 . This behaviour of ρ is crucially different than for a scalar field, for which ρ is frozen for nonrelativistic superhorizon modes. 44 The solution A L A L,0 is valid until H = m (i.e. a = a ), after which eq. (7) is approximated by (∂ 2 t + H∂ t + m 2 )A L = 0, with solution A L ∝ a −1/2 and ρ ∝ a −3 , which is the usual matter behaviour.
(2) On the other hand, the modes with k > k follow (∂ 2 t + 3H∂ t + k 2 /a 2 )A L = 0 when they enter the horizon at k/a = H (while still relativistic), and have solution A L ∝ a −1 so ρ ∝ a −4 . When they become nonrelativistic at k/a = m, they follow (∂ 2 t + H∂ t + m 2 )A L = 0, as before with solution A L ∝ a −1/2 and ρ ∝ a −3 .
Summarising, after they become nonrelativistic and subhorizon, all modes behave like matter. Approximating the transition between the different regimes as immediate, we have for low frequency modes A L /A L,0 (a /a) 1/2 , while for high frequency modes A L /A L,0 (a e /a nr )(a nr /a) 1/2 = (k /k) 3/2 (a /a) 1/2 , where k/a e ≡ H(a e ) and k/a nr ≡ m and we assumed radiation domination, for which H ∝ a −2 . This allows us to conclude that in radiation domination A L has a Gaussian power spectrum given by (using eq. (8)) where the exact form of F A L [x] (see Figure 2 left) can be by extracted by solving numerically eq. (7), but from the discussion above reported in the right hand side of eq. (9). Finally, using the fact that ∂ t A L and A L are independent Gaussian fields (and that P ∂tA L = m 2 P A L ), it is straightforward to get an expression for the density power spectrum [10], 44 As mentioned in the main text, the difference is due to the form of the mass terms of such particles, which controls the energy density of nonrelativistic superhorizon modes: for a scalar, 1 2 m 2 ϕ 2 does not change during the universe expansion, while 1 2 m 2 g ij AiAj ∝ a −2 for a vector.
where the second equality is the approximate expression reported in eq. (12). Note that the power P δ in [10] is a factor of 2 larger, because that reference neglected the first term in eq. (6).

B Evolution of Overdensities
The evolution of small density fluctuations during radiation and matter domination can be analysed perturbatively, with and without quantum pressure. Although this is not applicable to the large density perturbations at the scales k k (see Figure 3), the perturbative analysis still gives a useful hint towards the type of effects quantum pressure might have.
In the absence of quantum pressure, the evolution of a density perturbation δ(k) (in momentum space) as a function of time t, with scale factor a(t), is well known [103] δ where here f = Ω A /Ω M with Ω M the total matter content (this differs from e.g. f s used in the main text, which was normalised to Ω DM ). Eq. (56) is obtained by combining eqs. (16), (17) and (18) (with Φ Q = 0) in the limit of small overdensity and velocity. Defining y = a/a eq , eq. (56) simplifies to If f = 1 this leads to the usual growing solution δ = δ 0 (1 + 3 2 a aeq ), where δ 0 is the overdensity at an early time, deep inside radiation domination: overdensity are frozen during radiation domination, and grow linearly in matter domination. For f = 0.84 the (growing) solution is slightly modified, but the perturbation is still almost completely frozen before MRE and then grows. At large a/a eq we have δ ∝ (a/a eq ) √ 1+24f /4 . Quantum pressure leads to an additional term, and the perturbation evolves according to where k J is the quantum Jeans momentum, and k/k J ∝ y −1/4 . For k k J the solutions oscillate, both in radiation domination (y 1, when they would otherwise be frozen) and in matter domination (y 1, when they would otherwise grow). For k k J the evolution is as in the absence of quantum pressure, as expected. A detailed analysis, including a study of the applicability of the perturbative expansion and the effects of the next terms in the expansion can be found in [104].

C Solving the Schrödinger-Poisson equations C.1 At matter radiation equality
Here we summarise how the Schrödinger-Poisson (SP) system of equations eq. (14) can be written in a form suitable for the numerical evolution and how we implement a realisation of the initial conditions from inflation.
Around MRE the scale factor a(t) satisfiesȧ = aH(a) with where H 2 eq ≡ 8πGρ eq /3, so that (neglecting changes in the number of degrees of freedom, which is appropriate for the Standard Model) Note that the previous equation can be integrated exactly in conformal time. Similarly to [101,105,106], we rewrite the Schrödinger-Poisson system as where t = t/T , where T is any inverse mass scale. We also definedt = dt /a 2 (so that dt/dt = 1/a 2 ).
Given the form of the initial conditions, it is natural to choose the typical length to be L = 2π(ma /a eq ) −1 . As in Section 3.2, in the limit a /a eq 1, a /a eq = (H eq /( With this choice, eq. (60) simply becomes and the dependence on H eq (as well as any numerical input except for β) have dropped out from eqs. (61), (62) and (63). From now on we will set a eq = 1.
Next we discuss how we generate the initial conditions (at a a eq ). The field configuration ψ i is obtained from A i using the definition A i = 1 √ 2m 2 a 3 (ψ i e −imt + c.c.), which implies Re[ψ i ] = (a 3 /2) 1/2 (m cos(mt)A i − sin(mt)Ȧ i ) and Im[ψ i ] = −(a 3 /2) 1/2 (m sin(mt)A i + cos(mt)Ȧ i ). In the initial conditions the field A i andȦ i have only longitudinal component (A L andȦ L ). This is generated, using standard algorithms for a random field, according to the power spectra, P A L and P ∂tA L , defined by eq. (5). We use the expression of P A L in eq. (9). The same expression holds also for P ∂tA L , modulo a factor of m 2 .
The above procedure fixes the initial conditions except for the overall amplitude of ψ i , that depends on the amount of vector dark matter present (ultimately linked to H I ). We normalise it such that i |ψ i | 2 = f ρ eq /2, where as before f = Ω A /Ω M , 45 which we set to f = 0.84, so that the vector is all of the observed dark matter. 46 Consequently the amplitude of ψ is set such that Integrating eq. (62) with the initial configuration of ψ just described gives the corresponding initial configuration of Φ . 45 This is because 1 2Ȧ 2 i + 1 2 m 2 A 2 i = |ψi| 2 at the time when a = qeq = 1. 46 Note that given the power a −3/2 in the definition of ψ from A, i |ψi| 2 is constant in time; moreover, since a = 1 at MRE, i |ψi| 2 is the dark matter energy density at MRE.
Consequently, the SP system at around MRE boils down to eqs. (61) and (62) with a given by integrating eq. (63) with initial amplitude in eq. (64). 47 Note that the only free parameters on which these equations depend are g R and Ω A /Ω M . Additionally, from eqs. (61) and (62) and (63) with a eq = 1 it follows that Φ is of order one (which is the size of ∇ and ∂t when applied to Φ) when a 1/δ from eq. (62), and at this time the gravitational potential becomes relevant in the dynamics and ψ starts evolving nonlinearly (though, as we saw, perturbations do not collapse due to quantum pressure).

C.2 Numerical solution
We solve eqs. (61) and (62) (together with eq. (63)) on a finite lattice of constant comoving lattice spacing and in a periodic box of constant comoving length. We use a 6th order pseudo-spectral algorithm, developed in [108] and applied to study the cosmology of axion stars in [62] and fuzzy dark matter in [109] (see also e.g. [110] for other similar implementations). Given the fields ψ i at time t, the fields at time t + ∆t are obtained by evaluating where the product is ordered from right to left so the α = 1 part is applied to ψ i (t) first, and the evolution with the k 2 operator is understood to happen in momentum space. Here k denotes the momentum in units of 2π/L. The potential Φ is re-calculated after every step with the momentum operator, by solving the Poisson equation, eq. (62), numerically (by transforming to momentum space and back). As discussed in [62], this pseudo-spectral algorithm has several beneficial features and advantages over other approaches. It automatically conserves particle number and it has a high degree of stability (i.e. there are no spurious growing modes). Compared to lower order pseudo-spectral algorithms, such as used in [106], much larger ∆t are allowed without introducing significant systematic uncertainties, so simulations are faster (detailed analysis comparing algorithms of different orders can be found in [109]).

C.3 Systematic uncertainties
For our main data collection runs we use a lattice of N 3 = 1080 3 points, with a comoving box length of = 3.75λ . Evolution is started at a/a eq = 0.01, and stopped at a f /a eq = 7. We vary the timestep throughout a simulation as ∆t = 0.00025a eq /a. 48 We now show that these choices lead to negligible systematic uncertainties, which arise from various sources: 1. Finite lattice spacing. Once formed the solitons have constant physical size, whereas the lattice spacing increases ∝ a. The size of the resulting systematic uncertainties is controlled by the hierarchy of the physical distance between lattice points at a f and the size of the soliton cores. For negligible systematic uncertainties we expect that a f /N λ J (ρ s ) is required, where ρ s is the central density of a soliton. Since λ J is a decreasing function of ρ s , we require that this inequality is satisfied for the densest, i.e. largest mass, solitons that can form. λ J (ρ s ) is parametrically set by λ . Figure 13: Left: Comparison between the spherically averaged density profile around the densest object in a simulation at a/a eq = 14 with the space-step we use for our main runs ∆x 0 and with finer resolutions 0.66∆x 0 and 0.5∆x 0 (starting from identical initial conditions). Although resolution of the soliton core is starting to be lost in the ∆x 0 run, the central density and soliton mass are still accurate to few percent level. Right: The maximum density in a simulation obtained evolved with the space-step used for our main runs ∆x 0 , and with finer resolution 0.66∆x 0 , compared to that obtained with a 0.5∆x 0 (starting from identical initial conditions). At a/a eq our choice of space-step introduces less than a ∼ 10% systematic uncertainty.
To confirm that our choice of parameters leads to negligible systematics from this source, we compare results from a single simulation with our main value of ∆x = /N = ∆x 0 with results starting from identical initial conditions, but with ∆x a factor of 2/3 and 1/2 smaller. To have sufficient computing power to do the finer resolution tests we carry out this test in boxes with smaller . The initial conditions that we choose happen to have a largest soliton mass ∼ 0.17M J (a eq ).
Since k J (ρ s ) ∝ 1/M we run our lattice spacing tests to a f = 14, which gives the same value of a f / (N λ J (ρ s )) as in a large simulation with a soliton of mass 0.34M J (a eq ) at a f = 7, which covers the vast majority of the solitons that form in our main simulations. Results are shown in Figure 13.
One can see that the soliton core is just about resolved with our main value of ∆x 0 . The central density of the soliton, and therefore the inferred soliton mass, is still accurately reproduced to about 5% level, which is sufficient for our purposes.
2. Finite time step. As mentioned, owing to the gravitational potential, the evolution in eq. (65) is exact only in the limit ∆t → 0. To test the importance of systematic uncertainties with our choice of ∆t we compare results when identical initial conditions are evolved with our main choice ∆t 0 and with ∆t reduced by a factors of 2 and 3 smaller. Results are shown in Figure 14, where it can be seen that our choice of ∆t 0 is sufficient for percent level accuracy.
3. The initial time of simulations. Even though density perturbations are not expected to evolve much during radiation domination, there is still a numerical question of small a i /a eq must be so that significant systematic uncertainties are not introduced. An estimate of a suitable a i can be obtained from the f = 1, k J → ∞ perturbative prediction δ(a) = δ(0)(1 + 3 2 a/a eq ), which suggests that starting that a i /a eq = 0.01 is enough for percent level accuracy. To test this, we plot the maximum density in simulations as a function of time for different a i in Figure 15 (left). Because of the non-linear dynamics of the system after MRE small changes in the initial conditions can lead Figure 14: Left: Comparison between the spherically averaged density profile around the densest object in a simulation at a/a eq = 10 with the time-step we use for our main runs ∆t 0 and with a time-step of 0.33∆t 0 (starting from identical initial conditions). The almost perfect agreement indicates that the systematic uncertainties from this source are negligible. Right: The maximum density in a simulation obtained evolved with the time-step used for our main runs ∆t 0 , and with a smaller time-step 0.66∆t 0 , compared to that obtained evolving with a time-step of 0.33∆t 0 (starting from identical initial conditions). At a/a eq = 7 our choice of time-step introduces only much less than ∼ 1% systematic uncertainty compared to the smaller time-step.
to large changes at late times, so for each a i plot results averaged over a set of 5 different initial conditions (identical between different a i ). We see that the results with a i /a eq = 0.01 agree with those with a i /a eq = 0.005 (and even a i /a eq = 0.05 would be sufficient).
4. Finally, the finite box size can lead to systematic uncertainties. These are expected to be negligible as long as the most IR modes in a simulation have power spectrum P δ (k) 1, indicating that they are still perturbative, and dynamics on the scale of the box size are not affecting the evolution on smaller scales. With our choices of box size and a f in our main runs, the most IR modes have P δ (k) 0.5 (see Figure 16, right) at the final time, suggesting that this is not a source of major systematics.

D Further Results from Schrödinger-Poisson Simulations
In this Appendix we collect some additional discussion and results from our numerical simulations of the SP equations, supporting the analysis of Section 3.
First we consider the growth of density perturbations in more detail. In Figure 16 (left) we plot the evolution of the power spectrum P δ at different momenta (averaged over multiple simulations) and compare this to the prediction from Appendix B in the limit that quantum pressure is unimportant, k J → ∞. It can be seen that modes with k = 0.1k eq J evolve basically as they would in the absence of quantum pressure, and the growth of modes with k = 0.25k eq J at MRE is slightly hindered by quantum pressure. In contrast, as expected, quantum pressure prevents modes close to the peak of the initial power spectrum (see Figure 2 right) from growing until around MRE. The same dynamics can be seen in Figure 16 (right), where we plot the density power spectrum at different times (note the solitons are not screened when calculating this). At a given time, the growth of modes with k greater than the quantum Figure 15: The mean maximum density in simulations starting from different initial times. Due to the oscillation of modes k > k J even during radiation domination and the non-linear dynamics during gravitational collapse, simulations with initial conditions starting at different a i can differ significantly. However, averaged over multiple runs, there is no systematic effect from starting at a i /a eq = 0.01 (as we do in our main runs) rather than earlier.
Jeans momentum are suppressed. Since the quantum Jeans momentum is larger at larger density, an overdensity on a scale k −1 can grow provided k −1 < k −1 J (ρ) where ρ is the local dark matter density, which might be greater than the mean density. Indeed, at late times there is a clear UV cutoff at P δ k J (ρ max ), corresponding to the inverse size of the most massive solitons produced.
Next we clarify the choice made in Section 3.5 to determine the soliton masses in simulations from their central densities. In Figure 17 (left) we plot the average soliton mass as a function of time, grouping the solitons by the time when they formed a i , as in Figure 5 (right), but showing results for the masses calculated from the central densities M (ρ s ) and from integrating the density profile. In particular, we evaluate the mass inside the radius r s where the spherically averaged density ρ(r) drops a factor of 10 from the central value, i.e. M r = rs 0 ρ(r)4πr 2 dr, since over this region that profile is (on average) soliton-like, see Figure 7. The mass calculated from integrating the density profile is close to that inferred from the central density, as expected since we know that on average the profile is close to the soliton one. However, at least at early times there are a small differences and the mass measured by integrating the density profile initially increases slightly approaching that predicted from the central density. We interpret this as the soliton density profile taking some time to settle down to its eventual form. Meanwhile, apart from oscillations due to quasinormal modes, the solitons' central densities reach their late time values fast (see also Figure 4 left), so are a good way to infer the solitons' eventual masses.
We also note that the time of soliton formation is ambiguous. For the purposes of Figure 5 (right) and Figure 17 (left), we fix a i by the first time when the object has a central density ≥ 200ρ and the averaged density profile does not drop by more than 25% within one lattice spacing from the centre. Given that the gravitational collapse happens fairly fast (see Figure 4 left), the value of c M we extracted is not particularly sensitive to this choice.
In Section 3.7 we studied the fuzzy halo around solitons by considering the halos around relatively heavy solitons in Figure 7. 49 In Figure 17 (right) we show the analogous plot for medium mass solitons, Figure 16: Left: The average growth of different momentum modes as a function of time. The predicted growth in the absence of quantum pressure, k J → ∞, is also plotted. The growth of modes k/k J (ρ eq ) = 0.1 is only slightly affected, while the higher momentum modes have more suppressed growth. Modes k/k J (ρ eq ) 1 are prevented from growing almost entirely until the quantum Jeans momentum in some region k J (ρ) is comparable to k. Right: The density power spectrum in simulations at different times averaged over multiple realisations of the initial conditions, with statistical uncertainties. The quantum Jeans momenta corresponding to the mean dark matter density k J (ρ) and to the (averaged) maximum density in the simulations k J (ρ max ) are indicated. The IR modes grow unaffected by quantum pressure, following the standard expectation. Quantum pressure delays the growth of modes on smaller length scales. Once solitons form the power spectrum is peaked at scales k J (ρ max ), set by the density at their cores. The peak tracks the typical size of the solitons.
which form approximately between 5 a/a eq 7 (see Figure 6). We further restrict to objects that have a halo extending to at least r/λ J (ρ s ) = (0.95, 1.25, 1.6, 2.4) respectively at the four times. We also plot the NFW profile obtained fitting to the fuzzy halo around the relatively heavy solitons. Not surprisingly, at a fixed a/a eq the fuzzy halos extend less far from the medium mass solitons than from the heavy solitons. However, the halos are growing outwards from the medium mass solitons fast, and are approaching the NFW form of the fuzzy halos around the heavy solitons, with the same halo parameters. This is a nice feature, which makes it reasonable for us to assume that the halos around all solitons will eventually take this form, at least out to some cut-off. We do not know what this cut-off is, however fortunately this is not a major source of uncertainty. The typical central densities of the medium mass solitons is about 10 eV 4 , so from Figure 17 (right) the fuzzy halo is likely to be destroyed outside the radius (10 ÷ 15)λ J (ρ s ). This is approximately the distance that the fuzzy halos around the massive solitons already reach in simulations, so it is reasonable to estimate that fuzzy halos around all solitons will reach this size (which is what we assume when studying the rate of collision with an observer in Section 5.2).
In Figure 18 we plot the fraction of dark matter in solitons and the halos around them (with the edges density ρ ≥ 200ρ. Having identified such an object, we label the entire region out to the radius where the spherically averaged density profile drops to ρ = 20ρ as being part of the same object (and so not considered further in finding other compact objects). We also discard objects for which the spherically averaged density profile drops by more than 25% one lattice space away from the core. This removes a few isolated points, which would otherwise be identified as collapsed objects and are generally on the outskirts of a genuine collapsed object, just outside their cutoff. These choices do not make a substantial difference to our results. of the halos defined as the point where the spherically averaged density drops to 20ρ, as before). We also plot the fraction of the dark photon energy that is in non-longitudinal modes, which, due to gravitational interactions, increases from 0 at the start of the simulation and becomes a substantial fraction around MRE.

E Analytic Prediction of the Soliton Mass Distribution
As discussed in Section 3.6, our analytic estimate of the soliton mass distribution is based on two simple assumptions: 1) The solitons produced at any time have a mass M = M (a) given by eq. (23), with c M 0.45.
2) As soon as the comoving Jeans scale drops below the size of a fluctuation, this -if already of order one -should collapse into a soliton. Therefore, the number of solitons produced is expected to be proportional to the 'frequency' that the amplitude of such fluctuations in the initial conditions is larger than a critical value of order one, that we call δ c .
More precisely phrased: Between a 1 a eq and a, the fraction of dark matter collapsed into solitons is equal to the probability Π δ>δc (a 1 , a) that the field δ s smoothed between k J (a 1 ) and k J (a) is larger than is a window function that vanishes outside k 1 < k < k 2 , which therefore selects the component of the field within these momenta. Using 1), the number density of solitons produced per unit time is therefore simply We also plot the fraction of the dark photon energy that is in non-longitudinal modes, defined by 1 − f = (k.A(k)) 2 / A(k) 2 where . . . denotes the spatial average. The field's energy is initially in purely longitudinal modes, but gravitational interactions that become relevant around matter radiation equality transfer energy into transverse modes (e.g. the soliton solution is not pure longitudinal).
does not increase with D[a] ∝ 1 + (3/2)a/a eq , since -as mentioned -at those scales the overdensities oscillate without growing.
By integrating over a 1 < a < a 2 , a crude estimate for the soliton mass distribution at a = a 2 in eq. (36) is where the last equality is valid for M < M (a 2 ) (for M > M (a 2 ) the mass distribution vanishes). Here a M is the inverse of the function M (a) in eq. (23), and corresponds to the scale factor when solitons with mass M are produced. Note that, aside from the freedom in the choice of the window function W , the only free parameter in eq. (66) is δ c . Unfortunately it is not possible to compute the probability Π δ>δc (a 1 , a) rigorously starting from the power spectrum P δ alone, since, as mentioned, δ(x) it is not Gaussian at k k . We leave this for a future work. Here we limit ourselves to the Gaussian approximation. For a Gaussian distributed field, this probability is fully determined by the field's variance σ 2 δ (a 1 , a) ≡ δ 2 s (x) = dk/kP δ (k)W 2 (k; k J (a 1 ), k J (a)) (see eq. (5)) and is simply Π δ>δc (a 1 , a) = 1 a)]. Therefore, fixing a 1 = a eq (since the first overdensities collapse at MRE given the remarkable coincindence in eq. (22)) are less and less produced -this is related to the suppression of P δ at large momenta as 1/k). 50 Interpolating between these two behaviours, there is a peak whose position and amplitude depend on the value of δ c . These features can be seen in Figure 6, where we show the estimate in eq. (67) for δ c = 0.22 and c M = 0.45 with σ(M ) calculated using the smooth window function W (k, k 1 , k 2 ) = W s (k 1 /k)W s (k/k 2 ), with W s (x) ≡ 1 2 (1 − tanh(2(x − 1))). In Figure 19 we show the prediction for the comoving number density of soltions produced per unit time, dΠ δ>δc (a 1 ,a) da ρ(a)/M (a), for the same W and δ c and we compare it with the simulation data. As we observed, our analytic argument reproduces the data in Figure 6 remarkably well. Moreover, as mentioned, it predicts that at late times the soliton production is suppressed, with a peak in df s / log M at M/M eq J 0.1 (not captured by the numerical simulation, because of the limited time range). In any case, given that the field is not Gaussian and the exponential sensitivity to δ c in eq. (67) (present in the Gaussian approximation), we refrain from a precise fit of the parameter δ c , or of the window function, and regard the agreement as mostly qualitative.
To check the reliability of our analytic method, we have carried out simulations with the same initial conditions as the power spectrum in eq. (9) but with k → x 0 k with x 0 = 1/2. As expected, this change to the initial power spectrum leads to a different (smaller) soliton mass function, plotted in Figure 20, and production rate, plotted in Figure 21 (left). It also leads to a slightly smaller value of c M 0.4, which determines the masses of solitons produced at a given time. We find that our analytic approach, with the same window function, also matches the data well in this case, albeit with a slightly different value of δ c 0.3. This gives us confidence that our analysis is capturing the main aspects of the underlying physics, although a further refinement would be useful in the future. Figure 20: The soliton mass function when starting from artificially rescaled initial conditions k → k /2, to test our analytic analysis. We find that our analytic approach works well also with these initial conditions (with a slightly different value of δ c ).

F Details of N-body Simulations
We carry out our N-body simulations using the code Gadget 4 [111], modified to include the effect of radiation on the Universe's expansion. For simplicity, we follow the standard convention of including only radiation, dark matter and vacuum energy in simulations (i.e. neglecting baryons). We fix the cosmological parameters Ω 0 = 0.311, Ω Λ = 0.689 and Ω rad = 0.0000924 (corresponding to z eq = 3370) and the present day Hubble parameter H = 67.7 kms −1 Mpc −1 . Note however that most of our results can be rephrased in terms of a/a eq in which case these numerical choices do not matter. The absence of baryons (which, as discussed, do not collapse into clumps around MRE but might be relevant later) will affect the growth of density perturbations slightly, but to the level of accuracy that we are aiming for this is not a major uncertainty. As discussed in Section 4, we start our simulations prior to MRE, at a/a eq = 0.01 (as with our Schrödinger-Poisson simulations the particlar numerical value has no effect provided a/a eq 1). The simulations are pure N-body with just gravitational interactions, so cannot capture the effects of quantum pressure. We set Gadget's time-stepping parameter to 0.02, which has been found to be sufficiently accurate in similar simulations [72]. Halos are identified using Gadget's friend-of-friend algorithm, with dimensionless link length 0.2, so particles are linked if their spacing is less than 0.2 of the mean particle spacing, and fix the minimum group length for a set of particles to be classed as a halo to 32.
We use a box size of 80λ , and obtain our results averaging over three simulation runs. The comoving gravitational softening length is chosen to be 0.003λ . This is smaller than λ J in the cores of the densest objects that form, at which scales N-body simulation results will not reproduce the true physics anyway, see Section 4.3. E.g. at z = 70, corresponding to the density profiles plotted in Figure 10, the gravitational softening length is about 2 · 10 4 km(m/ eV) 1/2 and λ J (5 eV 4 ) 3 · 10 4 km. Simulations are run from a/a eq = 0.01 until z = 24 (i.e. a/a eq 140). At this time the most IR modes have P δ (k) 0.05, Figure 21: Left: The soliton production rate with artificially rescaled initial conditions compared to the analytic prediction, analogous to Figure 19 for the physical initial conditions. Right: The average mass of solitons produced in different time intervals compared to the predicted parametric dependence M J (a i ) with the rescaled initial conditions. The results are similar to those obtained starting from the physical initial conditions, as plotted in Figure 5. and still closely match the values predicted by a linear analysis, indicating that finite volume systematics are not yet becoming important (we could run to slightly later times, but the statistical uncertainty is already starting to limit our results anyway). 51 In Figure 22 we plot the mass function of compact halos from N-body simulations comparing results from initial conditions with a UV cut and from realistic initial conditions. The initial conditions with a UV cut are a realisation of Gaussian density field with the power spectrum of Figure 2 (right) for k < 0.5k and 0 otherwise, converted to particles (for the m of interest only the isocurvature part of the initial spectrum is relevant given our box sizes). This choice is made to exclude the density fluctuations on scales that will collapse to solitons around MRE (and which will be affected by quantum pressure), see Figure 16 (right). The full initial conditions are generated from a realisation of the density field, generated from the power spectrum of A in Figure 2 (left). Of course in this case the dynamics of the small scale density perturbations are not correctly reproduced. As discussed in the main text, the agreement of the mass functions for M 50M eq J is an indication that the mass function in this interval is unaffected by the dynamics of modes at small scales.

G Analytic Analysis of the Compact Halos
In the Press-Schechter analysis of Section 4.2 we smooth the field over a distance R using a window function W (k, R) that vanishes at k 2π/R, i.e. δ s (t, x) = (2π) −3 d 3 kδ(t, k) exp(i k · x)W (k, R). We choose the top-hat window function W (k, R) = W t (kR) with W t (x) ≡ (3/x 3 )(sin x − x cos x). As mentioned in the main text, since δ is Gaussian at scales larger than λ , the probability that the field smoothed over R = (3M/4πρ) 1/3 is larger than δ c is determined only by the variance Π δ>δc = 1 dk/k P δ (t t eq , k)W 2 (k, R), with P δ at t t eq as in eq. (12). Figure 22: Left: Comparison between the halo mass function obtained in N-body simulations starting from a realisation of the full dark photon initial conditions (solid) and from initial conditions including only perturbations on scales k/k < 0.5, which are unaffected by quantum pressure (dashed). The difference between the results from these two sets of initial conditions gives an estimate of the importance of dynamics at small scales for the system's evolution, and so an estimate of the scales on which the results of N-body simulations without quantum pressure can be trusted. The number of objects with small mass M/M eq J 50 differ widely, but properties of heavier objects broadly agree. Right: The mass function of compact halos at different redshifts from simulations and from the Sheth-Tormen inspired fit to the data (ST), which reproduces the results well.
To obtain eq. (41) we estimate σ 2 s (M ) 3 3/2 D 2 [a]M /(16π 3 M ) using a crude approximation P δ (k) ( √ 3/π)(k/k ) 3 at k < k and otherwise zero (since modes with k k are affected by quantum pressure, and form solitons). 52 As mentioned in the main text, although the PS result captures the most important features of the dynamics, it is not precisely accurate. Instead, to reconstruct the halo mass function at z 15, we fit the data with an empirical Sheth-Tormen (ST) inspired form, which reproduces the data at early times more accurately than the PS prediction in eq. (41) (a similar fit was carried out in the context of axion miniclusters in [72]). The ST form for the halo mass function is where ν(M ) = δ c /σ s (M ) as before, q and p are parameters to be fit, and A −1 (p) = 1 + 1 √ π 2 −p Γ 1 2 − p is required by the normalisation condition dM df h /dM = 1. We obtain q = 0.696 and p = 0.305 (with δ c = 1.7 fixed), which fit the N-body results at early times well (see Figure 22 right). These parameter values are close to those that come out of the original analytic analysis [79].
Finally, we give complete formula for the NFW parameters of the compact halo profiles, described in of the host inside the clump's orbit, ρ host (r orbit ) = M (r orbit )/ 4 3 πr 3 orbit , and ρ(r) = M clump (r)/ 4 3 πr 3 is the mean density of the clump within the radius r.

Dynamical friction
If a clump is bound in a larger object, dynamical friction leads to the clump's orbit decaying and falling into the centre of the host object, where the clump will be destroyed by tidal forces. The orbit decays on a timescale [25,113] where t orbit is the time the clump takes to orbit the larger halo. Putting typical values of M clump into eq. (72), dynamical friction is irrelevant for any of the objects of interest inside the Milky Way. However, it could be important for e.g solitons inside compact halos, or compact halos inside the lowest mass halos from the adiabatic perturbations, for which the ratio M host /M clump is not too large.

Collisions with stars
If it is in a galaxy, a DM clump can be disrupted by the energy transferred from a passing star. Typically such an encounter happens fast compared to the dynamical timescale of the clump. In this case, known as the impulse approximation, the energy transferred can be straightforwardly calculated. More precisely, the impulse approximation is valid if r clump /b v rel /σ, where r clump is the radius of the object, b is the impact parameter of the collision, v rel is the relative velocity of the collision, and σ is the velocity dispersion of the clump (i.e. the variance of the magnitude of the DM velocity in the clump is σ 2 ) [114].
The fuzzy halos around solitons and the compact halos are extended objects with densities that vary by five orders of magnitude between their centres and their outer edges. It is therefore plausible that the outer edges might be disrupted by a collision with a star while the centre remains. To account for this possibility we make the approximation that energy transferred to the mass in a shell at some distance from the centre of the clump can lead to this shell being lost from the object, independently of what happens to the mass closer to the centre. Compared to the alternative estimate of requiring that the entire object is destroyed by an encounter, this means that the outer edges of objects are more easily destroyed. We discuss the disruption of solitons in Section H.2.
We consider only the possibility that a shell is removed by a single close encounter rather than a series of encounters. In the second case the energy transferred to a shell might be redistributed throughout the object. Ignoring this effect, and assuming the object and perturbing stars have a Maxwellian distribution, the timescale for destruction by sequential heating and by a single catastrophic encounter happen to coincide [69]. Therefore our neglect of destruction by sequential heating will not have a major effect on our results, although further analysis would be worthwhile (see e.g. [114] for related studies in the context of WIMP dark matter).
We define b crit to be the critical impact parameter, so that the object is destroyed if the impact parameter b < b crit and is otherwise unaffected. Given that the probability of an encounter with impact parameter b grows linearly with b the majority of destructive encounters have b not too much smaller than b crit (the average destructive encounter has b = 2b crit /3). For all objects of interest b crit r clump , and also the conditions for the impulse approximation to hold are satisfied. The energy transferred to a shell of dark matter, of mass dm, at radius r < r clump from the centre of the clump by an encounter with a star of mass M s is 53 We approximate that the shell will be removed from the object if the transferred energy is greater than the gravitational binding energy of the shell, GM (r)dm/r, where M (r) = 4π r 0 r 2 dr ρ(r ). Consequently, the critical impact parameter is given by where, as before, ρ(r) is the mean density of the clump within the radius r.
We are particularly interested in objects that pass the neighbourhood of Earth on their orbits. Most such objects are on trajectories that do not entirely reside inside the galactic disk. The cross section for a destructive encounter, πb 2 crit , is proportional to the mass of the disrupting star. Therefore, the probability that a shell at radius r is removed from the clump after the clump crosses the disk n times is where S is the mass density encountered per unit transverse area of the disk, provided p dest 1 (this analysis mirrors that applied to axion miniclusters in [87]). Consequently, p dest = 0.4 n 100 0.05 eV 4 ρ (r) with the reference value S = 140M pc −2 a factor of 4 larger than the column density transverse to the Milky Way disk to account for the fact that on average trajectories are not perpendicular to the disk [87]. n 100 corresponds to roughly the number of disk crossings since the formation of Milky Way for a clump on a typical orbit [87], however this can vary significantly for different forms of orbit.
Tidal shocking by the disk A dark matter clump in a disk galaxy, could also be destroyed by the overall gravitational field of the disk, rather than by encounters with individual stars. As a clump moves through the disk, the energy of DM particles that are away from the centre of the clump in the direction perpendicular to the disk increases compared to those at the centre. We consider the dark matter of total mass dm in a small region that is a distance ∆z away from the centre of the clump in the perpendicular direction. The relative energy increase of this DM is [116] (see also e.g. [117,118]) ∆E = 2(2πGσ s (r)) 2 (∆z) 2 v 2 z dm(1 + a 2 ) −3/2 , (77) 53 In the limit b R a similar expression can be obtained, and for b R an approximate expression that interpolates between the two limiting behaviours can be used, see e.g. [115]. Note that the expression for the energy change of a particular particle gets further contributions, but these average to zero and are approximately the same magnitude as the term that remains anyway.
where σ s (r) is the surface mass density of the disk, and v z is the velocity of the clump perpendicular to the disk. At 8kpc from the Milky Way centre, σ s 10 8 M /kpc 2 . The final factor in eq. (77) is a suppression due to adiabaticity of the process [119]. Intuitively, the destruction is less efficient if the tidal force changes slowly compared to the time a DM particle takes to orbit the clump. In particular, a = t crossing /t internal where t internal is the time the particle takes to orbit the clump, and t crossing ∼ 5 · 10 13 s is the time the clump takes to cross the disk [117]. In the Milky Way, a 6 ρ/ eV 4 1/2 , which is relevant for the solitons, but only marginally important for the halos surrounding them.
Similarly to before, we compare the energy increase in eq. (77) to the binding energy of the dark matter particles in the clump. A single crossing disrupts regions of the clump for which ρ(r) 10 −6 eV 4 σ s 10 8 M /kpc 2 Even if the energy transferred adds up every disk crossing, the disruption through this process is negligible compared to that caused by collisions with stars.

Tidal shocking
During hierarchical structure formation, as a DM clump becomes bound in a larger object it experiences tidal forces on timescales shorter than the dynamical time of the clump. These are known as tidal shocks [120,121], and are distinct from the steady state tidal forces analysed previously. During a typical tidal shock, a shell of dark matter of mass dm at radius r from the centre of a clump is expected to get a relative energy increase compared to the rest of the clump of ∆E 4π 3 γ 1 Gρ host r 2 dm , where γ 1 is expected to be of order 1 [122] and, as before, ρ host is the mean density of the host within the orbit of the clump. We make the crude approximation that the energy transferred to the DM shell is not redistributed to DM in other parts of the clump. Then the rate of energy increase of the shell isĖ = γ 2 ∆E/t dynam , where t dynam (Gρ host ) −1/2 is the dynamical timescale of the larger object (parameterically the same as the time the clump takes to orbit the host t orbit once it is bound, defined previously), and γ 2 parameterises the number of tidal shocks per dynamical time, which is expected to be roughly order 1. Hence the shell is likely to be disrupted on a timescale t shock ρ(r) √ Gγ 1 γ 2 ρ However, we stress that this estimate is very rough and a detailed analysis would be valuable in the future. We have little control of the parameters γ 1 γ 2 , and we do not know how long tidal shocking continues until the clump reaches a steady state orbit. Also it is not even certain if it is accurate to compare the energy increase per tidal shock, eq. (79), to the binding energy. For instance [123] finds that disruption is not inevitable even if the energy increase is much larger than the binding energy. Moreover, as discussed, in a realistic history of structure formation, a DM clump is likely to be bound in a series of objects of increasingly large mass. Consequently, the clump might not need to survive too long inside a slightly higher mass clump until that itself is quickly destroyed when it is bound inside a larger object.

H.2 Solitons
It is not clear if our approach of allowing the disruption of shells is an accurate approximation for solitons, because the de Broglie wavelength is comparable to the size of a soliton. Also, if the outer layer of a soliton were to be lost the profile of the remainder might change, affecting the subsequent destruction probability. In particular, the analysis of [110] suggests that the remainder of the soliton might spread out, making later destruction more likely. For definiteness we make the approximation that the soliton will be disrupted if sufficient energy is transferred to remove a shell at radius r 0.7λ J (ρ s ) (where the density profile switched from soliton like to NFW-like). If we instead required that the entirety of the soliton inside 0.7λ J (ρ s ) was destroyed, the destruction probability would decrease. We have seen above that the probability that a shell at distance r from the centre of a clump is destroyed depends on the mean DM density inside the shell, ρ(r), rather than the density at the distance r, which is ρ(r). At 2λ J (ρ s )/3, we have ρ 0.2ρ s . Solitons with ρ 0.05 eV 4 are likely to survive destruction by the Milky Way central potential, collisions with stars, and dynamical friction in the Milky Way. Therefore, we conclude that these processes are likely to leave the majority of solitons produced around MRE undisrupted, since they have central densities ρ s (0.1 ÷ 100) eV 4 (see Figure 6, corresponding to ρ (0.03 ÷ 30) eV 4 ).
The possible disruption of solitons during hierarchical structure formation is less certain. Most solitons become bound in compact halos before being subsumed in an adiabatic halo. The probability that the soliton inside a compact halo is destroyed depends on the mean density of the compact halo inside the soliton's orbit ρ host (r orbit ). The compact halos have mean densities in the range ρ(r) host ∼ (10 −6 ÷1) eV 4 , see Figure 23 (right) and masses between M halo (10 2 ÷ 10 4 )M J (a eq ).
Because that fraction of dark matter bound in compact halos rapidly increases at late times, the majority of the solitons will be bound in the most massive and least dense compact halos. Given the low density of the these compact halos, tidal stripping will typically not destroy any solitons they contain. Additionally, the timescale for dynamical friction in eq.
where M s is the mass of the soliton and M host is the mass of the compact halo that it is inside. Consequently dynamical friction will also not destroy a typical soliton in a fairly massive compact halo. Meanwhile, the timescale for destruction by tidal shocking, eq. (80) might be relevant for the lower mass solitons if they are contained in a relatively low mass (i.e. high density) compact halo. However, many of the compact halos will be destroyed as adiabatic structures form, in which case a soliton need only survive inside the compact halo until that happens. Similarly, a typical soliton is unlikely to be destroyed if it is bound into a relatively low mass adiabatic halo prior to ending up in the Milky Way. A detailed analysis will be needed to draw definite conclusions, and this would be worthwhile in the future.
So far we have neglected the wave-like nature of the dark photon inside a soliton. Although we do not attempt a full analysis of the effect of this on the solitons' survival, we note that this does allow the soliton to lose mass by tidal stripping even from the region inside r t defined in eq. (71). This has been analysed in [57], and considered in more detail in [110], which accounts for the fact that the soliton becomes less dense after mass is removed, accelerating the subsequent destruction. Using numerical simulations the latter reference finds that a soliton of mass M orbiting a distance r orbit will be unaffected by tidal stripping by a central potential from a halo with mean density ρ host (r) if the central density of the soliton ρ s satisfies ρ s /ρ host (r orbit ) 100.
For a soliton in the Milky Way ρ(r orbit ) 5 · 10 −5 eV 4 , so given that the typical soliton central densities are 0.1 ÷ 100 eV they are unaffected by this effect. Likewise this will not affect a typical soliton that is bound in a compact halo.
Overall, although there are many uncertainties in our analysis, we conclude that it is plausible that the majority of the solitons survive to the present day.

H.3 Fuzzy halos
Next we consider the fuzzy halos around the solitons. If such objects lie on trajectories that cross the neighbourhood of Earth, we can see from the analysis in Section H.1 that destruction by stars is the most important of the effects of the Milky Way. From eq. (76), we expect that a halo is likely to survive out to a radius r such that ρ(r) 0.05 eV 4 . In Figure 23 (left) we plot ρ(r) as a function of the local density ρ(r) for the solitons and their surrounding fuzzy halos (both normalised relative to the soliton's central density ρ s , so that the plot applies all solitons and fuzzy halos, regardless of their mass). Values of ρ(r) 0.05 eV 4 in the range of interest corresponds to local densities of ρ(r) 10 −2 eV 4 . Outside these r the fuzzy halos are likely to be destroyed.
The halos could also be disrupted if it becomes bound in a compact halo prior to adiabatic structure formation (or in a low mass adiabatic halo prior to being bound in the Milky Way). Considering only the region with ρ(r) 0.05 eV 4 , which has a chance of surviving interactions with stars, this part of the fuzzy halo is unlikely to be outside the tidal radius given by eq. (71) when it is bound in a compact halo.
Meanwhile, the typical mass contained in a fuzzy halo inside ρ(r) 0.05 eV 4 is roughly a factor of 10 ÷ 100 greater than the soliton mass, so typically a factor of 10 4 ÷ 10 3 smaller than the relatively high mass compact halos. Therefore, from eq. (81) dynamical friction might potentially be relevant. Likewise, tidal shocking, see eq. (80), could be important. However, as previously discussed the compact halos are themselves likely to be quickly destroyed. Finally we note that, as with the solitons, wave-like effects could be relevant for the fuzzy halos, and a detailed analysis of such effects would be worthwhile in the future.

H.4 Compact halos
The values of ρ(r) as a function of ρ(r) for typical compact halos are plotted in Figure 23 (right). Like the halos around solitons, for a compact halo in the Milky Way the regions outside the radius r such that ρ(r) = 0.05 eV 4 are likely to be destroyed by collisions with stars. As a result, the majority of a compact halo with mass 10 3 M J (a eq ) is likely to be disrupted, possibly leaving a dense core behind. This means that the compact halos that contain the majority of the dark matter mass when adiabatic structure formation starts are likely to be mostly destroyed. Compact halos with masses ∼ 10 2 M J (a eq ) could survive to the present day, and these contain ∼ 10% of the dark matter when adiabatic structure formation starts (see Figure 8). A small mass compact halo might also survive as substructure once bound inside a larger mass compact halo and survive to the present day.
�� � ρ(�)/ρ ����� Figure 23: Left: The mean density inside a radius r, ρ(r), for a soliton and its surrounding fuzzy halo as a function of the local density at r, ρ(r) (normalised to the density at the soliton core ρ s ). These results are obtained by matching the soliton profile and NFW fit shown in Figure 7. On the upper axis we indicate the corresponding enhancement over the local dark matter density, fixed to 0.5 GeV/ cm 3 , which depends on ρ s . Regions of the fuzzy halo with ρ(r) 0.05 eV 4 are likely to survive undisrupted to the present day in the Milky Way. Right: The mean density inside a radius r, ρ(r), as a function of the local density ρ(r) for compact halos of different masses. These results are obtained from the NFW fit described in Section 3.7. The solid lines extend to densities where quantum pressure becomes relevant. The upper axis shows the corresponding enhancement over the local dark matter density. As with the fuzzy halos, regions with ρ(r) 0.05 eV 4 are likely to survive to the present day.