Warm and Fuzzy Dark Matter: Free Streaming of Wave Dark Matter

Abstract

Wave or fuzzy dark matter that is produced with relativistic wavenumbers exhibits free streaming effects analogous to warm or hot particle dark matter with relativistic momenta. Axions produced after inflation provide such a warm or mildly relativistic candidate, where the enhanced suppression and observational bounds are only moderately stronger than that from wave propagation of initially cold axions. More generally, the free streaming damping also impacts isocurvature fluctuations from generation in causally disconnected patches. As coherent spatial fluctuations free stream away they leave incoherent and transient superpositions in their wakes. These multiple wave momentum streams are the wave analogue of particle phase space fluctuations or directional collisionless damping of massive neutrinos or hot dark matter. The observable impact on both adiabatic and isocurvature fluctuations of fuzzy dark matter can differ from their cold dark matter counterparts due to free streaming depending on how warm or hot is their momentum distribution.

I Introduction

Wave dark matter refers to bosonic dark matter with masses $m\lesssim 30$ eV such that the occupation number is much greater than one, and can arise in a variety of theoretical contexts (see e.g. [1] for a review). One of the leading candidates is axion dark matter, where dark matter behaves as a classical wave below the de Broglie scale. The axion, originally proposed to explain upper limits on the neutron electric dipole moment and solve the strong CP problem [2, 3, 4, 5, 6], is also a viable dark matter candidate [7, 8, 9] and has stimulated much interest in dark matter physics and numerical experiment searches [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38]. The mass spectrum of axions extends beyond the original QCD axion for the strong CP problem [39], and we will use the term “axion” for any light (pseudo)scalar dark matter that has similar interactions to the QCD axion. Another interesting candidate is the dark photon dark matter that can be produced from a variety of mechanisms [40, 41, 42, 43, 44, 45, 46]. While we focus on ultralight axions in this work, similar physical phenomena often apply to the dark photon and axion dark matter in general.

For wave dark matter on the higher mass end, its de Broglie wavelength is much shorter than astrophysical scales, and laboratory experiments are necessary. Once the wave dark matter is ultralight (often called fuzzy dark matter [47]), cosmological measurements on the linear power spectrum or stellar kinematics of ultrafaint dwarfs can constrain lower mass ranges [47, 48, 49, 50, 51, 52, 53]. The wave nature of fuzzy dark matter can lead to rich phenomenology such as the formation of soliton cores at halo centers and interference effects [54, 55, 56, 57, 58, 59]. Therefore, fuzzy dark matter can also be probed by compact objects through its wave dynamics [60, 61, 62, 63, 64, 65, 66], and the nature of its couplings with visible matter can be constrained by various observables [67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 99, 101, 63, 64, 65, 66, 102, 103].

Though usually not thermally produced, fuzzy dark matter can still have a significant relativistic component, for instance post-inflation axions produced from relaxation of string networks. Ref. [104] pointed out that in such a relativistic regime, wave dark matter exhibits free streaming behavior much like the collisionless damping of warm or hot particle dark matter [105]. Such dark matter is thus warm and fuzzy simultaneously.

Ref. [104] highlights the apparent differences between free streaming behavior associated with the wavenumber distribution for wave dark matter and that with the particle momenta distribution for particle dark matter, and thus their respective effects on cosmological perturbations. These apparent differences are important to understand when applying free streaming considerations to bounds on the dark matter mass and the evolution of isocurvature fluctuations from post-inflation causal production as compared to cold dark matter isocurvature perturbations.

In this work, we further explore the relationship between the free streaming of wave and particle dark matter and resolve their apparent differences. We begin in §II by relating the particle and wave pictures of free streaming and the impact of wavenumber versus particle momentum distribution on the transfer function of density perturbations. We show that axion wave dark matter produced after inflation is warm in this sense and only moderately enhances the Jeans or free streaming damping already present for initially cold axions. In §III, we study with simulations the effect of free streaming on the causally produced isocurvature fluctuations of an even hotter, i.e. more relativistic, wave dark matter than axions, and resolve the paradox that the effective number density fluctuations do not damp even though the waves that compose them do – for particles, the initial number density fluctuations are averaged out over the free streaming volume; for waves, free-streaming damping causes the momentum or wavenumber distribution to become incoherent, effectively transferring power from spatial inhomogeneities to anisotropies in the momentum distribution. In §IV, we show that the impact of the incoherence of these fluctuations prevents their appearance in time- or spatially-averaged quantities, and should be thought of as the wave analogue of multiple streams in the phase space density. We discuss the implication of these results in §V and provide appendices on the computation of free streaming (App. A) for general wavenumber or momentum distributions and their impact on mass bounds (App. B). Throughout we employ units where $\hbar=c=1$.

II Free Streaming Duality

The free streaming of wave and particle dark matter shares the same underlying principles and can impact cosmological structure formation on scales smaller than the maximal free streaming scale. Given a particle mass $m$ and comoving momentum $q$, a non-interacting particle will stream with a velocity

$$v=\frac{q}{\sqrt{q^{2}+a^{2}m^{2}}},$$ (1)

where $a$ is the scale factor. Similarly, a free wavepacket of a field $\phi$ that obeys a relativistic wave equation such as the Klein-Gordon equation, with a dispersion relation $\omega=\sqrt{q^{2}+a^{2}m^{2}}$, will propagate at the group velocity

$$v=\frac{\partial\omega}{\partial q}=\frac{q}{\sqrt{q^{2}+a^{2}m^{2}}}$$ (2)

around a comoving wavenumber $q$. We use the terms momentum and wavenumber interchangeably throughout.

In both cases the free streaming length becomes

$$\lambda_{\rm fs}(q;a)=\int d\tau v(q,\tau)=\int\frac{d\ln a}{aH(a)}\frac{q}{\sqrt{q^{2}+a^{2}m^{2}}},$$ (3)

where $\tau=\int dt/a$ is the conformal time. Where no confusion should arise, we suppress the evaluation scale factor (“$a$”) in the argument of functions, e.g. $\lambda_{\rm fs}(q)\equiv\lambda_{\rm fs}(q;a)$.

For ultrarelativistic momenta $q\gg am$, $v\approx 1$ and $\lambda_{\rm fs}(q)\approx\tau$. For nonrelativistic momenta $q\ll am$, $v\approx q/am$, and the free streaming length grows logarithmically from its value of $\tau|_{a=q/m}$ during radiation domination and ceases to grow during matter domination. It is therefore convenient to scale $\lambda_{\rm fs}(q)$ to the comoving Hubble length at equality $a=a_{\rm eq}$

$$\lambda_{\rm fs}(q;a)\approx\frac{\sqrt{2}}{a_{\rm eq}H_{\rm eq}}F\left(\frac{q}{a_{\rm eq}m};\frac{a}{a_{\rm eq}}\right)\equiv\frac{\sqrt{2}}{a_{\rm eq}H_{\rm eq}}F(\hat{q};y),$$ (4)

where

$$\frac{\sqrt{2}}{a_{\rm eq}H_{\rm eq}}=\frac{1}{H_{0}}\sqrt{\frac{a_{\rm eq}}{\Omega_{m}}},$$ (5)

and carry the scaling behaviors in the various regimes with the dimensionless function $F$. The exact analytic form of this function and its scaling behaviors are given in Appendix A.

Of particular interest for viable dark matter models that mimic CDM on large scales are candidates that become nonrelativistic well before equality. The maximal scale for the impact of free streaming is the value that $\lambda_{\rm fs}$ achieves well after equality. Combining these two limits, we find the asymptotic approximation

$$F(\hat{q};y) \approx \hat{q}\ln(8/\hat{q})$$

$$= \frac{q}{a_{\rm eq}m}\ln\left(\frac{8a_{\rm eq}m}{q}\right),\quad q\ll a_{\rm eq}m,a\gg a_{\rm eq}.$$

The log term represents the logarithmic growth from the epoch that the momenta become nonrelativistic $a_{\rm nr}\sim q/m$ through $a_{\mathrm{eq}}$, and the $\hat{q}=q/a_{\mathrm{eq}}m$ prefactor likewise scales the comoving horizon at equality to $a_{\rm nr}$ given $\sqrt{2}aH(a)/a_{\rm eq}H_{\rm eq}=a_{\rm eq}/a$ in radiation domination.

