Decay of Ultralight Axion Condensates

TL;DR

This paper studies how axion self-interactions affect the stability, mass, and size of galactic-scale condensates, revealing a lower bound on stable fuzzy dark matter structures with lifetimes exceeding the universe's age.

Contribution

It extends previous non-interacting models by including interactions, analyzing stability, and establishing bounds on the size and lifetime of stable axion condensates.

Findings

Self-interactions reduce maximum condensate mass.

Stable condensates have a minimum size bound.

Lifetimes of stable condensates exceed the age of the universe.

Abstract

Axion particles can form macroscopic condensates, whose size can be galactic in scale for models with very small axion masses $m\sim10^{-22}$ eV, and which are sometimes referred to under the name of Fuzzy Dark Matter. Many analyses of these condensates are done in the non-interacting limit, due to the weakness of the self-interaction coupling of axions. We investigate here how certain results change upon inclusion of these interactions, finding a decreased maximum mass and a modified mass-radius relationship. Further, these condensates are, in general, unstable to decay through number-changing interactions. We analyze the stability of galaxy-sized condensates of axion-like particles, and sketch the parameter space of stable configurations as a function of a binding energy parameter. We find a strong lower bound on the size of Fuzzy Dark Matter condensates which are stable to decay, with lifetimes longer than the age of the universe.

Figures

Equations

$${\mathcal{A}}_N\to {\mathcal{A}}_{N-3}+a_p,$$

$$\Phi (r,t)=R(r)e^{-i{E}_{0}t}{a}_0+R(r)e^{i{E}_{0}t}{a}_0^{†}+{\psi }_f(r,t)+{\psi }_f^{†}(r,t),$$

$${\psi }_f(r,t)=\frac{1}{2{\pi }^2}\sum _{l,m}{Y}_l^m(\hat{r}){\int }_0^{\infty }\frac{dpp}{2{\omega }_p}{j}_l(pr)e^{-i{\omega }_{p}t}{a}_{lm}(p),$$

$$a_{lm}(p)=i^{l}p\int d{\Omega }_p{Y}_l^{m\ast }(\hat{p})a(\vec{p}),$$

$$[a_{lm}(p),a_{l^{\prime }{m}^{\prime }}^{†}(p^{\prime })]=2{\omega }_p(2\pi {)}^3\delta (p-p^{\prime }){\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }}.$$

$$\frac{1}{(2\pi {)}^3}\int \frac{d^{3}p}{2{\omega }_p}{e}^{i(\vec{p}\cdot \vec{r}-{\omega }_{p}t)}a(\vec{p}),$$

$$V(\Phi)=m^{2}\,f^{2}\,\Big{[}1-\cos\Big{(}\frac{\Phi}{f}\Big{)}\Big{]},$$

$${\mathcal{M}}_3$$

$$=-i{m}^{2}f\int dtdr{r}^2{J}_3\left[X(r)\right]e^{3i{E}_{0}t}\int d{\Omega }_r⟨0|{\psi }_f(r,t)|\varphi (p)⟩$$

$$\varphi (r,t)=⟨0|{\psi }_f(r,t)|\varphi (p)⟩=\sqrt{4\pi }\frac{e^{-i{\omega }_{p}t}\sin (pr)}{r}.$$

$${\mathcal{M}}_3=-4{\pi }^2\sqrt{4\pi }f\delta (3E_0-{\omega }_p)I_3(p),$$

$$I_3(p)$$

$$\approx \frac{m^2}{48}{\int }_{-\infty }^{\infty }drr{e}^{ipr}X(r{)}^3.$$

$$\Xi (q)=\frac{1}{(2\pi {)}^3}\int d^{3}rX(r)e^{i\vec{q}\cdot \vec{r}}.$$

$${\mathcal{M}}_3=-\frac{{\pi }^2\sqrt{4\pi }{m}^{2}f}{12}\delta (3E_0-{\omega }_p)\int {\delta }^3(\vec{p}-{\vec{q}}_1-{\vec{q}}_2-{\vec{q}}_3)\prod _{k=1}^3\Xi (q_k)d^3{q}_k.$$

$$\frac{A^{\prime }}{A}$$

$$\frac{B^{\prime }}{B}$$

$$X^{\prime \prime }$$

$$d{s}^2=-B(r)d{t}^2+A(r)d{r}^2+r^{2}d{\Omega }^2,$$

$$Y^{\prime \prime }(x)$$

$$a^{\prime }(x)$$

$$b^{\prime }(x)$$

$$Y^{\prime \prime }(x)+\frac{2}{x}{Y}^{\prime }(x)+\frac{1}{8}Y(x{)}^3=0.$$

$$Y(x)=\frac{8\rho }{x^2+{\rho }^2}-\frac{2}{3\rho }-\frac{1}{18{\rho }^3}(x^2+{\rho }^2)+\mathcal{O}([x^2+{\rho }^2{]}^2).$$

$$\Xi(q)=\frac{1}{(2\,\pi)^{2}\,i\,q\,m^{2}\,\Delta}\int_{-\infty}^{\infty}dx\,x\,\exp\Big{(}\frac{i\,q\,x}{m\,\Delta}\Big{)}\,Y(x).$$

$$I_3(p_0)\simeq \frac{32i\pi }{3}\frac{\rho }{\Delta }\exp \left(-\frac{2\sqrt{2}\rho }{\Delta }\right),$$

$$\Gamma_{3}=\frac{1}{T}\int\,\frac{dp}{(2\pi)^{3}2\,\omega_{p}}\,\Big{|}\mathcal{M}_{3}\Big{|}^{2}=\frac{2\,\pi\,f^{2}}{p_{0}}\Big{|}I_{3}(p_{0})\Big{|}^{2},$$

$$\frac{d\tau }{dN}\simeq m\frac{d\tau }{dM}\simeq -\frac{1}{3{\Gamma }_3}.$$

$$\tau =\frac{3y_M{\Delta }^2}{4096{\pi }^3{\rho }^{3}m}\exp \left(\frac{4\sqrt{2}\rho }{\Delta }\right),$$

