Relaxation times for Bose-Einstein condensation in axion miniclusters

We study the Bose condensation of scalar dark matter in the presence of both gravitational and self-interactions. Axions and other scalar dark matter in gravitationally bound miniclusters or dark matter halos are expected to condense into Bose-Einstein condensates called Bose stars. This process has been shown to occur through attractive self-interactions of the axion-like particles or through the field's self gravitation. We show that in the high-occupancy regime of scalar dark matter, the Boltzmann collision integral does not describe either gravitaitonal or self-interactions, and derive kinetic equations valid for these interactions. We use this formalism to compute relaxation times for the Bose-Einstein condensation, and find that condensation into Bose stars could occur within the lifetime of the universe. The self-interactions reduce the condensation time only when they are very strong.


I. INTRODUCTION
The composition of dark matter is one of the most longstanding problems in cosmology. The dominant model, known as Lambda Cold Dark Matter (ΛCDM), proposes that the dark matter is cold and has a low velocity dispersion. It has been successful at cosmological distance scales [1]. However, at galactic distance scales and smaller ( < ∼ 10 kpc) it has a number of problems. At these scales, the predicted density profiles disagree with observations and a higher abundance of dwarf galaxies is predicted than is observed [2][3][4]. While there are several proposed solutions to these problems [5][6][7][8][9], an attractive proposal considers the quantum properties of the dark matter particles. In this case, the large-scale predictions remain the same as in ΛCDM, but on scales less than the de Broglie wavelength the predictions change.
Among the proposed candidates for the dark matter are light bosons, such as the QCD axion [10][11][12][13][14][15] or fuzzy dark matter composed of ultra-light axions or other scalar fields [16][17][18]. The QCD axion is especially wellmotivated since it is hypothesized as a solution to the strong CP problem in QCD, and it has been shown that it could account for the correct dark matter abundance [19,20]. Both the QCD axion and ultralight scalars are the subjet of several ongoing experimental searches, most notably ADMX [21] and ABRACADABRA [22]. Experiments into the neutron electron dipole moment can also constrain models of axion physics [23].
These proposed candidates have in common that they thermalize to form compact, gravitationally bound solitons that can be described as Bose-Einstein condensation [24][25][26][27][28][29][30][31]. This phenomenon is proposed to occur over a wide region of the parameter space of masses and selfcouplings, with masses ranging from m ∼ 10 −22 eV in There are two main scenarios that can generate structure formation of axion miniclusters or halos. For the QCD axion, the Peccei-Quinn symmetry breaking could occur after inflation. This happens when the symmetry breaking scale f a is less than the energy scale of inflation. For QCD axions this is possible since its symmetry breaking scale is lower, f a ∼ 10 12 GeV, but for ultralight axion-like particles, the symmetry breaking scale is too high, f a ∼ 10 16 GeV [44].
In this scenario, the miniclusters are formed through the Kibble mechanism. Symmetry breaking causes the axion field to take random, uncorrelated values in different Hubble patches, resulting in density fluctuations of O(1) which decouple from the background Hubble expansion to form miniclusters [26,[45][46][47]. Compared to the dark matter halos of ΛCDM, these are quite small with masses on the order of 10 −13 M ⊙ , and radii on the order of 10 5 km, determined by the mass and size of the horizon at the QCD phase transition [48]. This distance scale also determines the characteristic wavenumber of the axions, since initially the de Broglie wavelength of the axions is the size of the Hubble patch. However, the wavenumber is redshifted after the QCD phase transition.
In the second scenario, CP symmetry can be broken during inflation. This scenario occurs when the symmetry breaking sale f a is higher than the energy scale inflation, which could be the case for QCD axions or ultralight scalars. In this scenario the axion field in our current universe originates from a single Hubble patch at the time of the symmetry breaking, and so does not exhibit O(1) fluctuations since these are inflated away. Fluctuations in the axion field that generate structure can still arise in a number of ways: by gravitational collapse due to the Jeans instability [44], by the growth of quantum fluctuations in the axion field enhanced by the axion's self-interactions [49], or by a phase transition in the sector determining the axion's mass [50]. Recently, Ref. [51] found that this scenario does not result in minicluster formation for the QCD axion, but can lead to minicluster formation for other axion like particles.
In this paper, we refer to any gravitationally bound structure of axion or ultralight scalar dark matter, formed pre-or post-inflation, as an axion minicluster. While these structures have different masses, sizes, and observational signatures, they all consist of scalar dark matter with the potential to form solitonic cores such as Bose stars, where the scalar field is in its ground state.

B. Gross-Pitaevskii-Poisson Equations
In this section, we review how the axions or scalar dark matter can be described by a classical complex field evolving under a system of equations known as the Gross-Pitaevskii-Poisson equations. Axion or scalar dark matter is a scalar field φ. In the case of the QCD axion, it arises as the Goldstone boson for a spontaneously broken symmetry, which in the instanton approximation results in the following potential, Here m is the mass of the axions and f a is the symmetry breaking scale of the axion. When φ ≪ f a , as is the case for cosmological axions, we can expand the potential to fourth order to obtain where is the attractive quartic self-interaction. For QCD axions, the potential is determined entirely by the symmetry breaking scale f a since the mass and symmetry breaking scale are related by m ∼ 6 × 10 −10 eV Thus the self-interaction λ and the mass are not independent parameters. For generic string theory axions, the potential is also a periodic function of the field with period 2πf a , so it has the same expansion to fourth order as in Eq. (2). However, the mass and axion decay constant are not related as in Eq. (4), so the potential has two independent parameters in this case. The density of axions is extremely high, compared to the characteristic particle volume set by the de Broglie wavelength λ dBr . For example in our galaxy [27] the density of QCD axions is estimated to be The de Broglie wavelength of virialized particles in an axion minicluster depends on the size of the minicluster, but for axions in our galaxy we have Thus the occupancy number is N ∼ n gal λ 3 dBr ∼ 10 26 .
Under these high occupancy conditions, the coherent state axion dynamics can be approximated by a classical non-relativistic field. The approximation as a classical field is valid because the quantum fluctuations in a coherent state depend inversely on the occupancy number. If we expand where δφ is the quantum fluctuations of φ about the mean field, then Equivalently, the timescale on which the quantum evolution of the field differs from the classical mean-field evolution is extremely long, orders of magnitude greater than the lifetime of the universe [52].
Since the field is well approximated by a classical field, we can write down a classical action that couples the field to gravity. We consider only minimal coupling, where g µν is the metric tensor. The classical equations of motion for this action are the Euler-Lagrange equations, In the Newtonian limit, the metric takes the form where U is the Newtonian gravitational potential (to linear order) and a is the Hubble scale factor. This results in an equation of motion for φ whereφ is ∂φ/∂t and the Hubble parameter is H =ȧ/a. Finally, we take the non-relativistic limit by writing the real scalar field φ in terms of a complex scalar field ψ as In the non-relativistic limit, the phase factors e ±imt oscillate rapidly and when we substitute Eq. (14) into Eq. (13), we can drop all terms with such a phase factor since they lead to only subdominant correction [53,54]. The result is the equations of motion for the complex field ψ, known as the Gross-Pitaevskii-Poisson equations (GPP) or the Nonlinear Schrödinger-Poisson equations In the above Poisson equation, we subtract the mean density n for consistency [55,56]. Here, U is the Newtonian gravitational potential, λ is the self-coupling, and n is the average density of axions in an axion minicluster.

C. Wigner distribution
Recently, Levkov, Panin and Tkachev [43] gave an argument a statistical ensemble of axions evolving under their self-gravity (without self-interactions) could not be treated as a standard Boltzmann collision process. They showed that since gravitational interactions are long range and interactions between distant axions are significant, the mean free path of the axions a gr is very small with respect to n 1/3 since long-range fluctuations are important. I.e., we have a dimensionless ratio a gr n 1/3 ≫ 1.
This implies that a Boltzmann collision process is inappropriate for modeling the gravitational interactions: the particles are too dense to treat collisions as a process involving only two particles.
We provide an additional argument that the Boltzmann collision integral is not valid, even for the short range self-interactions. Even though the mean free path of self-interactions is small with respect to n 1/3 , the de Broglie wavelength is not, and we have another dimensionless ratio of length scales, This is a restatement of Eq. (7) in terms of a dimensionless ratio of length scales. It tells us the quantum occupancy number N is high. This high occupancy number implies that the evolution of the statistical ensemble cannot be described by a standard Boltzmann equation because the axions cannot be localized to a definite position and momentum in phase space. Rather than describing an ensemble of particles by a phase space density, we can describe the ensemble by the Wigner function, The Wigner function is the closest mathematical object we have to a phase space description for an ensemble of waves. In the appropriate limit, when the occupancy number becomes low, it recovers the properties of a positive-valued probability density function of particles. In our case where the occupancy number is high, the Wigner function reflects the inability to localize particles by taking negative values on regions of phase space whose size is on the order of (i.e., on length scales set by the de Broglie wavelength and momentum scales set by the characteristic momentum). The negative values obtained by the Wigner function are the result of interference of the waves, a phenomenon that is neglected in a classical particle description. As a result, the Wigner function has been used to study quantum properties scalar fields during inflation [57,58]. The standard Boltzmann collision integral, which is developed for a localized collision of two or more particles, is not suited to describe the evolution of the Wigner function for the reasons stated above. Instead, we can systematically develop a kinetic equation by evolving the Wigner function by the GPP equations in Eq. (15), which relates the evolution of the Wigner function to a four-point correlation function. (U tot is defined in Eq. 23. For a derivation of this equation of motion see Appendix A.) Here, we have simplified notation by denoting The total derivative is By estimating the size of the four-point function, we obtain physical quantities such as the timescale for relaxation into the BEC state. Many numerical methods have been developed to study the evolution of the Wigner distribution in phase space [59].

A. Equation of motion
The gravitational potential U (x) is a functional of second order in ψ, the same order as the potential for nonlinear interactions, λ|ψ| 2 . These combine to form a single potential U tot and we rewrite the Gross-Pitaevskii-Poisson system as where and ∆ −1 x−x ′ is the Green's function for the Poisson equation, This can be expanded through Wick's theorem, Here the subscripts refer to the spatial arguments of the fields in the four-point function, the first two terms are Wick contractions, and the last term is the connected correlation function, which is nonzero whenever the distribution is not Gaussian. We assume that the initial distribution of the field in an axion minicluster is Gaussian, with randomly distributed phases. This is appropriate for an uncorrelated, but gravitationally bound system like axion miniclusters immediately after their formation through the Kibble mechanism [60]. When the initial state is a Gaussian distribution, the connected correlation function in Eq. (25) vanishes. As the ensemble evolves, the interactions cause non-Gaussianities to develop and the connected correlations grow at rates set by λ and G, as these are the coefficients of the nonlinear terms in the GPP equations.
Expanding the four-point function in Eq. (19) we obtain three factors from the Wick expansion in Eq. (25), where F 1 , F 2 are the contributions from the Wick contractions, To save space we have indicated the variables of integration as subscripts, so e.g., y = dy. The scattering integral I(f ) depends on the connected correlator and will be addressed in Sec. IIIB.
In the case of a Gaussian ensemble of random waves, all three of these factors vanish. The scattering integral I(f ) is proportional to the connected correlations, so it vanishes because the initial distribution is Gaussian. In the terms F 1 and F 2 , the Wigner functions, the Poisson Green's function, and the delta function are even in y, y ′ , but the imaginary part of the integral selects the odd component of the potential. Thus the integration over the spatial coordinates y, y ′ causes these terms to vanish as well. Under these initial conditions, the distribution of axions is initially static.
In order to obtain any timescale for the evolution of the ensemble of waves, we need to look at the second derivative of the Wigner function. We will differentiate the four-point function as before using Eq. (22), and expand the resulting six-point function through Wick's theorem. However, because the gravitational and self interactions operate on different distance scales, we must evaluate this second derivative for the two interaction terms in a different way. We can do this because the correlation function is linear, and the terms proportional to λ and G can be separated.

B. Landau scattering integral
The connected correlation function in the scattering integral is linear, so the terms governed by λ and G can be separated.
The term dependent on the gravitational potential U contains long range interactions. In the context of plasmas with Coulomb interactions, Landau first noted that it is dominated by fluctuations at long distances (compared to the de Broglie wavelength in this case) [61]. More recently in Ref. [43], Landau's analysis was adapted to Newtonian gravity. Near x the potential has a multipole expansion where y = |y|. In an axion minicluster of radius R, the the field and the potential are nearly homogeneous on scales much shorter than R. So when y ≪ R we can truncate the multipole expansion at first order. Now in the integrals that follow, U appears next to correlation functions of the form ψ(x + y/2)ψ * (x − y/2) , so the integrand is largest when the distance between the fields, y, is not much bigger than the correlation length 1/mv. But in the kinetic regime so the multipole expansion of U is valid to lowest order. Also, the equation of motion depends only on the part of U that is odd in y, so we get Following Landau, we write the gravitational scattering integral as a diffusion process in phase space, in terms of a Landau Flux, where F is determined by the four-point connected correlator In this equation, ψ ± are defined by Eq. (20), while we also introduce the shorthand, (throughout this paper, primes are shorthand for the arguments of functions, not derivatives). With these definitions, we see An evolution equation for the four-point function F can be obtained from the equations of motion Eq. (22) as before, This is a six-point function which we can expand into Wick contractions. It also has a connected component which we will neglect since it introduces additional factors of λ and G. We again use the multipole expansion to lowest order since the integral is dominated by short separations y, y ′ ≪ R. Finally, we solve the ODE to obtain an expression for the Landau flux valid at time t, In these equations, , p), and we use Einstein summation over repeated indices. Due to the logarithmic divergence from the Poisson Green's function, the integral over t ′ must be regulated at both long and short time scales. The axion minicluster has a radius R beyond which the axion field vanishes. Thus the long-time cutoff is R/v. Because the diffusion process is sensitive only to fluctuations at long distance, there is also a short-time cutoff 1/(mv 2 ), since the axion field lacks fluctuations at scales less than the de Broglie wavelength. This completely suppresses the contribution from the self-interactions, which operate only at short distance scales.
Physically, the relaxation rate due to gravity grows with time as fluctuations at further distances begin to interact. However, it does not grow without bound, since the minicluster has a finite size R. For times t > ∼ R/v we can apply these short-and long-time cutoffs to the integral in Eq. (41) to obtain the relaxation rate in this regime. As a result the integral depends on the Coulomb logarithm Λ = log(mvR).
The condensation timescale due to gravity is the inverse of this rate. From this expression for the Landau flux, we can estimate the rate of change for the scattering integral due to gravitational interactions. I.e. we have where C. Relaxation from self-coupling The self-interactions are not long range, so the scattering cannot be treated as a diffusion process in phase space like the scattering caused by gravitational interactions. We differentiate the equation of motion for the Wigner function in Eq. (19), replacing U tot with λ|ψ| 2 in that equation, since we have already treated the gravitational interactions in the previous section. This results in a six-point function which we consider to leading order in λ and G, (therefore neglecting its connected correlations). The Wick contractions yield six contributions. In the case of a homogeneous ensemble of random waves, most of these factors vanish. But we are left with one non-vanishing contribution, We recover a relaxation rate associated with the self interactions, Unlike the gravitational relaxation rate, the rate due to self-interactions does not depend on the distance between fluctuations, since the self-interactions are local. The relaxation rate does not grow with time until the connected correlations become significant. Instead we have a relation where d 2 f /dt 2 is directly proportional f , with proportionality constant γ. The constant γ has units of s −2 , so to obtain the relaxation timescale from self-interactions, we take a square-root,

IV. DISCUSSION
There are two qualitative differences between the two timescales for Bose condensation that arise from the non-local nature of the gravitational interactions. First, the gravitational relaxation time depends explicitly on the size R of the axion minicluster and on the axion's de Broglie wavelength 1/(mv), through the Coulomb logarithm Λ = log(mvR), while the self-interaction relaxation time does not. This is a straightforward consequence of non-locality. The coupling "constant" of gravitational interactions depends on the momentum k, so these two natural distance scales appear as cutoffs in the logarithmic integral in Eq. (41). Second, the gravitational relaxation rate is proportional to G 2 , while the self-interaction relaxation rate is only proportional to λ. This is also a consequence of non-locality. For the local self-interactions, we obtained a relation between the Wigner function and its second derivative in Eq. (46). The rate γ in that equation is proportional to λ 2 , but because this is an expression for the second derivative, we must take a square root to obtain the characteristic timescale associated with this process, resulting in a timescale that is inversely proportional to λ.
For the non-local gravitational interactions, it works out differently. Because the integral over t ′ in Eq. (41) does not go all the way to t but is regulated by the shortand long-time cutoffs, we do not get a relation between the Wigner function and its second derivative. Instead we have a relation between the Wigner function and its first derivative, and we do not need to take the square root of the rate in Eq. (43), so the timescale associated with this process is inversely proportional to G 2 .
These qualitative differences between the two timescales, as well as the difference in the strength of the coupling constants, leads to significantly longer relaxation times due to self-interactions than due to gravity. We find that for QCD axions, the relaxation timescale for gravity is substantially shorter than the timescale for self-coupling, but not to such an extreme degree as has been previously reported. Since λ only appears to first order in the relaxation rate while G appears to second order, the small self-coupling strength does not increase the relaxation time as strongly. We find In this case, Bose stars can form just within the lifetime of the universe due to their self-gravitation, while the self-interactions are too weak to have any effect during the formation process. For ultralight scalar dark matter, the de Broglie wavelength can be comparable to the size of the minicluster, mvR ∼ 1 which strongly affects the gravitational relaxation time sensitive to the Coulomb logarithm Λ. Moreover the mass and self-interaction are not determined by a simple relation like Eq. (4) for QCD axions, so there are more parameters which can vary. Assuming that the Coulomb logarithm is O(1), and taking the masses and interactions suggested by cosmological constraints in Ref. [62], we find that the condensation can occur much faster, though self-gravitation still dominates, τ λ ∼ 10 17 s.
When both gravity and self-interactions are present, the relaxation rate is simply the sum of the two rates, since at lowest order in λ and G there are no cross terms. Thus the total relaxation time is When either timescale vastly exceeds the other, this reduces to the more familiar form, As we have seen, the gravitational relaxation typically occurs much faster, so this expression reduces further to This proves that the formation process of Bose stars is dominated by gravitational interactions. By the time self-interaction have an effect on the fields evolution, gravity has already caused the field to condense. Finally, we note that while these calculations show that self-interactions of strength predicted for the QCD axion or most other scalar dark matter play a negligible role during the formation of the Bose-Einstein condensate, they can still play an important role in the phenomenology of the Bose stars. For example, Ref. [27] showed the sign of the self-interactions can determine whether longrange correlations are possible, with such correlations impossible under attractive self-interactions. Ref. [33] showed that the scattering length of self-coupling determines the mass-radius relation for Bose stars, as well as the maximum mass for which a stable equilibrium state exists. When this critical mass is exceeded, the axion star collapses and a number of phenomena can occur when the axions scatter under self-interactions [63][64][65]. Finally, we note that recent studies have shown that there is a second branch of solutions to the GPP equations known as "dense axion stars" in which self-interactions are significant and the full potential of Eq. (1) is needed [40,66]. Whether this state is the result of the collapse process of overcritical dilute axion stars is currently unknown.
In this appendix, we justify Eq.
This shows that df /dt is given by Eq. (19).