The distinction between various types of dark matter therefore mainly comes from their momentum or wavenumber distributions. For thermally produced dark matter, this comes from the distribution at the relevant temperature for production and kinetic decoupling. For non-interacting scalar wave dark matter $\phi$, the number density scales as $m\phi^{2}$ in the nonrelativistic regime, and thus the momentum spectrum is provided by the power spectrum of $\phi$ itself. Below we will use the terms power spectrum of field fluctuations and momentum distribution of the number density interchangeably where no confusion should arise.

For the axion, if the Peccei-Quinn symmetry breaking occurs after inflation, the axion field is uncorrelated across different horizon patches, resulting in white-noise field fluctuations above the horizon. At the critical time when axions acquire masses their potential becomes

$$V(\phi)=m^{2}f_{\rm ax}^{2}[1-\cos(\phi/f_{\rm ax})]$$ (7)

and once $H(a_{*})=m$, the potential energy stored in the random initial field $\phi/f_{\rm ax}\in[-\pi,\pi)$ will convert to locally coherent oscillations of the axion field, which produces mostly cold axions whose spatial number density varies from horizon patch to patch. This mechanism is known as vacuum misalignment production.

Furthermore, the post-inflationary axion also predicts the existence of topological defects such as axion strings due to the Kibble mechanism [106]. The string network evolves in such a pattern that the number of strings per horizon is nearly constant, and the energy stored in string cores is lost through the radiation of relativistic axion waves [107, 108]. This emission may contribute significantly to the axion relic density [109, 110, 111] and extend the axion momentum distribution to $q>a_{*}m$, providing a “warm” component.

In the post-inflationary case, the power spectrum of field fluctuations then gives the momentum spectrum of the average number density of axions after the relevant momenta become non-relativistic

$$\langle\phi^{2}\rangle=\int\frac{d^{3}q}{(2\pi)^{3}}P_{\phi}(q)=\int d\ln q\frac{q^{3}}{2\pi^{2}}P_{\phi}(q)$$ (8)

with

$$\langle\phi({\bf q})\phi({\bf q}^{\prime})\rangle=(2\pi)^{3}\delta({\bf q}+{\bf q}^{\prime})P_{\phi}(q).$$ (9)

This spectrum is white $P_{\phi}=$ const. for $q\ll q_{*}\equiv a_{*}m$. For $q\gg q_{*}$, the axions produced by misalignment and by decay of the scaling string network at higher redshifts $a<a_{*}$ dilute their number density as $n\propto a^{-3}$ (e.g. [112]) leading to the scaling expectation $dn/d\ln q\propto q^{-1}$ [107]. Following [104], we combine these behavior for $a\gg a_{*}$ as

$$q^{3}P_{\phi}(q)\propto\left(\frac{q}{q_{*}}\right)^{3}\theta(q_{*}-q)+\left(\frac{q_{*}}{q}\right)^{\alpha}\theta(q-q_{*}).$$ (10)

Simulations differ on whether the $q>q_{*}$ power law $q^{-\alpha}$ is strictly the scaling value of $\alpha=1$ and therefore how much of the string network energy is radiated at a given momentum, which can make a large change in the overall relic number abundance [109, 110, 111]. Notice however that as long as $\alpha>0$, the number density spectrum is still dominated by momenta around $q_{*}$ as is the energy density $\rho\approx mn$ after all momenta are non-relativistic (cf. Eq. 31).

On the other hand around $q_{*}$, we have simply joined the two asymptotic behaviors as a broken power law spectrum. In simulations of string dynamics, the spectrum around $q_{*}$ is smoother and can have transient plateau-like features before the asymptotic $q\gg q_{*}$ break [108, 109, 113]. In the main paper we will simply assume this broken power law form with $\alpha=1$ and in App. B we explore variations and their consequences (see also [104]v2, their App. B).

Since this number density spectrum of axions is peaked around $q_{*}$, we can expect that the net effect of averaging the free streaming of the momenta components in the spectra is dominated by these $\sim q_{*}$ modes, which are only quasirelativistic or “warm” at birth. Correspondingly, we would expect the impact of free streaming on density perturbations to occur at (see App. B and Fig. 10)

$$k_{\rm fs} \equiv \lambda_{\rm fs}^{-1}(q_{*}=a_{*}m)=\frac{\sqrt{a_{\rm eq}^{2}mH_{\rm eq}}}{2^{1/4}\ln(8a_{\rm eq}/a_{*})}.$$ (11)

Throughout, when we compute numerical values for such quantities we take a cosmology with matter density $\Omega_{m}h^{2}=0.142$ and 3 massless neutrinos.

This scale should be compared to the similar suppression of the transfer function at the Jeans scale in the case of Pecci-Quinn symmetry breaking before the end of inflation. Here the initial misalignment is coherent across the whole horizon volume today by the end of inflation and there is only an initially cold component to the axions. Curvature fluctuations of wavenumber $k$ then imprint field fluctuations $\phi(q)$ at wavenumber $q=k$ and the density fluctuations are carried by $\phi^{2}(k)\approx 2\phi(k)\langle\phi\rangle$. Throughout, “$k$” identifies the wavenumber of quadratic quantities such as $\phi^{2}$, number, and energy density, whereas “$q$” when different from $k$ distinguishes the wavenumber of field fluctuations that compose them.

Figure 1: Relative transfer function (15) due to the effect of free-streaming on the “warm” axions of the post-inflationary mechanism compared with that of the cold axions of the pre-inflationary case. Here $m=2.0\times 10^{-20}$ eV and the warm spectrum peaks at $q_{*}=a_{*}m$, the horizon wavenumber at the start of axion oscillations, and the evaluation epoch $a=0.2$ is chosen to reflect that of the Lyman-$\alpha$ forest.

The relevant free streaming scale for the pre-inflationary case is the comoving wavenumber whose associated free streaming length overtakes its wavelength

$$\lambda_{\rm fs}(k_{J})\approx\lambda_{J}\approx\sqrt{6}/k_{J},$$ (12)

and by employing Eq. (II), we have

$$k_{J}\approx 3^{1/4}\sqrt{\frac{a_{\rm eq}^{2}mH_{\rm eq}}{\ln(8a_{\rm eq}/a_{*})}}.$$ (13)

Note that the $\sqrt{6}$ in Eq. (12) is added so that Eq. (13) matches the definition in the literature [47] (their Eq. 9), modulo the log factor which we have here approximated at $k_{J}\sim a_{*}m=2^{-1/4}a_{\rm eq}\sqrt{mH_{\rm eq}}$ but is usually incorporated more precisely as a mass-dependent fitting factor to numerical calculations of the transfer function ([114] their Eq. 44).

We therefore expect that the free streaming of the quasirelativistic or “warm” axions in the post-inflationary scenario to scale in the same way as the Jeans scale of cold axions in the pre-inflationary scenario and differ only by the ratio

$$\frac{k_{J}}{k_{\rm fs}}={6^{1/4}}\sqrt{\ln(8a_{\rm eq}/a_{*})}.$$ (14)

We can improve on this estimate by averaging over the momentum spectrum $q^{3}P_{\phi}$ instead of evaluating at the peak to define an effective transfer function for density perturbations due to the free streaming effect following Ref. [104]

$$T_{\rm rel}(k)\equiv\frac{T_{\rm ax}}{T_{\rm CDM}}(k)\approx\frac{\int d\ln q\frac{q^{3}}{2\pi^{2}}P_{\phi}(q)\frac{\sin[k\lambda_{\rm fs}(q)]}{k\lambda_{\rm fs}(q)}}{\int d\ln q\frac{q^{3}}{2\pi^{2}}P_{\phi}(q)},$$ (15)

who derive this form from the WKB approximation for the amplitude of the free-streaming waves and linearizing adiabatic perturbations as a means of estimating where free streaming has an $O(1)$ effect. Here $T_{\rm CDM}$ is the cold dark matter (CDM) density transfer function (with no axions) and $T_{\rm ax}$ is the axion transfer function (with no CDM), with the same cosmological parameters otherwise. This ratio serves to isolate the free streaming effect by removing other cosmological effects from matter radiation equality and baryon acoustic oscillations. The amplitude reduction term can also be motivated from the treatment of wave propagation in §III where free streaming modifies the initial wave amplitude for each momentum $q$ according to the solution to the wave equation (see Eq. (20)). We will discuss the impact of the slowly varying phases of the momentum components in §IV.