$$\frac{\lambda_{dB}}{2\pi}=\frac{1}{m\,v}=1.92\text{ kpc}\Big{(}\frac{10^{-22}\text{ eV}}{m}\Big{)}\Big{(}\frac{10\text{ km/sec}}{v}\Big{)},$$

$$M(\kappa )\approx \frac{8.75}{\sqrt{\kappa }}\frac{M_{P}f}{m},R_{99}(\kappa )\approx 1.15\sqrt{\kappa }\frac{M_P}{fm}.$$

$$M(\Delta )\approx 1.75\Delta \frac{M_P{}^2}{m}{R}_{99}(\Delta )\approx \frac{5.75}{m\Delta }$$

$$M{R}_{1/2}\ge 3.925\frac{M_P{}^2}{m^2}\text{(non-interacting bosons)},$$

Full text

Decay of Ultralight Axion Condensates

Joshua Eby,1,2,​ 111Electronic Address: joshaeby@gmail.com   Michael Ma,1,​ 222Electronic Address: mam@ucmail.uc.edu  Peter Suranyi,1,​ 333Electronic Address: peter.suranyi@uc.edu   and L.C.R. Wijewardhana1,​ 444Electronic Address: rohana.wijewardhana@uc.edu

*1**Dept. of Physics, University of Cincinnati, 2600 Clifton Ave, Cincinnati, OH 45221, USA

2Fermi National Accelerator Laboratory, P.O. Box 500, Batavia, IL 60510, USA*

Axion particles can form macroscopic condensates, whose size can be galactic in scale for models with very small axion masses $m\sim 10^{-22}$ eV, and which are sometimes referred to under the name of Fuzzy Dark Matter. Many analyses of these condensates are done in the non-interacting limit, due to the weakness of the self-interaction coupling of axions. We investigate here how certain results change upon inclusion of these interactions, finding a decreased maximum mass and a modified mass-radius relationship. Further, these condensates are, in general, unstable to decay through number-changing interactions. We analyze the stability of galaxy-sized condensates of axion-like particles, and sketch the parameter space of stable configurations as a function of a binding energy parameter. We find a strong lower bound on the size of Fuzzy Dark Matter condensates which are stable to decay, with lifetimes longer than the age of the universe.

1 Introduction

Axions are hypothetical particles introduced to solve the Strong $CP$ problem in QCD [1, 2, 3, 4, 5, 6, 7, 8]. Axions are self-adjoint bosons, with no conserved discrete quantum numbers to guarantee particle number conservation. The axion potential can be written in terms of an angular variable with a $2\pi$ shift symmetry. Axion-like scalar particles also appear in a variety of models beyond QCD, especially in low-energy limits of string theories [9, 10, 11]. Those axions have properties similar to QCD axions, but their mass scales and decay constants are vastly different.

Axion-like particles are among the prime candidates for the composition of dark matter [12, 13, 14, 15, 16, 17]. Axions, being scalar bosons, can condense. Axion condensates have been discussed by numerous groups, with condensate sizes ranging from galaxy or galaxy cluster scale [18, 19], down to stellar size and smaller (termed “axion stars”) [20, 21, 22, 23, 24, 25], all the way to radii of a few meters [26, 27]. Also of interest is the cosmological evolution of the axion field, which has been worked on extensively in [28, 29, 30, 31, 32], but this will not be disussed further here.

Uncondensed QCD axions are not stable, as they can decay to photons through a process $a\to 2\,\gamma$, but the decay rate is slow enough such that most axions created after the Big Bang would survive many Hubble times [33]. Axion condensates, however, may also decay slowly due to the self-interaction of axions [34, 35]. The self-interaction term of their Lagrangian (for both the instanton [36, 37] and chiral [37, 38] cases) have only terms containing an even number of axion fields. Thus, disregarding the rare decay into photons, the axion number is conserved only modulo 2.

In a recent paper [34] we have investigated the decay of weakly bound axion stars due to the self-interaction of axions. The decay proceeds mostly through a sequence of processes,

[TABLE]

where ${\cal{A}}_{k}$ is an axion star, which is a condensate containing $k$ axions, and $a_{p}$ denotes an axion in a scattering state with the magnitude of the momentum $p$. The process (1.1) is the simplest of many possible decay modes responsible for the decay of axion stars.555The decay rates via processes ${\cal{A}}_{N}\to{\cal{A}}_{N-5}+a_{p}$ or ${\cal{A}}_{N}\to{\cal{A}}_{N-4}+a_{p_{1}}+a_{p_{2}}$ are significantly lower and unlikely to have any cosmological significance [34]. This process is allowed by energy-momentum conservation, provided the binding energy of a bound axion is small enough that a relativistic particle can be produced: $\delta E<2\,m\,/\,3$, where $m$ is the mass of a free axion.

In [34] we used an axion field operator, which was the generalization of the field proposed by Ruffini and Bonazzola [39]. To facilitate the decay process, terms creating and annihilating axions in the continuum of scattering states were included in the quantum field of axions, in addition to terms creating and annihilating bound axions [39]. The Ruffini-Bonazzola method is based on taking the expectation value of the quantum Einstein and Klein-Gordon equations in the condensate to derive equations of motion for the metric components and the scalar field. We solved the equations of motions numerically in the weak gravity and weak binding ($\delta E\ll m$) limits, to find the wave function of axions in the condensate [24].

The bound axions are not in momentum eigenstates. They have definite energies, but their wave functions extend over the size of the axion star. Accordingly, the bound axions have an extended momentum distribution as well. In the recent publications [35, 40], the authors questioned the validity of the decay mechanism proposed in [34], arguing that momentum is not conserved in (1.1), and that the decay rate through this process is exactly zero by the Optical Theorem. These authors have further suggested that one can show the rate to be zero by the classical equation of motion for the condensate. We will address these issues and explain our response in Appendix A.

