Analysis of Bose-Einstein condensation times for self-interacting scalar dark matter

We investigate the condensation time of self-interacting axion-like particles in a gravitational well, extending the prior work [arXiv:2007.07438] which showed that the Wigner formalism is a good analytic approach to describe a condensing scalar field. In the present work, we use this formalism to affirm that $\phi^4$ self-interactions will take longer than necessary to support the time scales associated with structure formation, making gravity a necessary part of the process to bring axion dark matter into a solitonic form. Here we show that when the axions' virial velocity is taken into account, the time scale associated with self-interactions will scale as $\lambda^2$. This is consistent with recent numerical estimates, and it confirms that the Wigner formalism described in prior work~\cite{Relax} is a helpful analytic framework to check computational work for potential numerical artifacts.


I. INTRODUCTION
Axions and axion-like particles have become some of the most well-motivated candidates for the dark matter. Axions across a wide mass range could have the correct abundance to compose the dark matter [2][3][4][5]. Additionally axions can address a number of longstanding problems such as the strong CP problem in QCD, and they can play a role in inflation and as dark energy [6][7][8][9][10][11][12].
As a dark matter candidate, axions can solve problems with small-scale structure formation which are difficult to explain under the dominant weakly interacting massive particle (WIMP) model [13][14][15]. While on cosmological distance scales the predictions of axion models agree with the predictions of WIMPs, axions behave differently on the length scale of their de Broglie wavelength. This not only addresses the small-scale structure problems of ΛCDM, but it could also lead to novel astrophysical structures which could one day be observed, for example through gravitational wave signatures from collisions with other objects [16][17][18].
One such structure is gravitationally bound solitons which are believed to form in regions of axion overdensity [19][20][21][22][23][24][25][26]. These solitons are often called Bose stars since the axion field condenses into a Bose Einstein condensate (BEC). While there have been many studies of these objects, especially their equilibrium properties, there are still questions about their formation process. Since relaxation into the BEC phase is driven by the interactions between particles, and the two interactions (gravitational and self-interactions) at play in the axion field are extremely weak, this process can take an exceedingly long time. Recently, it has become clear that gravity alone is sufficient to drive the condensation process [23][24][25][27][28][29][30], but there are still questions as to what role, if any, the axion self-interaction plays in this process, partic- * kkirkpat@illinois.edu † aem8@illinois.edu ‡ Chanda.Prescod-Weinstein@unh.edu ularly in axion species where these self-interactions are stronger than in QCD axions.
In this work we evaluate the effect of the initial conditions of the axions in virialized minclusters on the nonequilibrium condensation process. In Ref. [1], we previously showed that the evolution equation of the Wigner function describes the relaxation process, and found that gravitational scattering alone was sufficient to cause condensation. In this paper, in the opposite regime where self-interactions drive condensation, we analytically derive results that confirm previous estimates [25]. Specifically, using the analytic framework described in [1] we show here that when the limit on the axions' energy set by the virial velocity is taken into account, the relaxation rate scales with λ 2 . These new results update our prior ones.
This paper is organized as follows. In Sec. II, we evaluate the role of self interactions in the relaxation by expanding the equation of motion for the axion field's Wigner function. In Sec. III, we conclude and compare the impact of self-interactions and gravity on the condensation process.

A. Initial conditions in miniclusters
The QCD axion undergoes a phase transition when the universe cools below the Peccei-Quinn symmetry breaking scale, f a ∼ 10 12 GeV. Depending on the energy scale of inflation, this can either occur during inflation or after it has ended. In the post-inflationary scenario, which we consider here, f a is lower than the energy scale of inflation, and the result is the formation of gravitationally bound structures through the Kibble mechanism [12]. Symmetry breaking causes the axion field to take nonzero values. However, these values can differ in causally disconnected regions of the universe, resulting in a field with random, O(1) fluctuations on a distance scale set by the size of the Hubble volume during symmetry breaking.
The size of the miniclusters depends on the symmetry breaking scale f a and hence implicitly on the mass of the axion. For the QCD axion, the typical minicluster size is R ∼ 10 5 km [34]. Ultralight axions and other scalar dark matter can form dark matter clusters and halos through the same mechanism, though the size of these clusters depends on the model in question [12].
The typical velocity of the field is also determined by these initial conditions. Since the miniclusters are gravitationally bound objects, the velocity is determined by the virial theorem, which implies For the QCD axion, the virial velocity within the miniclusters is v ∼ 10 −10 c. The result of this process is that the axion field ψ(x, t) is only supported within a region of radius R. Similarly, the wavenumbers which are occupied by the axion field are limited. The de Broglie wavelength corresponding to the virial velocity provides an upper bound, and the result is that the field is occupied only for wavenumbers in the range