Figure 2: Mass scaling of the relative transfer function of warm axions relative to the cold axions in Fig. 1. As predicted from Eqn. (14), a warm axion mass of $5.6\times 10^{-19}$ eV produces a comparable scale of suppression to cold axions of $2.0\times 10^{-20}$ eV.

In Fig. 1, we compare the relative transfer function (15) with the usual Jeans suppression for cold axions from a numerical calculation using a modified version of CAMB[115] (see also [114] their Eq. 44, which closely match these results). We illustrate this with a mass of $m=2.0\times 10^{-20}$eV which is motivated by the bound on cold axions from the Lyman$-\alpha$ forest [116]. Correspondingly we take a redshift of $z=4$ or $a=0.2$ as the evaluation epoch. As expected, the free streaming transfer function for the “warm” axions gives a stronger suppression than the cold, pre-inflationary case but only by a log factor. In fact, Eq. (14) predicts $k_{J}/k_{\rm fs}\approx 5.3$, which captures most of the difference. This ratio can be then be used to approximately scale up any given Lyman-$\alpha$ bound on cold axions since $k_{J}\propto m^{1/2}$, here nominally $m\gtrsim 5.6\times 10^{-19}$eV. Fig. 2 demonstrates that warm axions of this mass give a transfer function comparable to the “cold” axions of $m=2.0\times 10^{-20}$ eV.

For heavier masses, where free streaming is negligible on observationally relevant scales, the random number of cold axions in each horizon patch at $a_{*}$ leads to so-called isocurvature fluctuations on large scales, which is also constrained by the Lyman-$\alpha$ forest observations. At large scales, the isocurvature perturbation is well described by the white-noise power spectrum and is not sensitive to the behavior of the power spectrum at $k\sim q_{*}$ which determines the free-streaming effect. Therefore, the ratio of the amplitude of the isocurvature fluctuations to that of the adiabatic is $f_{\rm iso}^{2}\propto 1/q_{*}^{3}$. A lower $q_{*}$, corresponding to a lighter axion in this range, gives a larger isocurvature to adiabatic ratio on large scales, imposing a lower bound on the axion mass in the post-inflationary scenario, $m>3\times 10^{-18}$ eV from Lyman-$\alpha$ forest constraints [117], which is stronger than the nominal free streaming bound above. For this mass, the free-streaming length $\lambda_{\rm fs}(q_{*})=0.014$ Mpc, which implies that the free streaming effect can be neglected when applying this isocurvature bound.

III Causally Coherent Patches

We have seen in the previous section that for axion dark matter produced around the beginning of the oscillation epoch $H(a_{*})=m$ with only a mildly relativistic or “warm” component, the free-streaming scale for the post-inflationary production mechanism remains close to that of cold axions from the pre-inflationary production mechanism. Thus, the main new large scale effect in the post-inflationary case is that the causally random patches on the horizon scale at $a_{*}$ produce a white noise spectrum of axion number or mass density fluctuations.

For a hypothetical wave dark matter candidate produced causally after inflation with more relativistic momenta, the free streaming scale can be larger; simultaneously, the white noise fluctuations can be smaller in amplitude on a fixed observational scale given the smaller horizon scale at production, and therefore weaken the isocurvature bounds on the mass until the free streaming limit comes to dominate [104].

More concretely as shown in App. B, if the number density spectrum peaks at $q_{\rm peak}=Rq_{*}=Ra_{*}m$, then the effective free streaming length increases as $\lambda_{\rm fs}(q_{\rm peak})\propto R$ up to a log correction using Eq. (II). For moderate increases where $R\sim{\cal O}(1)$, this can make free streaming more important and become competitive or stronger than the isocurvature bound. The combination of the two is therefore more robust for changes in the axion spectrum for masses in the $m\sim 10^{-18}-10^{-19}$eV range [104].

For $R\gg 1$, most of the scalar wave dark matter is ultrarelativistic at $a_{*}$. In this case, the dark matter is “hot” at birth. Note that this does not occur with axions, since at any given epoch after Peccei-Quinn symmetry breaking but before $H\sim m$, the Kibble mechanism establishes a coherent field with a random value of $\phi/f_{\rm ax}\in[-\pi,\pi)$ across the Hubble patch $1/aH$ with a self similar network of strings and their decay products.

In the more general case, it is also important to understand the evolution of the isocurvature density fluctuations from the causally random initial field fluctuations, since free streaming suppression and isocurvature enhancement work in opposite directions. In Ref. [104], it was shown that a certain characterization of these field fluctuations, namely fractional fluctuations in $\phi^{2}$, remains constant and white in power despite the free streaming of the waves that compose them. In this section, we seek to clarify the nature of these fluctuations from the standpoint of the free streaming of the causally coherent initial horizon scale patches.

To understand the impact of free-streaming on causally coherent patches, let us begin with an initial field profile that is a spherical tophat of radius $\tau_{i}$ and set the field normalization to unity at $r<R$. The Fourier transform of the tophat gives the momentum distribution of the patch

$$\phi_{p}(q;0)=\frac{3V_{p}}{(q\tau_{i})^{3}}\left[{\sin(q\tau_{i})-q\tau_{i}\cos(q\tau_{i})}\right],$$ (16)

where $V_{p}=4\pi\tau_{i}^{3}/3$ is the volume of the tophat patch. The spectrum is a constant $\phi_{p}(q;0)=V_{p}$ for $q\tau_{i}\ll 1$ and a random distribution of these causal patches would produce white noise on large scales as desired.

Figure 3: Free streaming of 4 coherent tophat field patches (see text) for the field $\phi$ (top panel) and the proxy for number density fluctuations $\phi^{2}/\langle\phi^{2}\rangle$ (bottom panel) for a series of snapshots in conformal time $\tau$. As the waves from different patches free stream, they individually damp in amplitude and spread in scale. When multiple wave fronts intersect and superimpose, fluctuations appear that are transient and reflect multiple momentum streams from different patches.

Figure 4: Power spectra of $\phi$ and $\phi^{2}/\langle\phi^{2}\rangle$ for the 4 patch case of Fig. 3. Free streaming damping in both cases scales as $k\propto\lambda_{\rm fs}^{-1}=\tau^{-1}$ and here the phase-coherent autocorrelation of momentum modes from individual patches dominates.

The Klein-Gordon equation in $q$-space for a free field

$$\ddot{\phi}+2aH\dot{\phi}+(q^{2}+a^{2}m^{2})\phi=0,$$ (17)

evolves the initial modes, where overdots denote conformal time derivatives. More specifically, each mode evolves via the growth function

$$\phi(q;\tau)=\phi(q;0)D(q;\tau),$$ (18)

where $D$ solves the Klein-Gordon equation (17) with an initially frozen field due to the Hubble drag, $\dot{\phi}(q;0)=0$. During radiation domination, $\tau=1/aH$ and $a^{2}m^{2}=\tau^{2}/\tau_{*}^{4}$, and the growth function of the field is given by

$$D(q;\tau)=e^{-i(\tau/\tau_{*})^{2}/2}{}_{1}F_{1}\left[\frac{3}{4}+i\frac{(q\tau_{*})^{2}}{4},\frac{3}{2},i\left(\frac{\tau}{\tau_{*}}\right)^{2}\right],$$ (19)

where ${}_{1}F_{1}$ a hypergeometric function. For $\tau\ll\tau_{*}$, $D(q;\tau)$ takes the simple form

$$\lim_{\tau/\tau_{*}\ll 1}D(q;\tau)=\frac{\sin(q\tau)}{q\tau},$$ (20)

(cf. Eq. 15). Here the field is frozen above the horizon $q\tau\ll 1$ and oscillates with amplitude $D\propto 1/\tau\propto 1/a$ below the horizon. This amplitude decay reflects the redshifting of relativistic waves inside the horizon. Furthermore, the number density associated with $q$ goes as $n\propto(\omega/a)\phi^{2}\propto a^{-3}$, and particle number in each mode is conserved in comoving coordinates.

With the initial tophat wavepacket, Fourier transforming the product of Eq. (16) and (20) gives in real space