In the present paper we will apply our method of discussing the decay of dilute axion stars [34] to condensates of cosmological size; such models have been referred to previously as Fuzzy Dark Matter (FDM) [41]. Condensates of galactic sizes have been considered by a number of authors, and typically correspond to a scalar particle mass of $m\sim 10^{-22}-10^{-21}$ eV [42, 43, 44, 41, 45, 46, 47, 48, 49, 50, 19].666For a brief but recent review of ultralight scalar field dark matter models, see [51] and references therein. When such condensates are formed from real scalars, a version of the decay analysis of [34] applies, and we will investigate whether interesting bounds can be placed on these models by taking decays into account. As we will explain in the next section, we will utilize the axion potential with a cosine dependence on the field; other proposals, for example a $\cosh$ potential [52, 53], have been investigated in the context of ultralight scalars as well.

It is also an aim of this paper to emphasize the inclusion of axion self-interactions in investigations of axion condensates. Although the self coupling $\lambda\sim m^{2}/f^{2}\sim 10^{-95}\lll 1$ for typical FDM models ($f$ is the axion decay constant), the astronomical number of axions in a condensate participating in these interactions could lead to large corrections to certain important physical quanties. We investigate the macroscopic properties of these condensates using the fully self-interacting analysis and emphasize the differences from the non-interacting limit.

In the next section, we give a more detailed explanation of how axion star decay through the process (1.1) can be calculated. In Section 3, we will outline the calculation of the wavefunction, following largely [24]. In Section 4, we apply the formulas for the macroscopic properties and decay rates to condensates formed from ultralight axion-like particles. We conclude in Section 5.

We will use natural units throughout, where $\hbar=c=1$.

2 Decay through self interaction

There is a variety of methods for the quantitative investigation of axion condensates, as classical [33, 21, 20, 54, 55], quantum mechanical [56, 57, 25, 58, 59, 60, 61], and field theoretic [26, 24, 27]. Our field theoretic discussion of the decay of an axion condensates into relativistic axions [34] was based on an extension of the Ruffini-Bonazzola operator [39], by the addition of scattering state contributions.777Appendix B we will discuss why a continuous spectrum of scattering state solutions can be added to the boson field operator. Furthermore, we will also discuss why using free spherical wave scattering states is quite sufficient in our calculations. Thus, we proposed to extend the expansion of the boson field using the form [34]

[TABLE]

where $E_{0}$ is the energy eigenvalue of a single bound axion, and where $R(r)$ and $a_{0}$ are the wave function and annihilation operator of the axions in the condensate (respectively). $\psi_{f}(r,t)$ is the annihilation part of a complete system of free axion operators expanded in scattering states,

[TABLE]

where $\omega_{p}$ and $a_{lm}(p)$ are the energy eigenvalue and the annihilation operator of the scattering state axion, with quantum numbers $l$ and $m$, respectively. The functions $j_{l}(x)$ and $Y_{l}^{m}(\hat{x})$ are spherical Bessel functions and spherical harmonics, respectively.

The annihilation operator in the spherical wave basis is defined by

[TABLE]

where $a(\vec{p})$ is the annihilation operator for a particle which is the eigenstate of the momentum operator with eigenvalue $\vec{p}$ (i.e. a plane wave). This annihilation operator and its adjoint, the creation operator, satisfy the commutation relation

[TABLE]

Note that (2.2) is exactly equal to the negative frequency part of a complete system of free axion states given by

[TABLE]

which was used in [34] to investigate the decay of QCD axion stars.

We will see later that the bosons produced by the decay of a weakly bound boson condensate are relativistic. Therefore, we have chosen to use free particle states to approximate the scattering states. For the purposes of this calculation, this choice is admissible, because the energy level of produced axions is sufficiently high compared to the effective depth of the potential created by gravitation and self-interactions. This is explained in greater detail in Appendix B. In a future work we will take into account corrections to this approximation.

Let us consider now (1.1) in the Born approximation. The axion self-interaction potential can be approximated by the so-called instanton potential [36, 37]

[TABLE]

where $m$ and $f$ are the axion mass and decay constant (respectively). In the Ruffini-Bonazzola paradigm, one finds the expectation value of eq. (2.6) transforms the cosine into a Bessel function $J_{0}$ [24].888More generally, the annihilation process for $k$ bound axions generates an effective potential proportional to the Bessel function $J_{k}$, as explained in the Appendix of [34]. One also finds that the transition matrix element for the process (1.1) is

[TABLE]

where $X(r)=2\,\sqrt{N}\,R(r)\,/\,f$ is the rescaled wave function of the condensate which, as we will see in Section 3, can be obtained by solving the equations of motion [24]. We are considering transitions of the form (1.1), where $\langle N|$ is the initial $N$ particle condensate (the left hand side of (1.1)), and $|N-3,\,\phi(p)\rangle$ is the direct product of the final state $N-3$ particle condensate and a scattering state $\phi(p)$ of momentum magnitude $p$ (the right hand side of (1.1)).

We restrict this work to non-rotating axion condensates; the reason for this is twofold. First, for the potential in eq. (2.6) and the parameters we use here, only the inner cores of galaxies composed of axion particles can be described as a condensate; outside of this inner region, the dark matter halo is described by a virialized gas of particles [19] and cannot participate in the decay processes we describe here.999For ultralight bosons with repulsive self-interaction, like those presented in e.g. [45, 48], the condensates can be much larger, and can even constitute the entire dark matter halo. Because the condensed core is small compared to the full radius of the halo, it carries at most a tiny fraction of the angular momentum of the galaxy, so as a first approximation we believe restricting to $\ell=0$ angular momentum states is appropriate. The second reason is that a full treatment of rotating axion condensates has not yet been done, though slowly rotating condensates were analyzed in a particular limit by [14]. This is a topic we hope to return to in the near future.

Because we work in the limit of zero angular momentum, annihilation processes of the form $a\,a\to G$, where $G$ is a spin-2 graviton [62, 63], have a rate of zero. Such an interaction would require the participating axions to have at least $\hbar$ of angular momentum each; and even if we accounted for the nonzero rotation speed of the galaxy, by our estimation the vast majority of axions in a galactic condensate would have far less angular momentum than what would be necessary for this process to occur.

For static condensates, note that the integration over $\Omega_{r}$ in eq. (2) vanishes for all but $s$-wave axions.101010Should we consider rotating axion stars, higher angular momentum scattering states would also contribute. In that case, the scattering state axion is described by the zero angular momentum contribution only, $|\phi(p)\rangle=a_{00}^{\dagger}(p)|0\rangle$. Then the wave function of the emitted axion is