B. Evolution due to self-interactions
The virial velocity v is much smaller than c, so the field is non-relativistic. In the non-relativistic regime, the real scalar field φ that describes the axion can be replaced with a complex scalar field ψ through (Throughout the paper, we suppress the timedependence of all fields and write, e.g., ψ(t) = ψ.) Because it is non-relativistic the field ψ evolves under Newtonian potential U due to the field's self-gravity. It also experiences self-interactions arising from the QCD axion potential, In the non-relativistic limit the field is much lower than the symmetry breaking scale, and the cosine potential can be expanded as with Ultralight axions and other scalar field dark matter also experience quartic self-interactions, though the relation between the coupling constant λ and the other parameters in the model may differ in those cases. The evolution of the non-relativistic field ψ under both these interactions is determined by the Gross-Pitaevskii-Poisson (GPP) equations, where the background density n 0 is removed as a source in the Poisson equation [22,27]. Here we have introduced a dimensionful coupling constant It is useful to consider the evolution of the Fouriertransformed field in momentum space, which we will write as ψ p . To obtain the momentum-space equation of motion for ψ p , one differentiates the Hamiltonian with respect to ψ * p . The Hamiltonian is where the kinetic energy is and the interaction Hamiltonian is, In general, there is an additional interaction term for the gravitational interactions between particles. In this paper we study the effect of the self-interactions, so we turn gravity off and compare condensation due to gravity with condensation due to self-interactions in the discussion. The momentum space equation of motion is obtained by differentiation where the self-interaction term is The Wigner function (sometimes called the Wigner quasiprobability distribution), which characterizes a statistical ensemble of fields, is f (x, p) = dy e −ip·y ψ * (x + y/2)ψ(x − y/2) The angle brackets · are an ensemble average over the distribution, which we take to be a distribution of Gaussian fields with random phases. The Wigner function describes the density of the axions in phase space. Unlike a distribution of classical particles, it can take on negative values in phase space regions with area less than . The phenomenon is the result of interference between waves. The equation of motion for the Wigner function can be obtained by applying Eq. (13) to each field in the integrand of Eq. (15). We obtain, So we have Assuming the fields are spatially homogeneous, we can take x = 0, and the first term vanishes. Introducing notation for the four-point correlator, we can write the equation of motion as Now under the assumption of an initial Gaussian distribution, the four-point function factors into a sum of two-point functions (the Wick contractions) at t = 0. If we assume this holds at all times, rather than only at t = 0, we obtain df dt (x, p) = n p Im dq gn q , where is the momentum space density. which vanishes since the interaction g is real. We see that under this assumption about the four-point function J, the derivative of the Wigner function vanishes, that is, the distribution appears constant.
In reality, the four-point correlator J changes with time, and since the distribution does not remain Gaussian, the connected correlations become relevant [35]. By applying the Gross-Pitaevsky equation (Eq. (13)) to each of the four fields in in the four point correlator J, we obtain an evolution equation for J, in the same way that we obtained the evolution equation for the Wigner function in Eq. (5). We get where A p1p2p3p4 is a six-point correlation function. This arises in the evolution equation of J due to the nonlinear term in the Gross-Pitaevskii equation. Due to our initial conditions of a Gaussian distribution with random phases, this six-point function can be factored into its Wick contractions, with negligible connected component, A p1p2p3p4 = n p3 n p4 (n p2 + n p1 ) − n p1 n p2 (n p3 + n p4 ).
(24) As with the four-point function J, the connected component of this six-point function becomes non-negligible as the system evolves. However, this process is suppressed by another factor of the interaction strength g so we will neglect these corrections to lowest order in g, and simply assume that the six-point correlation A is constant.
We can solve the evolution equation Eq. (22), obtaining We obtained this expression in order to estimate the rate of change of the Wigner function, and so we need to evaluate the integral in Eq. (19). The first term in J is an oscillating function. When t is large, its contribution to the integral over momenta becomes negligible due to the cancellation of phases. Specifically, we need t to be larger than the typical value of (∆ω) −1 , and therefore larger than the typical value of ω −1 p . Since the axions exist in virialized miniclusters, the typical frequency of the field is defined by the virial velocity v of the minicluster, namely So provided the time is much greater than the typical timescale for the variation of the occupation numbers, then the contribution from the oscillating term is negligible.
Physically, the fourth-order correlations J begin at zero since the initial distribution of axions in the minicluster is a distribution of Gaussian waves. Due to the nonlinear self-interactions, they cannot remain at zero but must grow as the self-interactions drive the system away from its initial Gaussian distribution. As we see in Eq. (25), they grow and oscillate around their steady state value gA/∆ω, on a timescale set by the typical energy of virialized axions in the minicluster, mv 2 , where v is the virial velocity. If we look at this process on a timescale in the regime of Eq. (27), then the dynamics of the fourth-order correlation function can be neglected, and we can assume that J has a constant, nonzero value from the beginning of the condensation process. For a realistic minicluster, this period is much shorter than the relaxation timescale that follows from this assumption, so the approximation is consistent.
Substituting our solution for J into Eq. (19) and neglecting the contribution from the oscillating term, we obtain From this kinetic equation, we can obtain a rate for the evolution of the Wigner function towards its steadystate.
In particular, the relaxation rate scales with the square of the self-interaction strength.