$$\phi_{p}(r;\tau)=\begin{cases}1,&r<\tau_{i}-\tau,\tau<\tau_{i}\\ 0,&r<\tau-\tau_{i},\tau\geq\tau_{i}\\ \frac{(\tau_{i}+\tau-r)(\tau_{i}-\tau+r)}{4r\tau},&|\tau_{i}-\tau|\leq r<\tau_{i}+\tau\\ 0,&r>\tau_{i}+\tau\end{cases}$$ (21)

which represents a spherically symmetric shell expanding radially where $r=|{\bf x}|$ with a front at $\tau_{i}+\tau$, reflecting causal propagation at the speed of light. As expected, the field amplitude damps as it spreads outwards and transfers its coherent fluctuations to the larger physical scales associated with the free streaming scale $\lambda_{\rm fs}=\tau$.

For $\phi_{p}^{2}$, its Fourier components are composed by a convolution of the field momenta

$$\phi_{p}^{2}({\bf k};\tau) = \int\frac{d^{3}q}{(2\pi)^{3}}\phi_{p}({{\bf q}};\tau)\phi_{p}({{\bf k}-{\bf q}};\tau)$$ (22)

$$= \int\frac{d^{3}q}{(2\pi)^{3}}\frac{\sin(q\tau)}{q\tau}\frac{\sin(|{\bf k}-{\bf q}|\tau)}{|{\bf k}-{\bf q}|\tau}$$

$$\times\phi_{p}({{\bf q}};0)\phi_{p}({{\bf k}-{\bf q}};0).$$

Since the initial profile is coherent at $r<\tau_{i}$, different field momenta $\bf q$ are correlated in their phase and coherently superimpose in this quadratic combination to produce the spatially coherent fluctuation $\phi_{p}^{2}({\bf x};\tau)$. As with $\phi_{p}$, the power spectrum of $\phi_{p}^{2}$ is strongly damped by free streaming and represents the dilution or averaging out of the coherent fluctuation in a given patch. Interpreted in the particle picture, the initial axion number fluctuation in a given coherent patch is averaged out over the free streaming scale.

On the other hand, the total $\phi$ is a sum over all horizon patches, each with a random amplitude, and $\phi^{2}$ receives contributions not just from the coherent propagation of modes of a single patch but also all of the phase-incoherent cross terms between patches. In this case, there is a sum over $N$ patches:

$$\phi({\bf q};\tau) = \sum_{\alpha=1}^{N}\phi_{\alpha}({\bf q};\tau),$$

$$\phi_{\alpha}({\bf q};\tau) \equiv A_{\alpha}\phi_{p}({\bf q};\tau)e^{i{\bf q}\cdot{\bf d}_{\alpha}},$$ (23)

where $A_{\alpha}$ is the field value at the center of patch $\alpha$ at spatial coordinate ${\bf d}_{\alpha}$. Correspondingly,

$$\phi^{2}({\bf k};\tau) = \sum_{\alpha\beta}\phi_{\alpha\beta}^{2}({\bf k};\tau),$$

$$\phi_{\alpha\beta}^{2}({\bf k};\tau) \equiv \int\frac{d^{3}q}{(2\pi)^{3}}\frac{\sin(q\tau)}{q\tau}\frac{\sin(|{\bf k}-{\bf q}|\tau)}{|{\bf k}-{\bf q}|}$$ (24)

$$\times\phi_{\alpha}({{\bf q}};\tau)\phi_{\beta}({{\bf k}-{\bf q}};\tau),$$

now has autocorrelation terms between modes from the same patch $\alpha=\beta$ and cross terms between patches $\alpha\neq\beta$. Both the evolution and the physical interpretation of the auto and cross terms differ. In particular, for the autocorrelation terms, the power spectrum of $\phi^{2}_{\alpha\alpha}=\phi^{2}_{p}$ contains phase correlations between differing momenta and damps with free streaming, whereas the cross terms carry incoherent phase shifts due to the displacements ${\bf d}_{\alpha}$, though these vanish when pairing the same momenta in the $\phi^{2}$ power spectrum, i.e. $\phi({\bf q})\phi(-{\bf q})$.

Figure 5: Free streaming evolution as in Fig. 3 but with pixel scale horizon patches and field values drawn from a uniform distribution. While the field $\phi$ free streams and becomes smoother and lower in amplitude with time, $\phi^{2}/\langle\phi^{2}\rangle$ becomes statistically time invariant with transient fluctuations at a given position.

To visualize the difference, consider explicitly a simple superposition of such patches. In Fig. 3 we set up 4 patches in a periodic box of length $64\tau_{i}$ on each side, with displacements $d_{\alpha}$ of $\pm 4\tau_{i}$ from the center of the box in the $x$ and $y$ directions, and positive and negative amplitudes respectively such that the total field has zero mean. The box is represented by $512^{3}$ discrete pixels, with $8$ pixels per unit $\tau_{i}$. We evolve this configuration until $\tau=8\tau_{i}$, again by employing Eq. (20) in Fourier space. We have tested the simulation procedure by verifying that the time evolution of one such patch is consistent with the analytic radial solution (21).

In Fig. 3 (top panels) we show the field profile itself in a $z=0$ slice of the box. Notice that the free streaming of the individual patches brings the wavefronts to intersect at the center of the box at around $\tau=4\tau_{i}$. Therefore the momentum distribution of the field fluctuations at the center is anisotropic, specifically quadrupolar. More generally, at any given time after the free-streaming intersection of fronts, the total field represents a transient superposition of waves composed of multiple field momentum streams at any given physical position. This is the same behavior as the particle free streaming of CMB photons after recombination or cosmic neutrinos after their decoupling: the initial particle number inhomogeneity becomes a phase space anisotropy after free streaming. In those cases, the total power in fluctuations of a given $k$-mode is conserved (e.g. [118], their Eq. 10), but its nature and effect on gravitational structure formation differ qualitatively. In the free-streaming damping context, this is known as the directional damping of collisionless components [105].

These free streaming considerations apply to $\phi^{2}$ as well. In Fig. 3 (bottom panels), we show $\phi^{2}$ normalized to its average in the box $\langle\phi^{2}\rangle$. This normalization removes the redshifting effect of the subhorizon modes and we can see the remaining strong effect of free streaming damping of the amplitude of $\phi^{2}/\langle\phi^{2}\rangle$ of each patch, from the change in scale of the panels with $\tau$. Moreover, even though $\phi^{2}({\bf x})$ itself is not a directional quantity, its local value reflects the directionally dependent propagation of the fronts of each of the 4 patches.

Figure 6: Power spectra of the pixel scale horizon patches of Fig. 5 as in Fig. 4. While the field power (top panel) continues to evolve and oscillate with free streaming, the power in $\phi^{2}/\langle\phi^{2}\rangle$ (bottom panel) approaches a constant power law behavior below the free streaming scale.

The corresponding power spectra for $\phi$ and $\phi^{2}/\langle\phi^{2}\rangle$ are shown in Fig. 4. Note that $P_{\phi^{2}/\langle\phi^{2}\rangle}=P_{\phi^{2}}/\langle\phi^{2}\rangle^{2}$. In this 4-patch case, the power spectra themselves are still dominated by the autocorrelation terms of each patch, and the total $\phi^{2}$ still strongly damps with free streaming. We can see that the turnover into the free streaming oscillation behavior scales with the free streaming length $\lambda_{\rm fs}=\tau$, $k\propto 1/\tau$, as expected.

On the other hand, as the number of patches $N$ grows, the number of cross terms grows as $N^{2}$. In fact, since the number of patches that fill a volume $V$ will be $N=V/V_{p}$, it is the cross terms that become the dominant source of fluctuations in $\phi^{2}$. At $\tau=0$ the auto and cross correlation contributions to the power spectrum of $\phi^{2}$ are comparable, but the autocorrelation terms free stream away at later times. These remaining contributions represent the incoherent superposition of fluctuations of different momenta ${\bf q}$.

From Eq. (III), we can see that cross terms have no time averaged effect on $\phi^{2}$ since $q$ and $|{\bf k}-{\bf q}|$ modes oscillate incoherently, but provide a source of instantaneous power

$$P_{\phi^{2}}(k;\tau) \approx 2\int\frac{d^{3}q}{(2\pi)^{3}}\left[\frac{\sin(q\tau)}{q\tau}\frac{\sin(|{\bf k}-{\bf q}|\tau)}{|{\bf k}-{\bf q}|}\right]^{2}$$ (25)