[TABLE]

The integration over $t$ also fixes the energy of the outgoing axion to $\omega_{p}=3\,E_{0}$. The matrix element takes the form [34]

[TABLE]

where the dimensionless integral is

[TABLE]

The symmetry of the integrand for the substitution $r\to-r$ has been used to extend the integration region to all real values of $r$, and to switch from $\sin(p\,r)$ to $e^{i\,p\,r}$ in the integrand. In the second equality, we expanded the Bessel function $J_{3}$ to leading order, an appropriate limit for dilute axion stars.

Now observe that for dilute axion stars the radius of the star $R$ is very large. In other words, $X(r)$ has a large coordinate uncertainty, $\delta r\sim R\sim(m\,\Delta)^{-1}$, where $\Delta=\sqrt{1-(E_{0}\,/\,m)^{2}}\ll 1$ [24, 34]. As a result, the range of $p$, as represented by the momentum uncertainty $\delta p\sim m\,\Delta\ll m$, is very small. Then due to the delta function in eq. (2.9), enforcing energy conservation, the emitted axion has a momentum peaked at a very large value, $p\simeq\sqrt{8}\,m$; as a result, the matrix element (2.9) is very small for weak binding. However, as the binding energy $\delta E$ increases, ${\cal M}_{3}$ will take larger values and the decay rate $\Gamma\sim|{\cal M}_{3}|^{2}$ also increases.

Now to bring out issues related to momentum conservation, we can define the momentum representation wave function as

[TABLE]

Then we can rewrite eq. (2.9) as

[TABLE]

Since for weakly bound condensates $p\simeq\sqrt{8}\,m$, the magnitude of the transition matrix depends crucially on the large $q$ tail of momentum distribution $\Xi(q)$. However, for calculational purposes, it is still advantageous to use (2.9) rather than (2.12). One can rely on the numerical solution of the equations of motion, as explained below, using a simple approach to estimate approximate behavior of $\Xi(q)$ at large $q$ [34].

Note that the process (1.1) is not the only possible channel through which decay can proceed; however, it is by far the dominant process. First, we have shown previously [34] that processes of the form ${\cal A}_{N}\to{\cal A}_{N-(2\,j+1)}+a_{p},$ are suppressed by higher powers of $\Delta$ for each higher $j>1$. In the weak-binding limit, where $\Delta\ll 1$, these corrections are completely negligible. On the other hand, this argument breaks down for dense axion stars [27, 64], where $\Delta=\mathcal{O}(1)$; we will return to this case in a future publication.

Second, there are processes of the form ${\cal A}_{N}\to{\cal A}_{N-k}+\mu\,a_{p},$ where $\mu>1$ axions are emitted at once. Unlike the process (1.1), the emission of $\mu>1$ axions from a condensate can proceed on-shell, implying that the corresponding decay rate has a weak dependence on $\Delta$. Nonetheless, as shown in [34], these processes are suppressed by the very small factor $m^{2}/f^{2}$ for each additional axion in the final state. Since in FDM $m^{2}/f^{2}\sim 10^{-95}\lll 1$, we can safely neglect these corrections as well. We conclude that (1.1) is by far the dominant contribution to the decay of axion condensates.

3 The calculation of the wave function of the condensate

We review here the calculation of the condensate wavefunction $X(r)$ [24]. The matrix elements of the $rr$ and $tt$ components of the Einstein equation, along with the Klein-Gordon equation, form a closed set of equations for the metric and the axion field, $X(r)$:

[TABLE]

where the metric is

[TABLE]

with $\delta=f^{2}/M_{P}{}^{2}$ and $M_{P}=1\,/\,\sqrt{G}=1.22\times 10^{19}$ GeV (the Planck mass). As above, we have taken only the leading contribution to the Bessel function which represents the self-interaction potential; doing so preserves the leading attractive $X(r)^{4}$ interaction term.

Assuming that $\delta\ll 1$, a condition satisfied in applications where gravity is weak (Newtonian limit), we can write $A=1+\delta\,a$ and $B=1+\delta\,b$, where $a,\,b=\mathcal{O}(1)$. Furthermore using the large radius approximation and the definition of the scale parameter $\Delta=\sqrt{m^{2}-E_{0}{}^{2}}/m$, we can introduce dimensionless radial coordinate as $x=m\,r\,\Delta$. In that case the axion field also scales with its engineering dimension, such that $X(r)=\Delta\,Y(x)$, leading to the following system of equations for $a$, $b$, and $Y$ in leading order of $\Delta$ and $\delta$ [24]:

[TABLE]

where111111This definition of $\kappa$ appears to differ by a factor of $8\pi$ compared with [34], because in that work we wrote $\delta$ in terms of the reduced Planck mass. In fact, the two definitions of $\kappa$ are equivalent. $\kappa=8\pi\delta\,/\Delta^{2}$. Since $b(x)$ is proportional to the Newtonian gravitational potential, $\kappa\sim G$ is the effective coupling constant of the field $Y(x)$ to gravity. Further details, and a more comprehensive justification of this double expansion of the equations, can be found in [24].

Solutions of the equations of motion (3) correspond to ground-state configurations of axions, which can be stable or metastable. In [24], we solved these equations and found a spectrum of solutions which were parameterized by $\kappa$ (or, equivalently, by $\Delta$). When applied to QCD axion parameters $m=10^{-5}$ eV and $f=6\times 10^{11}$ GeV, we found a maximum gravitationally stable mass of $M_{c}\sim 10^{19}$ kg.121212The maximum masses for attractive interactions were discussed by Stoof [68] in the context of condensed matter BECs, and in the context of boson stars by Chavanis and Delfini [56, 57]. By rescaling these solutions to values of $m$ and $f$ appropriate for FDM, we can analyze the properties of galactic-size condensates in a way that includes the self-interaction term in the potential. The physical interpretation of these condensates is that they form the cores of FDM halos; they are surrounded by a virialized distribution of axions which extend to the outer edge of the dark matter halo.