III. DISCUSSION
In Ref. [1], we introduced a formalism for analyzing the relaxation of the scalar field dark matter into its condensed state when both self-interactions and gravity are present, by analysing the evolution of the Wigner function. We previously estimated the relaxation time due to self-interactions without accounting for the upper limit on the momentum occupation in the axion field imposed by the virial velocity. When this natural cut-off is taken into account, we see that the growth of the four-point correlations is suppressed, never growing past a steady state value. The resulting relaxation rate is proportional to g 2 . Since the self-interactions are small, this is gives a much longer relaxation time than calculated in Ref. [1]. The timescale estimated here agrees with the timescale for condensation due to self interaction found in the recent simulations in Ref. [36].
When both gravitational and self-interactions are present the equation of motion for the Wigner function has the form where τ λ is the relaxation timescale due to selfinteractions alone, τ G is the relaxation timescale due to gravity alone (calculated in the section above), and the combined timescale is This simple form is seen in the Wigner function formalism. In particular, there is no term that drives the Wigner function of the form ∂ t f ∼ gGf . Since the gravitational and self-interactions operate on different distance scales, such cross terms do not arise at lowest order in the interaction strengths g and G, and the two interactions can be analyzed separately. In physical scenarios relevant for QCD or ultralight axions, τ G τ λ , so the self-interactions do not significantly speed up the relaxation process.
In Ref. [1], we previously predicted This equation occurred because the timescale for relaxation due to self-interaction arose out of a second-order equation for the Wigner function when we ignored the effect of the virial velocity in limiting the occupied momentum states in the axion field. When this effect is included, both self-interactions and gravity affect the evolution of the Wigner function at the same order, and the correct combined timescale is Eq. (32). As shown in Fig.  1, the two formulas approximate each other when either the relaxation time from self-interactions or from gravity greatly exceeds the other. While both expressions qualitatively fit the recent data obtained through numerical simulations in Ref. [36], we believe Eq. (32) fits the data better than Eq. (33), since it is more appropriate to this physical situation and since it reduces the apparent systematic difference in the fit to the data in Ref. [36].
Ref. [36] presented numerical simulations arguing that the correct time constant for condensation due to selfinteractions scaled with g −2 . This work reinforces that result, and is consistent (up to an O(1) constant) with relaxation time caused by a self-interaction cross-section, that was shown in the simulations and predicted in Ref. [28] within kinetic theory. Since these timescales for condensation into the Bose star are the scales associated with exponential growth of the correlations, they are only defined up to O(1). Finally, we note that this analytic work confirms that for QCD axions and other scalar dark matter with weak self-interaction coupling strength, the relaxation timescale is set by gravity, since the corresponding grav-itational timescale found in Ref. [28] and confirmed by simulations [30] is: For QCD axions in miniclusters, this gives a ratio τ λ /τ G > ∼ 10 12 .