$$\times P_{\phi}({{\bf q}};0)P_{\phi}({{\bf k}-{\bf q}};0),$$

since the square of the oscillating growth function is positive definite. Here we have dropped the autocorrelation terms that contain the phase coherence between the momentum modes, i.e. the connected pieces of the trispectrum of $\phi$.

This behavior can be seen directly in Fig. 5 where we take $\tau_{i}$ to be the grid spacing such that each pixel represents a horizon patch. The box size is $256\tau_{i}$ in this simulation and the field value at each pixel is a uniform random deviate with $\phi\in[-\pi,\pi)$. Beyond the axion context, we have also separately checked that a Gaussian random deviate produces qualitatively the same free streaming effects we describe below but with a larger number of rare high density peaks at the pixel scale. Again the top panels in Fig. 5 show the time evolution of $\phi$ and the bottom panels that of $\phi^{2}/\langle\phi^{2}\rangle$. Instead of the coherent propagation of discrete horizon scale patches, we are now dominated at late times by the cross terms between the $256^{3}$ initial pixel scale patches.

Notice also that the statistical properties of the $\phi^{2}({\bf x})/\langle\phi^{2}\rangle$ field become nearly time-invariant. This is in contrast to the 4 patch case even though both cases show the field $\phi({\bf x})$ continuing to evolve, especially in amplitude, as their respective patches free stream.

We quantify this in Fig. 6 for the respective power spectra. The initially white $q^{3}P_{\phi}(q)\propto q^{3}$ turns over to oscillate with a $\propto q^{1}$ scaling and decreasing amplitude which reflects the redshifting behavior of the growth function $D$ in Eq. (20) as with the 4 patch case. On the other hand, the power spectrum of $\phi^{2}$ gains a non-oscillating and nearly constant in time $k^{3}P_{\phi^{2}/\langle\phi^{2}\rangle}\propto k^{2}$ spectrum on scales below the free streaming scale. This reflects the scaling of Eq. (25) after normalization by $\langle\phi^{2}\rangle$ which removes the redshifting effect.

In particular, $\phi^{2}$ retains fluctuations across a range of scales below the free streaming scale but above the initial horizon scale $\tau_{i}$. To better see this, we low pass filter $\phi^{2}$ to retain only $k\tau_{i}<0.3$ and show two late time snapshots $\tau=16\tau_{i}$ and $24\tau_{i}$ in Fig. 7. Notice that even though the power in fluctuations on these scales is nearly constant, the $\phi^{2}$ field evolves on a time scale much shorter than the Hubble time reflecting the chance superposition of the underlying free-streaming modes. Again, in the particle picture, this reflects a phase space anisotropy in the distribution rather than a physical space inhomogeneity. Also, it is important to note that a low pass filter in $k$ for $\phi^{2}$ is not in general the same as a low pass filter in $q$ in $\phi$ since high momenta modes can still contribute to low wavenumbers ${\bf k}={\bf q}_{1}+{\bf q}_{2}$ if ${\bf q}_{1}\approx-{\bf q}_{2}$, which as we shall see below is the physically relevant case after all populated modes have become nonrelativistic.

Figure 7: Long wavelength fluctuations in $\phi^{2}$ at late times as in Fig. 5 but low pass filtered to $k\tau_{i}<0.3$. Although statistically the same, the two snapshots at $\tau=16\tau_{i}$ and $24\tau_{i}$ show strong time evolution that reflects the transient nature of incoherent superposition.

In fact, during the simulated epochs where the $q$-modes are still ultrarelativistic, this free streaming behavior fully parallels that of relativistic particles and the remaining phase space fluctuations would not contribute to the gravitational formation of structure.

The distinction between wave dark matter and the relativistic particle case is that $\phi$ is a massive field and even for initially relativistic $q$ modes, the redshifting due to cosmic expansion will eventually make the modes nonrelativistic, much like the massive neutrino component of dark matter in $\Lambda$CDM. For viable dark matter models, this occurs well before matter radiation equality and we must consider the impact of their further evolution until non-relativistic. Since these modes are ultrarelativistic at $\tau_{*}$ where $H=m$ by construction, this means that we must consider their evolution at $\tau>\tau_{*}$ and account for the change in the spectrum as progressively larger $q$ modes become nonrelativistic.

For $\tau>\tau_{*}$, even the superhorizon modes that are non-relativistic at $\tau_{*}$ evolve with $D(q;\tau)$ reflecting the coherent oscillations of the axion field due to the mass term. For modes that are above the horizon at $\tau_{*}$, $q\tau_{*}\ll 1$, the growth function becomes

$$\lim_{q\tau_{*}\ll 1}D(q;\tau)=\sqrt{\frac{2\tau_{*}}{\tau}}\Gamma\left(\frac{5}{4}\right)J_{1/4}(\tau^{2}/2\tau_{*}^{2}).$$ (26)

Note that $(\tau/\tau_{*})^{2}/2=mt$ where $t$ is coordinate time so that the Bessel function carries the mass scale oscillations for $mt=m/2H\gg 1$. The amplitude of these oscillations scale as $D\propto\tau^{-3/2}\propto a^{-3/2}$. This behavior again reflects the redshifting of non-relativistic matter $n\propto(\omega/a)\phi^{2}\propto m\phi^{2}\propto a^{-3}$. Notice that for modes that remain relativistic until after $\tau>\tau_{*}$, even though their number density $n\propto a^{-3}$, the extra redshifting of the frequency means that the field fluctuations decay less quickly as $a^{-1}$ instead of $a^{-3/2}$.

In general, the $k$-modes of $\phi^{2}$ are then constructed out of the $q$-modes of $\phi$ as

$$\phi^{2}({\bf k};\tau) = \int\frac{d^{3}q}{(2\pi)^{3}}D(q;\tau)D(||{\bf k}-{\bf q}|;\tau)$$ (27)

$$\times\phi({{\bf q}};0)\phi({{\bf k}-{\bf q}};0),$$

and the different behavior of $D$ as a function of momentum $q$ changes the weights of the field fluctuations that factor into a given $k$. On the other hand, since these subhorizon modes simply redshift, the number density simply dilutes with the expansion for all modes after $a\gg a_{*}$. We can therefore infer the later behavior of all relevant momenta modes directly from the relativistic simulations using number conservation instead of explicitly extending them to late times using computationally expensive evaluations of Eq. (19) for the general growth function $D$.

More specifically, during an epoch where the relevant $q$ modes are still relativistic and $\omega\phi^{2}/a\approx q\phi^{2}/a\propto q\phi^{2}/\tau$, we can define an effective number density

$$n_{\rm eff}({\bf k};\tau) \propto \int\frac{d^{3}q}{(2\pi)^{3}}\frac{\sin(q\tau)}{\sqrt{q\tau}}\frac{\sin(|{\bf k}-{\bf q}|\tau)}{\sqrt{|{\bf k}-{\bf q}|\tau}}\frac{1}{\tau^{2}}$$ (28)

$$\times\phi_{p}({{\bf q}};0)\phi_{p}({{\bf k}-{\bf q}};0),$$

which then reflects the spectrum of $\phi^{2}({\bf k})$ and their weights in $\phi({\bf q})$ after the waves have become nonrelativistic. Notice that given an initial white noise spectrum for $q\lesssim 1/\tau_{i}$, the field modes are no longer white after free streaming, but the integral in Eq. (28) is still dominated by $q\sim 1/\tau_{i}$ and these modes produce a white noise spectrum in $n_{\rm eff}/\langle n_{\rm eff}\rangle$ for $k\lesssim 1/\tau_{i}$ that is constant in time. Therefore, the white noise power at $k\ll 1/\tau_{i}$ still comes mainly from momentum modes $q\sim 1/\tau_{i}$ (see discussion of Fig. 7), which as we will see in the next section is analogous to the beat frequency from the superposition of closely spaced high frequency modes.

We illustrate this behavior and construction using the random pixel simulations of Fig. 6 in Fig. 8. Notice that the power spectrum now remains constant and white for the effective number density and now reflects $\phi^{2}$ after all of the relevant momenta are nonrelativistic.