To analyze the decay of these condensates through processes like (1.1), we investigate the singularity structure of solutions of (3). Now, (3) is a system of equations with two singular points, $x=0$ and $x=\infty$. Using boundary conditions we require that the solutions are regular at $x=0$ and decrease exponentially at $x\to\infty$. The solutions are even functions of $x$, so they also approach zero at $x\to-\infty$. Then they are real analytic functions at $-\infty<x<\infty$ and can be continued into the complex plane of $x$. As they fast vanish at infinity, the contour of integration can be moved up along the imaginary axis in the rescaled version of the integrals in eqs. (2.9) and (2), until we reach a singularity in the complex plane. The contribution of that singularity dominates the decay rate integrals at large momentum.

It is easy to show that the Klein-Gordon equation (3), in which the the leading order singular terms are retained, is131313In fact, this expression contains the next to leading order term proportional to $Y^{\prime}(x)$.

[TABLE]

Near the singular point $x=i\,\rho$, this has a solution of the form

[TABLE]

The parameter $\rho$ is an integration constant, having a one-to-one relationship with the rescaled central density of the axion field $Y(0)^{2}$, and in turn, with mass and the radius of the condensate. In fact, high order Taylor series expansion of equations around $x=0$ show that the singularity closest to the origin is indeed of the form (3.7), connecting the value of $Y(0)$ with $\rho$ [34]. In principle, gravitational interactions have an effect on the solutions $Y(x)$. In practice, however, the term in the equations of motion (3) which couple $Y(x)$ to gravity give a subleading contribution to the singularity. We can therefore solve the equations in the limit $\kappa\ll 1$, i.e. where gravity decouples. In that limit, the nontrivial solution has $Y(0)=12.268$, which implies a fixed value $\rho=.603156$.

Finally, we can rewrite the Fourier transform of eq. (2.11) as

[TABLE]

In this form, it is clear that at small $\Delta$ the singular term of (3.7) term dominates the integral. To calculate the decay rate, we follow the procedure of [34]: we take the leading order solution for $Y(x)$ near the singularity $x=i\,\rho$, given by eq. (3.7), and evaluate $I_{3}(p)$ in the matrix element of eq. (2.9). The result is

[TABLE]

where $p_{0}=\sqrt{9\,E_{0}{}^{2}-m^{2}}\simeq\sqrt{8}\,m$.

The decay rate for the process (1.1) is then

[TABLE]

where $T$ is the duration of the decay process. Then the lifetime of the condensate through this decay process is

[TABLE]

Further details regarding the evaluation of eq. (3.11) can be found in [34]; the result for the process (1.1) is

[TABLE]

where $y_{M}\simeq 75.4$ is determined by the relationship between $M$ and $\Delta$ in the large $\Delta$ region [24]. The lifetime is a monotonically decreasing function of $\Delta$ in the relevant range; in the case of QCD axions, we found in [34] that above a value $\Delta\simeq.05-.06$, axion condensates become very unstable to decay to relativistic axions, their lifetimes becoming shorter than the age of the universe. We will examine the consequences of this fact in the context of ultralight axions in the next section.

4 Stable spectrum of ultralight axion condensates

Very light axion fields can have de Broglie wavelengths as large as entire dark matter halos, possibly implying a connection between these two scales. Ultralight bosons appear often in theories of physics beyond the Standard Model, including those requiring compactification of extra dimensions. Such models, termed “Fuzzy Dark Matter” (FDM) [41], have been written about extensively [42, 43, 44, 41, 45, 46, 47, 48, 49, 50, 19]. While there are significant constraints on these models141414While this work was being finalized, a paper appeared suggesting a strong tension between the preferred mass scale for FDM, $m\sim 10^{-22}-10^{-21}$ eV, and data from Lyman-$\alpha$ forest simulations [65]. We will not comment here about whether such constraints could rule out FDM as a viable paradigm., they remain a viable alternative to WIMP or QCD axion models of dark matter.

We will consider an ultralight axion in this context, using the potential of eq. (2.6). The mass of the ultralight axion in question will be taken to be $m\sim 10^{-22}$ eV, which gives the right approximate scale for the size of dark matter halos [19],

[TABLE]

where $v$ is the velocity in the halo. This choice is also consistent with the known epoch of matter-radiation equality. Further, a decay constant of $f\sim 5\times 10^{16}$ GeV naturally leads to the correct relic density, and can thus account for the observed dark matter abundance [19]; however, to remain as general as possible, we allow $f$ to deviate from this value by a few orders of magnitude. At the upper limit of what we consider, $f=10^{18}$ GeV is still below the Planck scale, implying that the parameter $\delta=f^{2}/M_{P}^{2}\approx.007\ll 1$; thus, the weak-gravity approximation holds reasonably well over our entire range.

In [24], we found numerically the solutions to the system (3) over a wide range of $\kappa$. In the FDM picture, these solutions correspond to the cores of FDM halos discussed (most recently) in [19]. We found that there exists a maximum mass at $\kappa\approx.34$, above which axion condensates are gravitationally unstable. On the stable branch of masses $\kappa>.34$, the mass $M$ and radius151515We use the common convention that $R_{99}$, the radius inside which $.99$ of the mass of the condensate is located, represents the “size” of the condensate. $R_{99}$ of the condensate are fit by the functions [24]

[TABLE]

We observe in Figure 1 that at fixed $m$, the value of $f$ determines the position of the maximum mass, and thus the turnaround of the function $M(\Delta)$. A similar structure can be observed in Figure 2 for the radius, where the position of the maximum mass corresponds to a slight dip in the otherwise straight line representing $R_{99}(\Delta)$. Trading $\kappa$ for $\Delta$ in eq. (4.2), we see that the lines

[TABLE]

bound the full set of solutions from above.

For the specific choice of FDM parameters $m=10^{-22}$ eV and $f=5\times 10^{16}$ GeV we find the maximum mass $M_{c}\approx 6\times 10^{10}M_{\odot}$; this is lower than the value found in the non-interacting limit of $8\times 10^{11}M_{\odot}$ [19] by about an order of magnitude, due to the inclusion of attractive self-interactions. Our estimate of the maximum mass also agrees well with the recent analysis of [66], which also includes the leading attractive self-interaction. More generally, the mass and radius of FDM halo cores over a wide range of the scale parameter $\Delta$ and at different values of $f$ are illustrated in Figures 1 and 2.

We can also analyze the relationship between $M$ and $R_{99}$, which were investigated for both attractive and repulsive self-interactions in [56, 57]. In a recent paper [19], the authors present a bound on the product

[TABLE]

where $R_{1/2}$ is the radius inside which $.5$ of the mass of the condensate is contained; the inequality is saturated for stationary, ground state configurations, i.e. for the condensates considered here. In our calculation, on the stable branch of solutions (where $\kappa>.34$), we find the product

[TABLE]

using eq. (4.2). To find the relationship between $R_{99}$ and $R_{1/2}$, we calculated their ratio numerically and found that $R_{99}/R_{1/2}\approx 3.65$ holds within 1%, in a range of $1/\Delta$ extending over many orders of magnitude. This implies

[TABLE]

This product is below the lower bound (4.4) presented in [19] due to our inclusion of self-interactions.

It is worth noting also that in the limit $f\sim M_{P}$, the mass-radius relationship for axion condensates approaches that of a black hole. This is easy to see using eq. (4.2):

[TABLE]

Near $\kappa=\mathcal{O}(1)$ (the position of the maximum mass) and $f\sim M_{P}$, we find $G\,M/R_{99}\sim 1$, implying that $R_{99}\sim R_{S}$, the Schwarzschild radius.

We must also ensure that the weak-binding approximation, on which our analysis [34] and the classical one of [19] depends, is also valid.161616A stability analysis for very strongly-bound condensates, with $\Delta=\mathcal{O}(1)$, is a task we plan to undertake in the near future. We observe in Figure 2 that cores of radius $R\lesssim 1$ pc have $\Delta\gtrsim.3$, and become relatively strongly bound. Such cores would not be well-described by our weak-binding analysis.

An estimate of the decay rate, obtained from the expression derived in [34], is given in eq. (3.12); it is a one-to-one function of the binding energy parameter $\Delta$ in the region of interest. Because the condensate mass $M$ is determined by the value of $\Delta$, it is easy to connect $\tau$ to $M$ as well. Following the analysis of [34], we find that axion condensates with $m\sim 10^{-22}$ eV which have $\Delta\gtrsim.1$ have lifetimes shorter than the age of the universe. This region is represented by the shaded regions in Figures 1 and 2. For $f\lesssim 10^{18}$ GeV, the transition to decay instability occurs on the gravitationally unstable branch of solutions; however at $f\gtrsim 10^{18}$ GeV, the bounds from decay are as strong or stronger than those coming from gravitational stability. This can be an important constraint on bound structures originating in theories of Planck scale axions.

In Figure 2, it is particularly striking that almost regardless of the value of $f$, condensates with $R\lesssim 2$ pc lie in the unstable, shaded region. This implies a fundamental limiting radius of $R_{min}\sim 2$ pc for FDM cores composed of axions. Observe also that, in Figure 1, it is easy to read off the maximum mass of FDM condensates for each value of the decay constant $f$. For any axion theory with $f\ll M_{P}$, no stable condensate exists with a mass larger than about $M_{max}\approx 10^{12}M_{\odot}$.

5 Conclusions

In a previous publication [24] we established scaling relations for the mass and radius of weakly bound condensates of interacting axions, as functions of the mass of the axion, its decay constant, and the particle energy (or alternatively the central density). We also found the maximum mass and size of the bound states as functions of those parameters. In this paper we have applied those results to condensates of axions forming FDM, providing corrections to similar calculations which neglect the self-interaction of axions [44, 41, 45, 47, 49, 50, 19].

In another publication [34] we developed methods to calculate the lifetime of axion condensates due to their self-interaction, through the four-axion interaction term in which three bound axions produce a single free relativistic axion. Here we have applied those results to estimate the lifetime of condensates formed from FDM. We have found that, provided the decay constant of FDM axions satisfies $f\lesssim 0.05\,M_{P}$, all condensates having binding energy smaller than that of those of maximal mass have lifetimes greater than the age of the universe making them viable candidates for forming central regions of galactic halos. We have also explained in details the decay mechanism described in [34] and further clarified the justification of its validity.

The methods we have described here, based on previous work in [24, 34], rely on a double expansion to leading order in $\delta$ and $\Delta$. This is appropriate for so-called dilute axion stars, which are weakly bound. However, it has been pointed out that an effective short-distance repulsive interaction in the axion potential also gives rise to dense axion stars [27, 64, 73], which are at least energetically stable. We plan to extend our methods to this regime to analyze the properties of these states in the near future.

Collapsing boson condensates have been investigated by a number of groups [67, 68, 69, 70, 64, 71, 72, 73]. Recently, we found that supercritical QCD axion condensates, having masses larger than the maximal allowed stable mass, collapse towards the global minimum of the effective axion potential [64, 73]. Similar arguments indicate that FDM axion condensates which exceed the maximum mass $M_{c}$ will also collapse in this way. Such supercritical condensates can form during galactic collisions, in a manner similar to the mechanism outlined in [74]; such events would lead to collapse, causing the condensate to emit a large number of relativistic particles. Consequences of such events will be studied in a future publication.

Acknowledgements

We thank M. Amin, P. Argyres, E. Braaten, J. Brod, P. Fox, R. Harnik, A. Kagan, A. Mohapatra, K. Schutz, G. Semenoff, M. Takimoto, H. Zhang, and J. Zupan for fruitful discussions. The work of J.E. was supported by the U.S. Department of Energy, Office of Science, Office of Workforce Development for Teachers and Scientists, Office of Science Graduate Student Research (SCGSR) program. The SCGSR program is administered by the Oak Ridge Institute for Science and Education for the DOE under contract number DE-SC0014664. J.E. also thanks the Fermilab Theory Group and the Weizmann Institute Department of Physics for their hospitality. L.C.R.W. thanks Mainz Institute for Theoretical Physics for their hospitality, and the participants of the Quantum Vacuum and Gravitation program, especially M. Bartelmann, A. Mazumdar, and T. Prokopec, for discussions.