On the other hand, this constancy of the power spectrum of $n_{\rm eff}/\langle n_{\rm eff}\rangle$ should not be equated with a particle number density in real space since its spectrum is still composed of field modes with different $q$, both in magnitude around $1/\tau_{i}$ and in direction, due to the incoherent superposition of contribution from individual horizon scale patches, i.e. the wave analogue of a phase space distribution.

Figure 8: Effective number density power spectrum (see Eq. 28) for the simulation case of Fig. 5. Since number density is conserved during the transition of each momentum mode from relativistic to nonrelativistic, the effective number density weights the momentum modes according to the final dark matter density, and its fractional power spectrum remains constant and white at late times despite the free streaming of the underlying modes themselves.

IV Incoherent Superposition

We have seen in the previous section that the effective number density fluctuations in a free scalar field $\phi$, with an initial white noise spectrum from causal production, do not damp by free streaming. On the other hand, the initial field fluctuations themselves strongly evolve due to free streaming on scales smaller than the free streaming length $\lambda_{\rm fs}$ in Eq. (3).

As alluded to through visualizations of simulations in the previous section, the resolution of this apparent paradox is that the effective number density fluctuations constructed from $\phi^{2}({\bf x})$ with different field momenta should instead be considered as a phase space number density fluctuation, just as it would be for particle dark matter.

For the case considered in the previous section where the field fluctuations at some $q\gg am$ are relativistic, the correspondence to relativistic particles or classical waves is direct. In the classical limit of high photon occupancy, the radiation associated with the photons would be characterized by its electric field ${\bf E}({\bf x})$ and the power in radiation by $|{\bf E}|^{2}({\bf x})$, analogously to $\phi^{2}({\bf x})$. Despite the fact that the electric field contributed by electromagnetic waves of different momenta $q$ always superimpose, observational quantities like the specific intensity do not carry the cross terms of different $q$. The reason is that the cross terms average away over many cycles of their respective oscillations. Moreover, the two-point correlation of the specific intensity $\langle|{\bf E}|^{2}({\bf x})\rangle\langle|{\bf E}|^{2}({\bf x}^{\prime})\rangle$ does not carry either the time averaged power of the individual $q$ modes that $\langle|{\bf E}|^{2}({\bf x})|{\bf E}|^{2}({\bf x}^{\prime})\rangle$ would: the power spectrum of the time average is not the time average of the instantaneous power. Phrased in terms of the highly occupied particle states, the rapidly varying fluctuations in the field that are characterized by the $q$-spectrum represent fluctuation in the fine grained phase space or photon occupancy of momentum states.

While the analogy to photons is direct when the $q$-modes of $\phi$ are relativistic, this absence of a time averaged effect on $\phi^{2}$ is also manifest after the $q$-modes of $\phi$ have become non-relativistic but before equality. Consider the temporal frequency of two non-relativistic modes with $q_{1}\neq q_{2}=|{\bf k}-{\bf q}_{1}|$

$$\omega_{1,2}\approx ma+\frac{1}{2}\frac{q_{1,2}^{2}}{am}.$$ (29)

The cross term contribution to $\phi^{2}$ between the two modes evolves as

$$e^{i\int d\tau(\omega_{1}-\omega_{2})}\approx e^{i\int d\tau\frac{{\bf k}\cdot{\bf q}_{1}}{am}},$$ (30)

where the approximation assumes $k\ll q_{1},q_{2}$. As with $\lambda_{\rm fs}$, this integral grows logarithmically until equality. In fact, unless the wavelength exceeds the free streaming length $k\lambda_{\rm fs}(q_{1};\tau)\ll 1$, the contribution to $\phi^{2}$ of the interference between this pair oscillates in time and would prevent the interference from enabling the growth of dark matter density perturbations below the free streaming scale.

Figure 9: Density evolution and interference of two field momentum modes $q_{1}/a_{*}m=22,25$ that were relativistic at $\tau_{*}$ but nonrelativistic at $\tau\sim 100\tau_{*}$. The beating of these high frequency modes produces density fluctuations on longer scales and larger wavelengths, but the Hubble time averaged density reflects the incoherent sum of the density contributions of the individual modes $\rho_{1}+\rho_{2}$. Normalization is arbitrary.

To see this explicitly, in Fig. 9 we plot the time evolution of the total density as constructed from just two $q$-modes, $\phi=\phi_{1}+\phi_{2}$ with $q_{1}/a_{*}m=22$ and $q_{2}/a_{*}m=25$ using the full growth function (19) for their time evolution from initially equal amplitudes. The total density,

$$\rho=\frac{1}{2}\left(\frac{d\phi}{dt}\right)^{2}+\frac{1}{2a^{2}}(\nabla\phi)^{2}+\frac{m^{2}}{2}\phi^{2},$$ (31)

contains both quadratic combinations from the same $q$ and the cross combination or impact of the superposition of the two. Here $\nabla\phi$ is the spatial field gradient in comoving coordinates, though its impact on $\rho$ is subdominant for these non-relativistic modes where $\langle(d\phi/dt)^{2}\rangle\sim m^{2}\langle\phi^{2}\rangle$ on time average. We compare this to the quadrature sum of the individual modes $\rho_{1}+\rho_{2}$ that omits the superposition term. Even though each individual field mode oscillates with time according to Eq. (29), with phase evolution given by $\omega\tau\sim ma\tau=(\tau/\tau_{*})^{2}\sim 10^{4}$, $\rho_{1}+\rho_{2}$ only evolves on the Hubble time scale since energy is covariantly conserved between the kinetic and potential terms of each term. On the other hand, the superposition causes a beat contribution that oscillates faster than the Hubble time scale and much slower than the mass scale, but time averages to $\langle\rho\rangle_{\tau}=\rho_{1}+\rho_{2}$. Moreover, fluctuations for a given beat frequency ${\bf k}$ are composed of all possible pairs of high frequency momenta that satisfy ${\bf k}={\bf q}_{1}+{\bf q}_{2}$ and each pair contributes with a random temporal phase.

After equality $a\gg a_{\rm eq}$ and for $a\gg a_{*}$, we can remove the fast but $q$ independent mass scale oscillations in Eq. (29) by recasting the scalar field with the Schrödinger wavefunction $\psi$

$$\phi=\frac{1}{\sqrt{2}}(\psi e^{-imt}+\psi^{*}e^{imt})$$ (32)

as done with simulations of ultralight dark matter [119]. In this case, the temporal oscillations between $q_{1}\neq q_{2}$ components of $\psi$ are slow but the spatial phases embedded in $\psi$ are still incoherent.

Well above the de-Broglie scale where $k\ll q$, the field $\psi({\bf x})$ encodes the full phase space of the collisionless dark matter through the Husimi representation. This effectively assigns the spatial variation induced by $q$ on scales that are smaller than some spatial smoothing scale $\eta$ to the momentum distribution at ${\bf p}={\bf q}$:

$$\Psi({\bf x},{\bf p}) = \left(\frac{1}{2\pi}\right)^{3/2}\left(\frac{1}{\pi\eta^{2}}\right)^{3/4}\int d^{3}r\psi({\bf r})$$

$$\times\exp\left(-\frac{({\bf x}-{\bf r})^{2}}{2\eta^{2}}-i[{\bf p}\cdot({\bf r}-{\bf x}/2)]\right).$$

The analog of the phase space distribution function

$$f({\bf x},{\bf p})\equiv|\Psi({\bf x},{\bf p})|^{2}$$ (34)

obeys the collisionless Boltzmann equation. In particular, the directional variation of contributions from incoherent sources visualized in the previous section becomes an anisotropy in the phase space distribution $f({\bf x},{\bf p})=f({\bf p})$ rather than an inhomogeneity in the spatial distribution $f({\bf x},{\bf p})=f({\bf x})$. This construction is similar to the WKB approximation used in Ref. [104] to derive the transfer function (15) in that the small scale variations in $\psi$ are assigned to a spatially varying phase, which carries the phase space momentum, whereas the large scale variations in amplitude carry the positional phase space information.

We can see this explicitly by reversing the construction and embedding a particle phase space distribution $f_{p}$ into the wave $\psi$ in the same way as the random field simulations on a pixelized grid of the previous section. In general, given the discrete Fourier transform momenta indexed by $i$, one assigns $\psi$ at grid points indexed by $\iota$ as [119]