Appendix A: The ${\cal A}_{N}\to{\cal A}_{N-3}+a_{p}$ decay

The results of our paper, ”The Lifetime of Axion Stars” [34], have been called into question in recent publications [35, 40]. Before discussing this issue in detail, we review the premises of our work. We proposed a way to discuss the decay of axion stars through the repeated elementary decay mode

[TABLE]

where ${\cal A}_{k}$ denotes an axion star, which is a condensate of $k$ axions and $a_{p}$ denotes a (relativistic) free axion labeled by its momentum $p$. This calculation was performed using an extended axion field operator, eq. (2.1), which included both bound and scattering state contributions.

First of all, we need to establish the fact that there is no conservation law that would forbid (A.1). Axions are real bosons, and consequently the axion number is not conserved. Axions, being coupled to the electromagnetic field as $\Phi\,E\cdot B$, can decay to photons through the slow process $a\to 2\,\gamma$, though this decay process does not significantly affect the lifetime of axion stars [33]. Disregarding axion decay to photons, the axion number is conserved modulo 2, as the self-interaction terms of the axion potential contains only even powers of the axion field, through a dependence of $\cos(\Phi\,/\,f)$ in eq. (2.6). Momentum and energy conservation would allow the decay process to proceed even if the condensates were in momentum eigenstates: a decay process is always allowed if the sum of the masses of the decay products is smaller than the mass of the decaying object, unless the conservation of discrete quantum numbers prevents it.

At any rate, the condensates represented in (A.1) are not in momentum eigenstates. They are quantum objects, which have extended wave functions localized on a large radius $R$. Consequently, though they have mean momentum of zero, their momentum distributions are smeared. In fact, the momentum distributions extend from zero to infinity, albeit with fast decreasing amplitudes.

We have to emphasize that a Gross-Pitäevskii approach, being the non-relativistic limit of the relativistic quantum field theory, cannot be used to discuss particle number violating processes. In the non-relativistic limit the axion number is conserved and process (A.1) is forbidden. Note however that a method to adapt the Gross-Pitäevskii approach to these processes was formulated by the authors of [75]. We also point out Gell-Mann’s totalitarian principle [76], borrowed from T. H. White’s ”The Once and Future King”, which when it is applied to physics is: ”Everything not forbidden is compulsory.” This principle implies that the process (A.1) should exist. In the present context, this principle also implies that taking the nonrelativistic limit discards important contributions that make such a transition possible.

Now the questions raised in [35, 40] have two sides: (a) microscopic, and (b) macroscopic. The microscopic argument (a) pertains to the question of momentum conservation in the elementary process of three bound axions at rest turn into a moving axion. The macroscopic argument (b) concerns the whole condensate, namely how momentum is conserved in the overall process (A.1), i.e. how the momentum carried away by the produced scattering state axion is transferred to the outgoing axion star, ${\cal A}_{N-3}$.

The argument of [35, 40] pertaining to the microscopic side (a) is that in the overwhelmingly dominating Born approximation, the process around which the decay process (1.1) is built is

[TABLE]

where $a_{c}$ represents a bound axion and $a_{f}$ represents the final state axion, is inadmissible due to energy-momentum conservation. The total energy of the three bound axions is $E_{tot}\simeq 3\,m$, which is certainly sufficient to produce a free relativistic axion. However, if the three axions were in momentum representation, having zero momentum in the rest system of the axion star, then the axion in the final state would not have the required momentum of $p\simeq\sqrt{8}m$.

To resolve this problem, observe that the three axions contributing to the decay process are not at rest. The notion of ”rest” is tied to particles in momentum representation, and for the axions in the axion star, this is not so. They have momentum space wave functions, which, though peaked near zero, extend to large momenta (albeit with very small probability). As we discussed above, the uncertainty of the momentum of each axion is $\delta p\sim R^{-1}$. The probability that three of the bound axions have sufficient total momentum to create a free axion, $p=\sqrt{9\,E_{0}{}^{2}-m^{2}}$, decreases rapidly with the size of a condensate, though it is never exactly zero. Consequently, local momentum conservation always allows the decay of condensate. The question of whether this decay process affects cosmology, due to the survival or non-survival of the condensate, is a question of numerical calculations and depends on the parameters of the axion theory, as well as the size and mass of the condensate.

Another argument presented in [35, 40] to the local momentum conservation is based on the optical theorem. The argument invokes diagram

[TABLE]

claiming that the propagator of the axion $a_{f}$ in the intermediate state does not have an imaginary part, because presumably its denominator, $E_{tot}^{2}-p^{2}-m^{2}\simeq 8\,m^{2}\not=0$. The diagram having no imaginary part, the decay rate must vanish. However, this argument is based again on the premise that the momentum of the condensed axions in the initial and final states of the process (A.3) is zero. This is not a valid assumption, since those particles are not in momentum eigenstates and with a tiny probability they can produce sufficient momentum to allow the axion in the intermediate state to go on mass-shell. A comprehensive discussion of how the imaginary part appears in a process like (A.3) for particles in a condensate can be found in [75].

The macroscopic argument (b) claims that even if the above discussed elementary process is possible, there is no mechanism by which the three axions participating in the decay process transmit the momentum to the axion star as a whole; so (the argument goes) contrary to the requirement of momentum conservation in (A.1) the recoil of the axion star ${\cal A}_{N-3}$ in the final state of (A.1) is not possible.

Consider, however, that the condensate itself is not in momentum representation either, unlike condensates in a container, for which the boundary conditions are set by the wall of the container. Still, one could evoke the valid argument that just like the expectation value of the coordinate, the expectation value of the momentum of the condensate must be zero. Consequently, the only constraint we can impose on the decay process is the conservation of the expectation value of the momentum. Now, the condensate of the final state of (A.1) is not in momentum representation either, though its average momentum is zero as well. Consider now the created scattering state axion. It is produced as a zero angular momentum spherical wave, going with the same probability into every direction. Though the magnitude of the momentum is sharply peaked at a particular value, the spherical wave also has a vanishing average momentum. Thus, the average momentum is conserved in (A.1).

One should not confuse the emission of the axion with its detection. Suppose we detect an emitted axion, the decay product of (A.1). By performing a measurement we alter the system. Just like in the case of the collapse of a simple wave packet by performing a measurement, the conservation of the average momentum is valid only if we include the measuring device, which absorbs the appropriate amount of momentum, to make the average momentum of the complete system zero.

We turn now to the last critique presented in [35, 40], namely that the classical equation of motion for the axion star precludes a linear coupling to a scattering state axion. In [35], the authors write about our previous work [34] as follows:

They expanded the scalar axion field $\phi$ around the classical field $\phi_{0}$: $\phi=\phi_{0}+\tilde{\phi}$, where $\tilde{\phi}$ is the quantum fluctuation field. The Hamiltonian includes a term proportional to $\phi_{0}{}^{3}\,\tilde{\phi}$ from the axion interaction potential. They claimed that this term produces transitions of $3$ condensate axions into one relativistic axion of energy $3\,m_{a}$. However, the sum of all terms in the Hamiltonian that are linear in $\tilde{\phi}$ is zero by the classical equations of motion for $\phi_{0}$. Thus the matrix element for producing a single relativistic axion is [math].

On the contrary, we have emphasized (and explain further in Appendix B) that a complete set of states satisfying the Klein-Gordon equation for the axion field contains both bound and scattering state components, $\Phi=\Phi_{b}+\Phi_{s}$. Our expansion does not distinguish a classical and a quantum component, and as such, one cannot identify $\Phi_{b}$ as $\phi_{0}$, a purely classical field. By direct calculation we see that the matrix element for (A.1) is nonzero and is proportional to $\Phi_{s}\langle N|\Phi_{b}{}^{3}|N-3\rangle$.

We can see this more quantitatively by analyzing the equation of motion. The standard procedure to calculate the equation of motion is to take the matrix element

[TABLE]

of the Klein-Gordon operator $KG[\Phi]$. This is precisely what we have used in eq. (3.3) and is equivalent to the Gross-Pitaevskii equation used by numerous other authors, e.g. [25, 27, 56]. In an “exact” field theory, any matrix element of the form $\langle\psi|KG[\Phi]|\psi^{\prime}\rangle=0$, including

[TABLE]

However, because this exact theory is not known, one is restricted to using some ansatz for the field $\Phi$. We have chosen the Ruffini-Bonazzola ansatz, expanded to include the scattering state solutions, which was given in eq. (2.1). With this choice, the wavefunction $R(r)$ (which is equivalent to the Gross-Pitaevskii wavefunction), does not satisfy eq. (A.5), and in this same parameterization, the matrix element for the process (A.1) is nonzero as well.

We are working to extend the Ruffini-Bonazzola paradigm so that eq. (A.5), as well as higher-order expressions in the operator Klein-Gordon equation, can be simultaneously satisfied. In this extension, we find that the rate for the process (A.1) is still nonzero, and is equal at leading order to the result we obtain here. We will present these and related findings in a future publication.

Appendix B: The continuous spectrum in the Ruffini-Bonazzola equations

The most general bound solution, discussed in [39], for a real scalar field in spherically symmetric metric of eq. (3.4) and satisfying the non-interacting Klein-Gordon equation is

[TABLE]

where $c_{nlm}$ are arbitrary amplitudes of the $n$th bounds state solution with angular momentum $l$. The wave functions $R_{nl}(r)$ satisfy wave equations

[TABLE]

To be able to describe condensates of bosons, Ruffini and Bonazzola [39] introduced second quantization by promoting coefficients $c_{nlm}$ to creation and annihilation operators $c_{nlm}\to a_{nlm}$ and $c_{nlm}^{*}\to a_{nlm}^{\dagger}$; these operators satisfy

[TABLE]

Now notice that $\Phi_{b}$ is not the most general solution of the wave equation. Because the attractive gravitational potential, which makes bound states solutions possible, vanishes at large $r$, (B.7) has scattering state solutions as well. After quantization these states can be written as

[TABLE]

where $k=\sqrt{\omega_{k}^{2}-m^{2}}$ is the momentum and

[TABLE]

The complete set of solutions to the Klein-Gordon equation is $\Phi=\Phi_{b}+\Phi_{s}$.

The Ruffini-Bonazzola method was used by Barranco and Bernal [26] and by us [24] to describe axion stars. It is sufficient to use of quantized field $\Phi_{b}$ in leading order to describe static axion stars. However, the term $\Phi_{s}$ becomes significant if we notice that in first order of the expansion of the axion potential using $\Phi$ results in an operator

[TABLE]

This has been described in detail in [34] and also applied in the present paper.

The second subject we discuss in this Appendix is the reason why free spherical waves, $f_{l}(k)\to j_{l}(k\,r)$, can be used to calculate the decay rate of axion stars. This is the consequence of the fact that dilute axion stars, those we consider in this work, are weakly bound. In fact, only axion stars whose particle energy satisfies $|E_{0}-m|\lesssim 0.002\,m$, corresponding to a value of $\Delta\approx.05$, could survive from the big bang until the present epoch.

To further justify our use of free spherical waves in the final state, consider (3.3) and (3). Using $Y(x)=X(r)/\Delta$ and $x=m\,r\,\Delta$, where $\Delta^{2}\simeq-2\,(E_{0}-m)/m$, the Klein-Gordon equation in (3) can be rewritten as the nonlinear Schrödinger equation

[TABLE]

where we defined an effective potential

[TABLE]

Now considering that $b(r)<0$, $b^{\prime}(r)>0$ and $X^{\prime}(r)<0$ over the whole range of $r$, as shown by the numerical calculations [24], the effective potential $V_{\rm eff}$ satisfies

[TABLE]

The quantity $3\,(E_{0}-m)\gtrsim-0.006\,m$ at $\Delta\lesssim.05$. This shows that the depth of the effective potential is of $\mathcal{O}(10^{-3}\,m)$ or smaller, implying that produced relativistic axions, which have energy $E\simeq 3\,m$, can well be regarded as free. The situation would be different in a discussion of dense axion star states, where the binding energy is much more significant; we will return to this topic in a future work.