$$\psi({\bf x}_{\iota})=\sum_{i}\sqrt{f_{p}({\bf x}_{\iota},{\bf p}_{i})}e^{i{\bf x}_{\iota}\cdot{\bf p}_{i}}e^{i\alpha_{i}},$$ (35)

where $\alpha_{i}$ is a random spatial phase for each momentum. By explicit substitution we can verify that the Husimi construction (IV) returns the phase space distribution $f\approx f_{p}$ since the cross terms between momenta $i\neq j$ average away due to spatial phase incoherence. Notice that the effective averaging in Eq. (IV) occurs before the squaring in Eq. (34) and is analogous to the temporal averaging considered in Fig. 9 for incoherent temporal phases. These cross terms are negligible so long as the spatial scale of interest is much longer than the smoothing scale $\eta$ and $f_{p}$ is smooth or averaged over momenta scales $1/\eta$.

More specifically in our case where the real space number density fluctuation of interest is on a scale $k\ll q$, we can take $f({\bf x},{\bf p})=f({\bf p})$ and carry the field fluctuations on small scales that enter into the effective number density as momentum space fluctuations only. As with particle dark matter, the free streaming of waves produces a random distribution of momentum streams at each spatial point rather than just spatial fluctuations of a cold distribution. Since CDM isocurvature fluctuations grow in the matter dominated epoch, this would cause the relative transfer function to evolve even if the phase space density fluctuations are conserved. For Schrödinger-Poisson simulation based tests of this construction, see Ref. [120].

V Discussion

We have elucidated the relationship between the free streaming of particle dark matter and wave dark matter and shown how to map the properties of the former onto the latter.

For axion wave dark matter where Peccei-Quinn symmetry breaking occurs after inflation, axion field fluctuations behave like a warm component of particle dark matter in the sense of possessing a mildly relativistic wave spectrum originated from misalignment and axion string emission. Correspondingly, the free streaming length and its impact on curvature fluctuations is only larger than that of cold axions from the pre-inflationary scenario by a logarithmic factor.

We quantify these scaling in terms of the free streaming scale as a function of the characteristic momentum, $\lambda_{\rm fs}(q_{*}=a_{*}m)$, that corresponds to the horizon wavenumber when axions begin their oscillations $H(a_{*})=m$, and compare this to the cold axion Jeans scale where the free streaming of wave fluctuations from curvature perturbations overtake their own wavelength. For axions, free streaming bounds on cold axions or fuzzy dark matter can be roughly translated to the warm case with these scaling relations. However, for warm axions from the post-inflationary scenario, isocurvature fluctuations from the random misalignment on the Hubble scale at $a_{*}$ provide the stronger bound.

If wave dark matter is born ultrarelativistic, then free streaming can have a larger effect as with “hot” dark matter. We provide closed form expressions for the free streaming length $\lambda_{\rm fs}(q)$ for an arbitrary momentum in Appendix A that can be used to assess its impact in any given model with its specific momentum distribution. Generally in such scenarios, the isocurvature fluctuations from causal generation in horizon scale patches at birth can also be affected by free streaming. We illustrate the effect on phase coherent patches and show that they also rapidly damp via free streaming, leaving only phase incoherent transient fluctuations from the superposition of waves of different patches in the effective number density. Despite the free streaming damping of these waves, these incoherent effective number density fluctuations are constant and white at late times when all modes are non-relativistic.

However, these effective number density fluctuations are not spatial number density fluctuations in the wave dark matter, but rather the wave analogue of phase space density fluctuations. At a given spatial position, these fluctuations are composed of multiple field momentum streams from the incoherent sources and the impact of free streaming is similar to the directional damping of collisionless particles. While relativistic, the process is the same as the creation of CMB anisotropy out of plasma inhomogeneities before recombination. As the wave momenta become non-relativistic, the process is analogous to the free streaming damping of fluctuations in massive neutrinos.

More specifically, we show that these free streamed effective number density fluctuations do not behave like real space number density fluctuations over a Hubble or dynamical time in the background or spatially averaged on scales much larger than the de Broglie wavelength of the momentum components. Observables that evolve over a dynamical time respond to the time average of these fluctuations. During radiation domination, these fluctuations oscillate at the beat frequency of the combination of field momenta that compose them; during matter domination, the effective or Husimi phase space representation of the wave dark matter explicitly maps them into multiple momentum streams of the phase space, much like warm or hot dark matter.

Therefore, the astrophysical effects of warm or hot fuzzy dark matter isocurvature fluctuations may also differ from those of cold dark matter isocurvature fluctuations below their free streaming length in a manner that depends on the initial momentum distribution and observable in question. Beyond the axion case, where the free streaming scale is relatively small, we leave the evaluation of specific models and observables to future simulation work.

Appendix A Free Streaming Scalings

In Eq. (4), we defined the free streaming length $\lambda_{\rm fs}$ in units of the comoving Hubble length at matter radiation equality $F(q/a_{\rm eq}m,a/a_{\rm eq})=a_{\rm eq}H_{\rm eq}\lambda_{\rm fs}/\sqrt{2}$ of a momentum component $q$ through the free streaming integral in Eq. (3). In $\Lambda$CDM before dark energy domination, we can explicitly evaluate this integral to obtain

$$F(\hat{q},y)=\frac{\hat{q}\left[\mathcal{F}(\varphi(\hat{q},y),\mu(\hat{q}))-\mathcal{F}(\varphi(\hat{q},0),\mu(\hat{q}))\right]}{\left(1+\hat{q}^{2}\right)^{\frac{1}{4}}},$$ (36)

where $\mathcal{F}(\varphi,m)$ is the incomplete elliptical integral of the first kind with arguments

$$\varphi(\hat{q},y) = \arccos\left(\frac{\sqrt{1+\hat{q}^{2}}-1-y}{\sqrt{1+\hat{q}^{2}}+1+y}\right),$$

$$\mu(\hat{q}) = \frac{1}{2}\left(1+\sqrt{\frac{1}{1+\hat{q}^{2}}}\right).$$ (37)

Although the exact free streaming integral is used in all numerical computations in this work, it is useful to examine the approximate scaling behavior of this solution in various regimes of interest and provide a global approximation that is simple to evaluate.

For $y=a/a_{\rm eq}\ll 1$ the integral is manifestly simple to evaluate and becomes

$$F(\hat{q},y)\approx\hat{q}\sinh^{-1}(y/\hat{q}),\qquad y\ll 1.$$ (38)

Notice that for $\hat{q}\gg y$, the wave is ultrarelativistic $q/am\gg 1$, $F\rightarrow y$, and the free streaming length is the horizon length $\lambda_{\rm fs}=\tau$ as expected. This growth continues until $a\sim q/m$ and thereafter $\lambda_{\rm fs}$ grows logarithmically from its value at $\tau|_{a=q/m}$,

$$F(\hat{q},y)\approx\hat{q}\ln\left(2y/\hat{q}\right),\quad\hat{q}\ll y\ll 1.$$ (39)

This logarithmic growth halts at matter-radiation equality and brings the free streaming scale for modes that are non-relativistic at equality to

$$F(\hat{q},y)\approx\hat{q}\left(\ln\frac{8}{\hat{q}}-\frac{2}{\sqrt{y}}\right),\quad\hat{q}\ll 1\ll y,$$ (40)

which we provided in Eq. (II) in its leading order ($y\gg 1$) form. For axions this limit (40) applies to essentially the entirety of the number density spectrum as we have explicitly verified by comparing its use to the full expression (36) to calculate the transfer function as in Fig. 1.

For evaluating the small contribution from axion momenta that are still relativistic at equality $\hat{q}>1$ or in more general models, it is useful to examine late time contributions to the free streaming integral. Here the free streaming length continues to grow as $\lambda_{\rm fs}\approx\tau\propto a^{1/2}$ during matter domination until $a=q/m$. Taking the matter only scaling for $H(a)=H_{0}\Omega_{m}^{1/2}a^{-3/2}$ in the integral (3), we find

$$\lambda_{\rm fs}=\frac{2a^{1/2}}{H_{0}\Omega_{m}^{1/2}}G(q/am),\quad\hat{q}\gg 1$$ (41)

where

$$G(f) = {}_{2}F_{1}[1/4,1/2,5/4,-f^{-2}]$$

$$\approx \frac{K_{1/2}f^{1/2}}{[1+10f^{1/2}/3+K_{1/2}^{4}f^{2}]^{1/4}},$$

where ${}_{2}F_{1}$ is a hypergeometric function, with $K_{1/2}\approx 1.854$ for the complete elliptic integral of the first kind $K_{m}$ and the approximation is good to a few percent for all $f$. Notice that the $q\gg am$ limit again returns $\lambda_{\rm fs}=2a^{1/2}/H_{0}\Omega_{m}^{1/2}=\tau$ as expected. On the other hand for modes that have become non-relativistic between $a_{\rm eq}$ and $a$, the free-streaming length approaches a constant value

$$\lambda_{\rm fs} \approx 2K_{1/2}\frac{\sqrt{q/m}}{H_{0}\Omega_{m}^{1/2}},\quad a_{\rm eq}\ll q/m\ll a,$$ (43)

associated with the horizon length at the epoch the momentum becomes non-relativistic.

In summary, to a few percent accuracy we can approximate the whole solution Eq. (36) in the $y\gg 1$ regime by joining these approximations

$$F(\hat{q},y)=\begin{cases}\hat{q}\ln(\frac{8}{\hat{q}}),&\hat{q}<1\\ 2y^{1/2}[G(\frac{\hat{q}}{y})-G(\frac{1}{y})]+\ln 8,&\hat{q}\geq 1\end{cases}$$ (44)

and using the simple approximation for $G$ in Eq. (A) such that all of the various scaling regimes are manifest.

Figure 10: Approximations for the relative transfer function (15) [$T_{\rm rel}$ warm], using $k_{\rm fs}=\lambda_{\rm fs}^{-1}(q_{*})$ in Eq. (45) [$T_{\rm rel}^{k_{\rm fs}}$], and an approximation from Ref. [104] v2 using $k_{\rm AM}$ [$T_{\rm rel}^{k_{\rm AM}}$] of Eq. (47) instead. The former captures the scale at which free streaming occurs whereas the latter changes this scale to be orders of magnitude smaller in $k$ and would cause this mass $m=10^{-14}$eV to be inappropriately ruled out.

Appendix B Averaging over Momenta

In the main text Eq. (15), we averaged the effect of free streaming over the momentum distribution of the axions $q^{3}P_{\phi}(q)$ to approximate the net impact on the transfer function. For the spectrum of Eq. (10) which has a sharp peak at $q_{*}$, the impact is similar to evaluating an effective $k_{\rm fs}=\lambda_{\rm fs}^{-1}(q_{*})$ and taking

$$T_{\rm rel}^{k_{\rm fs}}(k)=\frac{\sin(k/k_{\rm fs})}{k/k_{\rm fs}}$$ (45)

in the region $k\sim k_{\rm fs}$ where the damping starts to have its effect. In Fig. 10, we compare this approximation to Eq. (15). Notice that this approximation does provide the correct scaling for the wavenumbers where free streaming starts to become important but underestimates the effect at higher $k$. Mathematically, this comes about because Eq. (15) is an integration over an oscillating quantity. Even in this case where the spectrum is peaked near $q_{*}$, the phase $k\lambda_{\rm fs}(q)$ varies over an increasingly large range as $k$ increases. Note however that the derivation of Eq. (15) itself in Ref. [104] is not ensured to be physically valid for $k\lambda_{\rm fs}\gg 1$ and should be considered as an estimate for the half power point.

In Ref. [104] v2, this effective $k_{\rm fs}$ approach was adopted but instead of weighting the impact of free streaming by the number density spectrum, they equated the Taylor expansion of Eq. (15)

$$\lim_{k\rightarrow 0}T_{\rm rel}(k)\approx 1-\frac{1}{{6}}\frac{k^{2}}{k_{\rm AM}^{2}}$$ (46)

where

$$\frac{1}{k_{\rm AM}^{2}}=\frac{\int_{0}^{am}d\ln q\frac{q^{3}}{2\pi^{2}}P_{\phi}(q)\lambda_{\rm fs}^{2}(q)}{\int d\ln q\frac{q^{3}}{2\pi^{2}}P_{\varphi}(q)},$$ (47)

to that of Eq. (45)

$$\lim_{k\rightarrow 0}T_{\rm rel}^{k_{\rm fs}}(k)\approx 1-\frac{1}{{6}}\frac{k^{2}}{k_{\rm fs}^{2}}$$ (48)

to imply that $k_{\rm fs}\rightarrow k_{\rm AM}$ in Eq. (45) and $T_{\rm rel}^{k_{\rm fs}}\rightarrow T_{\rm rel}^{k_{\rm AM}}$ (their Eq. 5 v2). In Fig. 10, we compare this transfer function to Eq. (15) and Eq. (45). Notice that this weighting scheme overestimates the effect of free streaming by three orders of magnitude, with $k_{\rm fs}\approx 3.2\times 10^{3}$ Mpc^-1 while $k_{\rm AM}\approx 3.1$ Mpc^-1. The overestimate is so large that this $m=10^{-14}$eV example would be inappropriately ruled out. In Ref. [104] v2, this approximation was used to place a bound of $m>2\times 10^{-12}$eV for the spectrum considered here and $m>10^{-12}$eV for their axion parameterization with mildly relativistic modes from string decay.

Figure 11: Relative transfer function for the warm axions at the same $m=10^{-14}$ eV mass as Fig 10 but with variations to the power spectrum which respectively increase the peak momentum by $R=1,2,4$ and alter the high momentum slope $\alpha=3/4,1,5/4$ using Eq. (49). We see that the variation of $R$ shifts the scale of the free streaming damping for the same $\alpha=1$, whereas altering $\alpha$ for the same $R=1$ only makes small changes in the shape of the damping.

This overestimate is tied to the behavior of the high momentum tail with the spectrum in Eq. (10). In Eq. (43) we show that for waves that become nonrelativistic in the matter dominated regime, $\lambda_{\rm fs}\propto q^{1/2}$. Thus for the spectrum $q^{3}P_{\phi}\propto q^{-1}$, the integral in Eq. (47) receives nearly constant contributions per $d\ln q$ up to the $q\approx am$ limit where the waves are still relativistic at the evaluation epoch despite the highly suppressed number of axions with such momenta. The result is that the estimate of the effective free streaming wavenumber $k_{\rm AM}$ produces a suppression of the transfer function to much smaller wavenumbers or much larger scales than calculated from Eq. (15) or estimated by Eq. (45). Mathematically, the Taylor expansion (46) is not a good approximation at $k\sim k_{\rm AM}$ since the dominant momenta near $q_{*}$ have vanishingly small free streaming effect as $k\rightarrow 0$, causing the second derivative of $T_{\rm rel}$ to run strongly with scale.

The source of this discrepancy is the difference in the weighting scheme. Since $\lambda_{\rm fs}$ grows as $\tau$ for relativistic momenta, the weighting in Eq. (47) allows a very small number density in high momentum waves to dominate the effective free streaming length $k_{\rm AM}$ of the whole population, whereas physically free streaming implies that instead this small component is smooth across scales where the dominant component remains clustered, similar to the small admixture of massive neutrinos and cold dark matter in $\Lambda$CDM. That both momenta can be represented by the single field $\phi({\bf x})$ is also related to the Husimi phase space construction discussed in §IV. The spatially smooth and clustered components are embedded in the inferred momentum distribution.

Since the key quantity that controls the free streaming effect is the shape of the momentum distribution, we also explore variations from Eq. (10) that adjust the position of the peak in the spectrum and the power law decline from the peak, parameterized by $R$ and $\alpha$ as follows:

$$q^{3}P_{\phi}(q)\propto\left(\frac{q}{Rq_{*}}\right)^{3}\theta(Rq_{*}-q)+\left(\frac{Rq_{*}}{q}\right)^{\alpha}\theta(q-Rq_{*}).$$ (49)

Fig. (11) shows the corresponding change in the transfer function. In the top panel, we fix $\alpha=1$ and increase $R$ from $1$ to $2$ and $4$. The damping wavelength increases nearly linearly with $R$ in accordance with the expectation that the free streaming length scales as $\lambda_{\rm fs}(Rq_{*})$ discussed in §III. Varying $\alpha$ in the range where most of the particles still have momenta $\sim q_{*}$ has a much smaller effect since only the small tail of high momenta waves are affected. These variations in $\alpha$ encompass the full range found in current state of the art axion simulations [109, 110, 111].