Evolution of linear wave dark matter perturbations in the radiation-dominant era

TL;DR

This paper analyzes the evolution of wave dark matter perturbations during the radiation era, revealing distinct phases and quantum effects that cause deviations from cold dark matter predictions, with implications for axion models.

Contribution

It provides a detailed analysis of $\psi$DM perturbations in the radiation era, identifying four evolutionary phases and quantum effects causing mode suppression, and compares axion models with $\psi$DM.

Findings

Late-stage $\psi$DM perturbations resemble CDM but with subtle differences.

Quantum effects lead to significant mode suppression in $\psi$DM.

Axion models with near-top initial angles predict higher power spectra and cutoff frequencies.

Abstract

Linear perturbations of the wave dark matter, or $\psi$ dark matter ($\psi$DM), of particle mass $\sim 10^{-22}$eV in the radiation-dominant era are analyzed, and the matter power spectrum at the photon-matter equality is obtained. We identify four phases of evolution for $\psi$DM perturbations, where the dynamics can be vastly different from the counterparts of cold dark matter (CDM). While in late stages after mass oscillation long-wave $\psi$DM perturbations are almost identical to CDM perturbations, some subtle differences remain, let alone intermediate-to-short waves that bear no resemblance with those of CDM throughout the whole evolutionary history. The dissimilarity is due to quantum mechanical effects which lead to severe mode suppression. We also discuss the axion model with a cosine field potential. The power spectrum of axion models are generally almost identical to those of $\psi$DM, but in the extreme case when the initial axion angle is near the field potential top, this axion model predict a power excess over a range of wave number and a higher spectral cutoff than $\psi$DM as if $\psi$DM had a higher particle mass.

Equations

$${\Psi }^{{}^{\prime \prime }}+2H{\Psi }^{{}^{\prime }}+m^2{a}^2\Psi =0,$$

$${\psi }^{{}^{\prime \prime }}+2H{\psi }^{{}^{\prime }}+(k^2+m^2{a}^2)\psi =4{\Psi }^{{}^{\prime }}{\varphi }^{{}^{\prime }}-2m^2{a}^2\Psi \varphi .$$

$${\delta }_{\gamma }^{{}^{\prime }}-\frac{4}{3}(k^2{\theta }_{\gamma }+3{\varphi }^{{}^{\prime }})=0,$$

$${\theta }_{\gamma }^{{}^{\prime }}=-\frac{{\delta }_{\gamma }}{4}-\varphi ,$$

$$\begin{array}{rl}& -k^2\varphi -3H({\varphi }^{{}^{\prime }}+H\varphi )=\\ & 4\pi G\{[{\Psi }^{{}^{\prime }}{\psi }^{{}^{\prime }}+m^2{a}^2\Psi \psi -({\Psi }^{{}^{\prime }}{)}^2\varphi ]+a^2{ϵ}_{\gamma }{\delta }_{\gamma }\},\end{array}$$

$${\varphi }^{{}^{\prime }}+H\varphi =4\pi G[{\Psi }^{{}^{\prime }}\psi -a^2({ϵ}_{\gamma }+P_{\gamma }){\theta }_{\gamma }].$$

$$a^2{ϵ}_{\psi }{\delta }_{\psi }={\Psi }^{{}^{\prime }}{\psi }^{{}^{\prime }}+m^2{a}^2\Psi \psi -({\Psi }^{{}^{\prime }}{)}^2\varphi$$

$${ϵ}_{\psi }\equiv [({\Psi }^{\prime }{)}^2+m^2{a}^2{\Psi }^2]/(2a^2).$$

$${ϵ}_{\alpha }{\Delta }_{\alpha }={ϵ}_{\alpha }{\delta }_{\alpha }-3H(P_{\alpha }+{ϵ}_{\alpha }){\theta }_{\alpha },\alpha =\gamma ,\psi ,$$

$${\theta }_{\psi }=-\frac{{\Psi }^{{}^{\prime }}\psi }{a^2(P_{\psi }+{ϵ}_{\psi })},$$

$$P_{\psi }\equiv [({\Psi }^{{}^{\prime }}{)}^2-m^2{a}^2{\Psi }^2]/(2a^2),$$

$$-k^2\varphi =4\pi G{a}^2({ϵ}_{\psi }{\Delta }_{\psi }+{ϵ}_{\gamma }{\Delta }_{\gamma }).$$

$$\delta P_{\psi }=a^{-2}{\Psi }^{{}^{\prime }}{\psi }^{{}^{\prime }}-m^2\Psi \psi -a^{-2}({\Psi }^{{}^{\prime }}{)}^2\varphi ,$$

$${\delta }_{\psi }^{{}^{\prime }}+[ln(a^3{ϵ}_{\psi }){]}^{{}^{\prime }}{\delta }_{\psi }+3H\frac{\delta P_{\psi }}{{ϵ}_{\psi }}-(1+\frac{P_{\psi }}{{ϵ}_{\psi }})(k^2{\theta }_{\psi }+3{\varphi }^{{}^{\prime }})=0.$$

$${\theta }_{\psi }^{{}^{\prime }}+\{ln[(P_{\psi }+{ϵ}_{\psi })a^4]{\}}^{{}^{\prime }}{\theta }_{\psi }=-\frac{\delta P_{\psi }}{P_{\psi }+{ϵ}_{\psi }}-\varphi .$$

$$-i{\hat{\psi }}^{{}^{\prime }}+\frac{k^2}{2ma}\hat{\psi }=-ma\varphi .$$

$${\delta }_{\psi }=2\Re [\hat{\psi }]+\frac{3H}{ma}[2{\sin }^2(mt+\alpha )-1]\Im [\hat{\psi }]$$

$$\begin{array}{rl}{\Delta }_{\psi }& ={\delta }_{\psi }-\frac{6H}{ma}\Im [\hat{\psi }]{\sin }^2(mt+\alpha )\\ & =2\Re [\hat{\psi }]-\frac{3H}{ma}\Im [\hat{\psi }].\end{array}$$

$$\frac{⟨{\delta }_i⟩}{⟨1+w_i⟩}=\frac{⟨{\delta }_j⟩}{⟨1+w_j⟩}$$

$${\Delta }_{\psi }=-\frac{3H}{ma}\Im [\hat{\psi }{]}_h.$$

$$⟨\delta P_{\psi }⟩=\frac{k^2}{2m^2{a}^2}{ϵ}_{\psi }(\Re [\hat{\psi }])\approx \frac{k^2}{4m^2{a}^2}{ϵ}_{\psi }⟨{\delta }_{\psi }⟩.$$

$$-i\frac{d{\hat{\psi }}_h}{d\eta }+\frac{k^2}{4H_m^2}{\hat{\psi }}_h=0.$$

$$\left.\begin{aligned} \Delta_{\psi}\sim|g|\Big{\{}&2\sin\Big{[}{{k^{2}}\over{4H_{m}^{2}}}(\eta-\eta_{k})\Big{]}+\\ &{{3H}\over{ma}}\cos\Big{[}{{k^{2}}\over{4H_{m}^{2}}}(\eta-\eta_{k})\Big{]}\Big{\}}.\end{aligned}\right.$$

$${\psi }^{{}^{\prime \prime }}+2H{\psi }^{{}^{\prime }}=-2m^2{a}^2\Psi \varphi ,$$

$${\psi }^{{}^{\prime \prime }}+k^2\psi =4{\Psi }^{{}^{\prime }}{\varphi }^{{}^{\prime }}$$

$${\Delta }_{\nu }=\{\begin{array}{ll}{\Delta }_{ph},& \tau \le {\tau }_1,\\ 0,& \tau >{\tau }_1,\end{array}$$

$$w_{ph}\equiv{{P_{ph}+P_{b}}\over{\epsilon_{ph}+\epsilon_{b}}}={{a_{eq}}\over{3\Big{(}a_{eq}+a{{\Omega_{b}}\over{\Omega_{m}}}{{\Omega_{ph}+\Omega_{\nu}}\over{\Omega_{ph}}}\Big{)}}}.$$

$${\delta }_{ph}^{{}^{\prime }}-\frac{4}{3}(k^2{\theta }_{ph}+3{\varphi }^{{}^{\prime }})=0,$$

$${\theta }_{ph}^{{}^{\prime }}=-\frac{\delta P_{ph}}{{ϵ}_{ph}+P_{ph}}-\varphi +a{n}_e{\sigma }_T({\theta }_b-{\theta }_{ph}),$$

$${\delta }_b^{{}^{\prime }}=k^2{\theta }_b+3{\varphi }^{{}^{\prime }},$$

Full text

I Introduction

A new form of dark matter, so-called $\psi$DM, wave dark matter or fuzzy dark matter, that consists of extremely light particles of typical particle mass, $10^{-22}$eV Hu et al. (2000); Schive et al. (2014); Marsh and Pop (2015); Calabrese and Spergel (2016); Hui et al. (2016), has been noted to exhibit nice, yet peculiar, features for explaining recent observational results. For example, it produces finite density cores for Milky Way dwarf spheroid galaxies to reside in Schive et al. (2014); Lora and Magaña (2014), an observational fact Moore (1994); Goerdt et al. (2006); Gilmore et al. (2007); Walker and Peñarrubia (2011); Amorisco et al. (2013) that has been puzzling cold dark matter (CDM) proponents for some time. This core in the $\psi$DM scenario results from the uncertainty principle of these extremely light particles, which yields quantum pressure and avoids the central density singularity predicted by CDM Navarro et al. (1997). Other peculiar features of $\psi$DM arise from the linear matter spectrum described below.

The $\psi$DM spectrum is known to be such that it resembles CDM spectrum in long wavelengths and becomes heavily suppressed in short wavelengths Hu et al. (2000). The spectral suppression, due also to the quantum pressure, not only occurs already in the radiation-dominant era but occurs sharply at some transition wavenumber, much sharper than the suppression of the short-wavelength spectrum of warm dark matter (WDM) Bode et al. (2001). This spectral feature on one hand leads to paucity of Milky Way satellite galaxies Klypin et al. (1999); Moore et al. (1999), which is closely related to the too-big-to-fail problemBoylan-Kolchin et al. (2011); Papastergis et al. (2015), and on the other hand postpones first galaxy formation thereby resulting in delay reionization Marsh, and Silk (2014); Schive et al. (2016). These nice features are not attainable for the WDM model Macciò et al. (2012). The transition wavenumber, which we call the critical wave number $k_{c}$, is generally smaller than, but of the same order of, the slowly evolving Jeans wave number $k_{J}$ of $\psi$DM in the matter-dominant era Hu et al. (2000); Woo, and Chiueh (2009). It is therefore of importance to understand how such a spectral transition occurs dynamically and what makes the connection between $k_{c}$ and $k_{J}$.

The evolution of the linear $\psi$DM perturbation in the matter-dominant era has been analyzed in Woo, and Chiueh (2009), in which it gave insights as to why the long-wavelength power resembles to the CDM model and how the short-wave power is suppressed due to the presence of quantum Jeans length. While the perturbation dynamics of matter-dominant regime is straightforward, the dynamics in the radiation-dominant regime is much more complex and must be solved numerically Hlozek et al. (2015); Marsh (2016). Unfortunately, the numerical solution offers limited insights to the dynamics. A fluid approach Marsh (2016) for sub-horizon perturbations has also been explored, which can be simple enough to permit detailed analyses and sheds some lights on the perturbation dynamics. However, important dynamics that shapes the final spectrum turns out to occur when perturbations are super-horizon where the fluid approach fails. Hence the fluid approach can only partially answer the questions we intend to address.

Note that $\psi$DM consists of many free bosons condensed in a quantum ground state, a Bose-Einstein condensate, for which they remain phase coherent over an astronomical distance to share the same wave function. However, there is a problem. Lacking causal contacts, free bosons are not capable of achieving phase coherence. Nonlinearity is needed to make these bosons interact and become phase-locked. In this regard, the axion mechanism offers a plausible solution to the phase-locking problem, for example Chiueh (2014); Hui et al. (2016). The axion model can provide strong boson coupling early on and sets free these bosons at a later time; in other words, the phase-lock mechanism takes place in the very early time and since then these bosons remain phase-coherent even after they become free. Due to this connection to the initial condition of $\psi$DM, it is also important to extend the free-particle model of $\psi$DM and examine the axion model.

The organization of this paper is as follows. Sec. (II) introduces the relevant equations for later analyses. We also give the fluid formulation from field equations and show that it generally cannot be evolved. Sec. (III) contains our main results, where four asymptotic phases are identified, analyzed and approximate solutions obtained. We make clear the distinction between the passive evolution and the evolution of full treatment in Sec. (IV), where solutions of both will be presented and compared. We identify a special solution common at the boundary of the four asymptotic phases and define a critical wavenumber $k_{c}$ which marks the spectral transition in Sec. (V). The adiabatic condition for super-horizon modes is examined in Sec. (VI) and found to be valid instantaneously for this $\psi$DM problem without the need of time average over the fast mass oscillation. We examine the axion model and identify its similarities to and differences from the free-particle model in Sec. (VII). Conclusions are made in Sec. (VIII). In Appendix (A), we give full solution of passive evolution. In contrast to the main text in which we obtain solutions to the approximate equations in various asymptotic phases, we derive in Appendix (A) approximate solutions to the full solution for the four asymptotic phases. In Appendix (B), we turn to addressing equations and analyses for the full treatment of perturbations, including neutrino, baryon and photon perturbations.

Throughout the paper, we take the fiducial particle mass to be $10^{-22}$eV and adopt the Newtonian gauge for the perturbation. We also let the speed of light, $c$, and the Planck constant $\hbar$ be equal to $1$.

II Governing Equations

Standard derivation of the governing equations is given in Appendix (A). Here, we provide notations and essential equations to be used in this paper.

Let $a_{*}$ and $H_{*}$ be the scaling factor and Hubble parameter at some given epoch in the radiation era, and define $d\tau=dt/a$ and the Hubble parameter $H\equiv d\ln a/d\tau$. Since $H\propto a^{-1}$, we have $Ha=H_{*}a_{*}$, a fixed constant, and it follows $\tau=(H_{*}a_{*})^{-1}a=H^{-1}$.

Moreover, we decompose the $\psi$DM wave function into a time-dependent background $\Psi(\tau)$ and a space-time dependent perturbation $\psi$ as usual. The metric perturbation in the Newtonian gauge is denoted as $\phi$. The quantities $\delta_{\gamma}$ and $\theta_{\gamma}$ are the dimensionless energy density perturbation and the velocity potential of the photon fluid. These perturbed quantities in the comoving coordinate are Fourier transformed in space into plane-wave eigenmodes with the comoving wavenumber $k$ as eigenvalues.

Denote the prime to be $d/d\tau$. The zeroth-order field $\Psi$ obeys

[TABLE]

and the perturbed field $\psi$ obeys

[TABLE]

The radiation perturbation equations are

[TABLE]

and

[TABLE]

where $\delta_{\gamma}\equiv\delta\epsilon_{\gamma}/\epsilon_{\gamma}$ with $\epsilon_{\gamma}$ being the energy density of the radiation fluid, and $\theta_{\gamma}$ is the perturbed velocity potential of the radiation fluid. Neutrinos can be regarded as a part of the radiation fluid for super-horizon perturbations. But upon entering horizon, neutrino perturbations die out sharply as a result of their collisionless nature where collisionless damping prevails. In Sec. (IV) we have more discussions on this issue.

The equations for the metric perturbation $\phi$ read

[TABLE]

and

[TABLE]

One can identify the perturbed energy density as

[TABLE]

on the right-hand side of Eq. (5), where

[TABLE]

One can further substitute Eq. (6) into Eq. (5) to obtain a simplified equation. Recognizing the gauge covariant energy perturbation

[TABLE]

and the momentum potential

[TABLE]

with

[TABLE]

it follows that Eq. (5) becomes

[TABLE]

The Poisson equation with gauge covariant sources is recovered.

Although the fluid description of $\psi$DM is not useful in general, it is illuminating to find out its difference from that of CDM. We multiply Eq. (2) by $a\Psi^{{}^{\prime}}$ and add Eq. (1) multiplied by $a\psi^{{}^{\prime}}$. As the perturbed pressure as

[TABLE]

we find the perturbed energy equation

[TABLE]

For the perturbed momentum equation, we multiply Eq. (1) by $\psi$ and obtain

[TABLE]

Equations (14) and (15) are actually quite general and valid for any standard field potential. They become the radiation fluid perturbations when $(\delta P,P)=(1/3)(\epsilon\delta,\epsilon)$ with $\epsilon\propto a^{-4}$ and also the CDM perturbations when $\delta P=P=0$ with $\epsilon\propto a^{-3}$, where the ”$\psi$” index is replaced by the photon index and the CDM index, respectively. Note that Eqs. (14) and (15) need equations of state to relate $P$ to $\epsilon$ and $\delta P$ to $\epsilon\delta$. These relations are readily available for standard fluids, but generally far from trivial here. In fact, one must solve the zero-order and perturbed field equations, Eqs. (1) and (2), to construct $P_{\psi}$ and $\delta P_{\psi}$ expressed in terms of $\epsilon_{\psi}$ and $\delta\epsilon_{\psi}$. Therefore without solving for the field, the fluid description of the field is generally of no practical use.

III Passive Evolution and Asymptotic Solutions

Passive evolution of dark matter can be a good approximation of the dark matter dynamics in the radiation-dominant era, where the metric perturbation is governed entirely by the radiation fluid, which includes neutrinos. Neutrinos are actually collisionless particles, and their perturbations behave identical to radiation perturbations before entering horizon but die out rapidly once they enter horizon Ma, and Bertschinger (1995). We approximate neutrinos always as a part of the radiation fluid for passive evolution, and will show in the next section that this crude treatment is still a good approximation for dark matter perturbations of main interest.

The right-hand side of Eq. (2) can thus be approximated as external sources independent of $\psi$ for passive evolution. This assumption greatly simplifies the following asymptotic analyses, bringing out the details of the differences between the CDM perturbation and the $\psi$DM perturbation caused by the introduction of a finite particle mass $m$ in Eq. (2). In the limit of $m\to\infty$, the $\psi$DM perturbation should recover the CDM result. With this approach, analytical expressions for source terms in Eq. (2) can actually be obtained and therefore the solution $\psi$ can be integrated by the Green’s function method. The full solution however involves confluent hypergeometric functions and is not illuminating to be given in the main text, and hence we present it in Appendix (A).

Plotted in Fig. (1) is the passive time evolutions of two distinct modes with mode numbers $k=0.1k_{c}$ and $k=10k_{c}$, where $k_{c}$ is the wave number of the critical mode that enters horizon at the onset of mass oscillation. (More discussions on $k_{c}$ are given in Sec. (V).) There exist four asymptotic phases of evolution as shown in Fig. (1), and this plot sets the stage for discussions to follow.

The evolution of $\psi$DM perturbations has a natural dividing line, the onset of mass oscillation. Prior to the mass oscillation, $\psi$DM is like an inflaton and after that, $\psi$DM is like CDM influenced to various degrees by the additional quantum pressure. Since the latter is closely resembled the familiar CDM, we shall discuss the second phase first. Little surprise arises in the second phase, given that we know the perturbed $\psi$DM dynamics in the matter-dominant era Woo, and Chiueh (2009). The new features of $\psi$DM is actually given by the earlier phase prior to mass oscillation, resulting in a sharp transition of the $\psi$DM spectrum.

(i) After mass oscillation $2H,k\ll ma$

For long wave ($k\ll ma$) perturbations well after the onset of mass oscillation, $ma\gg 2H$, the dynamical system is like a driven damped oscillator near resonance, with the main frequency $m$ and detuning frequency $k$. Equation (2) can be simplified by extracting the common mass oscillation factor in all terms and is approximated by a complex amplitude equation, for which the detuning frequency $k$ appears in the equation. The simplification procedure goes as follows. Let $\psi\equiv A\Psi_{r}+B\Psi_{s}$, where $\Psi_{r}$, the regular solution of $\Psi$, and $\Psi_{s}$, the singular solution that diverges at $\tau\to 0$, are given in Appendix (A), $A$ and $B$ are real fields with slow time variations, and the zero-order field $\Psi=\Psi_{r}$. To the accuracy of $O(2H/ma)$, $(\Psi_{r},\Psi_{r}^{\prime}/ma)\sim(2\epsilon_{\psi}/m^{2})^{1/2}(\cos(mt-q+\alpha)+O(H^{2}/m^{2}a^{2}),-\sin(mt+7q+\alpha)+O(H^{2}/m^{2}a^{2}))$, with $\alpha$ being a constant phase and $q=3H/16ma\propto 1/t$, and $(\Psi_{s},\Psi_{s}^{\prime}/ma)\sim(2\epsilon_{\psi}/m^{2})^{1/2}(\sin(mt-q+\alpha)+O(H^{2}/m^{2}a^{2}),\cos(mt+7q+\alpha)+O(H^{2}/m^{2}a^{2}))$. One may further define a complex $\hat{\psi}$, where $\hat{\psi}=A+iB$, and therefore $\psi=Re[\hat{\psi}(\Psi_{r}-i\Psi_{s})]\to(2\epsilon_{\psi}/m^{2})^{1/2}\Re[\hat{\psi}\exp(-i(mt+\alpha))]$ in the limit $H/ma\to 0$.

The amplitude $\hat{\psi}$ equation can be straightforwardly derived from Eq. (2),

[TABLE]

This equation is the familiar linearized Schrödinger equation for the perturbed field.

Note that in order for Eq. (16) to be correctly derived, it is essential to keep the small $q$ in $\Psi_{r}$ and $\Psi_{s}$. Otherwise, there would be additional incorrect terms of order $O(H^{2}\hat{\psi})$ appearing in Eq. (16), which would have otherwise given an erroneous solution for Phase (ia).

Eq. (16) is driven by a source proportional to $\phi$, and the solution consists of a particular integral and a homogeneous solution. The adiabatic perturbation, which is of interest to us, corresponds to the particular integral of the full equation, Eq. (2). For the asymptotic regimes to be discussed below, certain limits have been taken to simplify the full equation, but this simplification comes with a cost, in decomposing the particular integral of the full equation into the particular integral and the homogeneous solution of the simplified equation. The homogeneous solution can be regarded as the initial condition of the solution to the simplified equation and assumes a undetermined amplitude, which needs to be fixed by some means, for example, matching to the solution in the previous phase.

(ia) $k\ll 2H\ll ma$

In this regime, the term proportional to $k^{2}$ in Eq. (16) can be ignored to the leading order. Now, as the super-horizon metrics perturbation $\phi=\phi_{0}[1-(1/30)(k/H)^{2}+O(k^{4}/H^{4})]$ where $\phi_{0}$ is a constant (shown in Appendix (A)), the driving source in the real part of Eq. (16) grows as $a$. Thus the particular integral solution $\Im[\hat{\psi}]_{p}=-(ma/2H)\phi_{0}\propto a^{2}$, and $\Re[\hat{\psi}]_{p}=0$. On the other hand the homogeneous solution, $\hat{\psi}_{h}$, of Eq. (16) is a complex constant.

To the leading order, we have

[TABLE]

(c.f. Eq. (7)), where $\Im[\hat{\psi}]^{{}^{\prime}}=-ma\phi_{0}$ has been used. The quantity $\Re[\hat{\psi}]$ is an arbitrary homogeneous solution $\Re[\hat{\psi}]_{h}$, and $\Im[\hat{\psi}]$ is a combination of the particular integral and a constant homogeneous solution. Finally, the covariant energy perturbation

[TABLE]

Note that while the gauge-dependent $\delta_{\psi}$ oscillates, the gauge-covariant $\Delta_{\psi}$ manifestly does not.

It requires some guidance to fix the homogeneous solution $\hat{\psi}_{h}$ in Eq. (18). Conventionally, matching of solutions in different asymptotic regimes can provide such guidance. However we adopt a different approach here. For super-horizon adiabatic perturbation, $\delta_{\psi}$ obeys the relation:

[TABLE]

for any species $i$ and $j$, where we denote $<...>$ to represent a short-time average to filter out the fast mass oscillation. (However, see Sec. (VI).) It follows $\langle\delta_{\psi}\rangle=(3/4)\delta_{\gamma}$. From Appendix (A), we also know that $\delta_{\gamma}\sim-2\phi_{0}$ and hence $\langle\delta_{\psi}\rangle\sim-(3/2)\phi_{0}$. Therefore the constant $\Re[\hat{\psi}]_{h}=-(3/4)\phi_{0}$. Given $\Im[\hat{\psi}]_{p}=-(ma/2H)\phi_{0}$, we find

[TABLE]

The $\phi$ dependence in $\Delta_{\psi}$ cancels and only the arbitrary homogeneous solution $\Im[\hat{\psi}]_{h}$ survives.

If $\Im[\hat{\psi}]_{h}\sim O(\Re[\hat{\psi}]_{h})$, the remaining term of $\Delta_{\psi}$ is a small quantity of order $2H/ma$ and decays as $a^{-2}$, thus negligible. From the definition of $\hat{\psi}$, the imaginary part $\Im[\hat{\psi}]$ is to be multiplied by $\Psi_{s}$ and the real part $\Re[\hat{\psi}]$ by $\Psi_{r}$ to yield the original field $\psi$. The ratio $\Psi_{s}/\Psi_{r}$ diverges in early epoches when $2H/ma\to\infty$, and the ratio $\Im[\hat{\psi}]/\Re[\hat{\psi}]$ in early time must approach zero. Therefore at the onset of mass oscillation $2H=ma$ immediately before the perturbation enters the present phase, which is the main contribution of the homogeneous solution, the quantity $\Im[\hat{\psi}]_{h}$ is at most of the same order of $\Re[\hat{\psi}]_{h}$. In Appendix (A), we show that $\Im[\hat{\psi}]_{h}$ is in fact of higher-order smallness compared to $\Re[\hat{\psi}]_{h}$. Therefore the gauge-covariant $\Delta_{\psi}$ cancels itself to the leading order ($O(1)\phi_{0}$), and we must consider the next-order contributions to $\Delta_{\psi}$.

In this phase, there are two types of high-order terms, $O(2H/ma)$ and $O((k/H)^{2})$. The former decays as $a^{-2}$ and the latter grows as $a^{2}$. We shall consider the latter. The next order contributions to $\Delta_{\psi}$ arise from all other terms in Eq. (2) but are neglected in Eq.(16).

A straightforward but lengthy calculation by keeping all terms of order $O(k^{2}/H^{2})$ in Eq. (2) reveals that $\Re[\hat{\psi}]_{p}=-(7/40)(k/H)^{2}\phi_{0}$ and $\Im[\hat{\psi}]_{p}=-(ma/2H)\phi_{0}(1-(1/60)(k/H)^{2})$. The homogeneous solution $\Re[\hat{\psi}]_{h}=-(3/4)\phi_{0}(1+O((k/H)^{2}(H/ma)))$ and remains the same as before to the order in question; so is $\Im[\hat{\psi}]_{h}$. Substituting the above findings to the second equality of Eq. (18), we now have $\Delta_{\psi}=-(3/8)(k/H)^{2}\phi_{0}=(9/16)\Delta_{\gamma}\propto a^{2}$, where the Poisson’s equation (Eq. (12)) has been used to bring out $\Delta_{\gamma}$.

The $\psi$DM gauge-covariant energy perturbation in this phase is identical to the CDM counterpart, for which the growth is independent of the particle mass $m$. Note that the perturbed field goes through this phase only for sufficiently low-$k$ modes. Modes with sufficiently high-$k$ skip this phase and directly enter Phase (ib) to be discussed below.

(ib) $2H\ll k\ll ma$

In this regime, the pressure perturbation can be of dynamical importance despite $\langle P_{\psi}\rangle=0$. Here, $\langle\delta P_{\psi}\rangle=(\epsilon_{\psi}/ma)\Re[i\hat{\psi}^{{}^{\prime}}-ma\phi]$ from Eq. (13). Compared with the real part of Eq. (16), one readily recognizes that

[TABLE]

The second equality of Eq. (21) holds because the term proportional to the metric perturbation $\phi$ in $\delta_{\psi}$ is a high-order term and can be ignored. Now that with $\langle\delta P_{\psi}\rangle$ available, the dynamics of perturbation in the sub-horizon regime can be described by fluid equations.

As $\epsilon_{\psi}\langle\delta_{\psi}\rangle\gg\langle\delta P_{\psi}\rangle$, we seem to recover CDM. But it is not so even when $k^{2}\ll m^{2}a^{2}$, since this inequality is only an indication of $\psi$DM becoming non-relativistic. Non-relativistic dark matter can nevertheless have sufficiently large pressure to counter the gravity for short waves. Hence only in the long wave limit of this regime $k^{2}\ll(m^{2}a^{2})(2H/ma)$ can $\psi$DM perturbations resemble CDM perturbations, as we shall see below.

In discussions to follow, we shall stick to the field equation for consistency, despite that the fluid equations are also well-defined. In this regime the metric perturbation $\phi\sim\phi_{0}(H/k)^{2}\cos(k(\tau-\tau_{k})/\sqrt{3})$ as discussed in Appendix (A), where $\tau_{k}$ is the instant for mode $k$ to enter horizon. Therefore the driving source ($\propto\phi$) of Eq. (16) is of order $O(a^{-1})$. It is straightforward to find that the particular integrals $\Re[\hat{\psi}]_{p}\sim\phi_{0}(3/2)(H/k)^{2}\cos(k(\tau-\tau_{k})/\sqrt{3}))$ and $\Im[\hat{\psi}]_{p}\sim-\phi_{0}(\sqrt{3}maH^{2}/k^{3})\sin(k(\tau-\tau_{k})/\sqrt{3})$, which decay as $a^{-2}$ and $a^{-1}$, respectively, and can be ignored.

On the other hand, the homogeneous solution derived below exhibits constant-amplitude oscillations and will dominate the particular integral. It is convenient to choose the reference epoch $a_{*}=a_{m}$ and $H_{*}=H_{m}$, where $a_{m}$ and $H_{m}$ refer to the expansion factor and Hubble parameter at $2H=ma$, or $H_{m}\equiv ma_{m}/2$. Multiplying Eq. (16) by $a$, we can change the time variable to a dimensionless $\eta=\ln(a/a_{m})\equiv\ln(H_{m}\tau)$. Equation (16) can be cast into

[TABLE]

This equation admits the homogeneous solution $\hat{\psi}_{h}=g\exp[i(k^{2}/4H_{m}^{2})(\eta-\eta_{k})]$ with $g$ being a complex constant and $\eta_{k}=\ln(a_{k}/a_{m})$ denoting the duration between the onset of mass oscillation and when mode $k$ enters horizon. The homogeneous solution is oscillating with a constant amplitude and dominates the particular integral for $a>>a_{k}$. Thus the perturbed field $\hat{\psi}$ decouples from the gravity in this regime and becomes a free matter wave with a constant frequency $k^{2}/4H_{m}^{2}$ in the $\eta$ space.

In the long-wave limit, where $k^{2}/4H_{m}^{2}\ll 1$, we can estimate the complex constant $g$ by matching the solution to Phase (ia), which has a constant $\Re[\hat{\psi}]\sim O(1)\phi_{0}$ and a growing $\Im[\hat{\psi}]\sim O(ma/H)\phi_{0}(>>\Re[\hat{\psi}])$. We therefore expect that in the limit of small $\eta-\eta_{k}$ in the present phase, $\Re[\hat{\psi}]\ll\Im[\hat{\psi}]$. This is only possible when $\Re[\hat{\psi}]\propto\sin[(k^{2}/4H_{m}^{2})(\eta-\eta_{k})]$ and $\Im[\hat{\psi}]\propto\cos[(k^{2}/4H_{m}^{2})(\eta-\eta_{k})]$. That is, $g$ is pure imaginary. Moreover, we can also estimate that $|g|\sim O(H_{m}^{2}/k^{2})\phi_{0}$, obtained from the small argument expansion of $\sin[(k^{2}/4H_{m}^{2})(\eta-\eta_{k})]$ for $\Re[\hat{\psi}]$. This estimate of $|g|$ is also consistent with $\Im[\hat{\psi}]$ at the end of Phase (ia) when $k\to H$, since the factor $ma/H=maH/H^{2}\to H_{m}^{2}/2k^{2}$.

To fix the value of $|g|$, we use the CDM limit in Appendix (A) with $H/ma\to 0$. The CDM covariant energy perturbation is $\Delta_{CDM}\to(c_{1}-9(\eta-\eta_{k}))\phi_{0}$, when $c_{1}$ is a constant of order unity. On the other hand, the expression of the covariant energy perturbation $\Delta_{\psi}$ is the same as Eq. (18),

[TABLE]

Matching $\Delta_{CDM}$ and $\Delta_{\psi}$ in the limit $k^{2}/4H_{m}^{2}(\eta-\eta_{k})\to 0$ but $\eta-\eta_{k}\gg 1$, we find that $|g|=-18(H_{m}^{2}/k^{2})\phi_{0}$ from the sine term of $\Delta_{\psi}$. We stress that this result is valid only for long waves, and for short waves $g$ is generally complex.

Note that all modes, except for the very long-wavelength ones that have not yet entered horizon at the radiation-matter equality, must go through this final phase of the radiation era. Equation (16) is also valid in the matter era, and the solution characteristics deviates from the above description since the self-gravity of $\psi$DM becomes important. This equation in the matter era has been discussed previously and an analytical solution been obtained Woo, and Chiueh (2009).

(ii) Before mass oscillation

Before mass oscillation, the $\psi$DM energy perturbation is always much smaller than the radiation energy perturbation unless particle mass $m>10^{-28}$ eV, and hence passive evolution can always be a good approximation throughout this regime. In this regime, $\Psi$ is almost a constant, $\Psi^{{}^{\prime}}\sim-(m^{2}a^{2}/5H)\Psi$, and hence the equation of state of the zero-order field is like that of an inflaton with $P_{\psi}\sim-\epsilon_{\psi}$.

(iia) $ma,k\ll 2H$

The super-horizon metric perturbation is again almost a constant, i.e., $\phi\sim\phi_{0}(1-(1/30)(k^{2}/H^{2})+O(k^{4}/H^{4}))$. The perturbed field equation Eq. (2) can be approximated by:

[TABLE]

to the leading order. The three terms originally in Eq. (2) ignored in Eq. (24) are of order $O(k^{2}/H^{2})$ or $O(m^{2}a^{2}/H^{2})$ compared with other terms retained. Similar to Phase (ia), it turns out that the solution of Eq. (24) also yields a $\Delta_{\psi}$ that cancels itself to the leading order.

To see the cancellation, we shall find the solution for Eq. (24). The particular integral can be easily obtained by substituting $\psi^{{}^{\prime}}\propto(m^{2}a^{2}/H)\Psi\phi$ into Eq. (24) and it follows that $\psi_{p}\sim-(m^{2}a^{2}/10H^{2})\Psi\phi_{0}$. The homogeneous solution decays as $\psi^{\prime}\sim a^{-2}$ or $\psi^{\prime}=0$, and can be ignored. The energy perturbation for the particular integral becomes $\delta_{\psi}\sim(-3/25)(ma/H)^{2}\phi_{0}$. On the other hand from Eq. (10), $\Psi^{\prime}\psi\sim(1/25H)(ma/H)^{2}(m^{2}a^{2}\Psi^{2}/2)\phi_{0}$, and thus the gauge covariant energy perturbation $\Delta_{\psi}$ exactly cancels to the leading order of $a^{4}$.

All corrections of to $\psi_{p}$, and therefore $\Delta_{\psi}$, are integer powers of $(k/H)^{2}$ as only $k^{2}$ appears in the original equation, Eq.(2). Thus, $\Delta_{\psi}\sim O((k/H)^{2}\delta_{\psi})\sim O((k/H)^{2}(ma/H)^{2}\phi_{0})\sim O((ma/H)^{2})\Delta_{\gamma}\propto a^{6}$. Here we have again employed the Poisson equation to relate $\Delta_{\psi}$ to $\Delta_{\gamma}$.

The $a^{6}$ rapid growth in this phase is drastically different from the CDM perturbation $(\propto a^{2})$, and it is the main cause of the sharp cutoff in the matter power spectrum near a critical $k$ to be discussed in the next section. Moreover, this $a^{6}$ rapid growth occurs regardless of whether $k>ma$ or $ma>k$ and modes of all $k$ must go through this initial phase.

(iib) $ma\ll 2H\ll k$

In this regime, the background field is still like an inflaton, where the background energy density is dominated by the field potential energy. Same as Phase (ib), the sub-horizon metric perturbation undergoes a damped oscillation, i.e., $\phi\sim\phi_{0}(H/k)^{2}\cos(k(\tau-\tau_{k})/\sqrt{3})\propto a^{-2}$. This source drives the perturbed field to also oscillate with the same frequency. However, unlike Phase (ia), the driving frequency $k/\sqrt{3}$ is different from the natural frequency $k$ of the $\psi$DM perturbation, resonance is impossible and hence the treatment is different from Phase (ia).

The source now is dominated by $\Psi^{\prime}\phi^{\prime}\sim(k/5\sqrt{3}H)(m^{2}a^{2})\Psi\phi_{0}\sin(k(\tau-\tau_{k})/\sqrt{3})$ and is $\propto a$. Hence the particular integral to the approximate equation

[TABLE]

is $\psi_{p}\sim 6\Psi^{\prime}\phi^{\prime}/k^{2}\sim O((m^{2}a^{2}/kH)\phi\Psi)\propto a\sin(k(\tau-\tau_{k})/\sqrt{3})$. The energy perturbation $\delta_{\psi}\sim 2(\psi^{\prime}\Psi^{\prime}-(\Psi^{\prime})^{2}\phi)/m^{2}a^{2}\Psi^{2}\sim O((ma/H)^{2}\phi)\sim O((ma/k)^{2}\Delta_{\gamma})\propto a^{2}$. The gauge-covariant energy perturbation $\Delta_{\psi}$ differs from $\delta_{\psi}$ by a high-order term $O(H/k)$ smaller, and therefore $\Delta_{\psi}\sim\delta_{\psi}$, growing as $a^{2}$ in the oscillation amplitude.

Note that only for sufficiently large $k$ can a mode enter Phase (iib), and after that it skips Phase (ia) to go directly to Phase (ib). For a mode of low $k$, it initially goes through Phase (iia) and skips this phase to enter Phase (ia) and then (ib).

IV Evolution of Perturbations in Full Treatment

Full treatment of perturbations in the radiation-dominated era must take into account several effects beyond the passive evolution, namely, the decoupling of neutrino after entering horizon, the relative drag between baryons and photons, the non-negligible baryon component in the photon fluid energy density, and the self-gravity of matter. The last three become significant only near the radiation-matter equality.

(i) Neutrino Decoupling

After entering horizon $k\geq 2H$, photon fluid perturbations oscillate with a constant amplitude indefinitely until the effect of baryon-photon drag sets in. On the other hand, neutrino perturbations oscillate along with photon fluid perturbations in the first half cycle of oscillation, and after crossing the first null neutrino perturbations are abruptly damped out due to free streaming Ma, and Bertschinger (1995). We therefore set neutrino density perturbations to zero immediately after its first oscillation null, as an approximation to solve for perturbations of the rest of surviving species. Thus, we let

[TABLE]

where $\Delta_{\gamma}\equiv\Delta_{ph}+\Delta_{\nu}$, and $\Delta_{\nu}$ and $\Delta_{ph}$ are neutrino and photon covariant perturbed energy densities, respectively; $\tau_{1}$ is the conformal time of the first oscillation null of $\Delta_{\gamma}$.

(ii) Photon Fluid Equation of State

In between the end of radiation-dominant era, $a=a_{eq}$, and the onset of photon-electron decoupling, $a\approx 3a_{eq}$, baryons can be non-negligible in modifying the equation of state parameter for the photon fluid, i.e.,

[TABLE]

According to the current values $\Omega_{b}\approx 0.05$, $\Omega_{m}=\Omega_{dm}+\Omega_{b}\approx 0.316$ and $(\Omega_{ph}+\Omega_{\nu})/\Omega_{ph}\approx 1.7$ Planck Collaboration et al. (2016), the value of $w_{ph}$ approximately equals $20/77$ and $20/111$ when evaluated at the radiation-matter equality and at immediately before photon-electron decoupling, respectively.

In practice, we should not consider a single coupled photon$+$baryon fluid, but should keep track of photon and baryon perturbations separately for the reason given below. Their respective density equations and coupled momentum equations through Thomson scattering are

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Here $n_{e}$ is the electron number density and $\sigma_{T}$ is the Thomson scattering cross section.

(iii) Baryon-Photon Drag

Thomson scattering that couples photons and baryons can be characterized by a scattering mean-free-path $l_{T}\equiv(n_{e}\sigma_{T}a)^{-1}$. After the photon perturbation enters horizon, the wave number $k$ becomes a relevant scale. Here, we can have a length ratio, $kl_{T}\propto a^{2}$. When $kl_{T}\ll 1$, the photon perturbation is efficiently coupled to baryons, and photons are well-described by a fluid. However, there is another length ratio $(k/2H)(kl_{T})\propto a^{3}$ that is larger than $kl_{T}$ and becomes greater than unity earlier than $kl_{T}$ does after perturbations entering horizon. The new length scale ratio characterizes when the phase lag between the oscillation of baryon and the oscillation of photon becomes so severe that both photon and baryon perturbations are rapidly damped when $(k/H)(kl_{T})>1$. Physically, photons undergo random walks due to Thomson scattering, and in a Hubble time $H^{-1}$ the random walk distance is $(l_{T}/H)^{1/2}$. When this distance is larger than the wavelength, the drag damping can occur Hu, and Sugiyama (1996). In the wave number range of interest to the present work, $k\sim O(k_{c})$, we find $(k/2H)(kl_{T})\sim O(1)$ near the radiation-matter equality and the drag damping just begins to take effect. At this moment, we still have $kl_{T}\ll 1$ and the photon fluid is a good description and it justifies the photon momentum equation Eq. (29).

(iv) Matter Self-Gravity

In the full treatment, the source of Poisson equation contains perturbations of all species, to be contrasted with passive evolution where only radiation perturbations are the only source of Poisson equation. The matter gravity only becomes important near the radiation-matter equality, where the self-gravity prolongs the matter-wave oscillation cycle and slightly modifies the $\psi$DM power spectrum.

A 4-th order Runge-Kutta scheme is adopted for integrating the background equations (Eqs. (1) and Friedmann equations described in Appendix (B), i.e., Eqs. (B1) and (B2), with equations of states) and the perturbed equations (Eqs. (2), (28), (29), (30), (31) and perturbed Einstein equations described in Appendix (B), i.e., Eqs. (B3) and (B4), with the neutrino contribution given in Eqs. (B5) and (B6)). Plotted with dotted lines in Fig. (2) are the numerical solutions of the covariant energy perturbations $\Delta_{\psi}$ for full treatment as functions of the scaling factor $a/a_{eq}$ for particle mass $m=10^{-22}$eV. We have chosen modes with $k\gg k_{c}$ and $k\ll k_{c}$ to illustrate the difference. The vertical solid line marks the onset of mass oscillations ($a=a_{m}$). The three asymptotic phases discussed in the last section are clearly shown in Fig. (2) for the low-$k$ mode and the other three for higher-$k$ mode. Solutions of passive evolution are also plotted as solid lines in Fig. (2) for comparison. One notices that the full treatment and the passive evolution deviate from each other only near the radiation-matter equality, indicative of the slowing down of matter-wave oscillation due to the matter self-gravity. It is not surprising to find that treating neutrinos as a fluid for passive evolution does not differ from a more sophisticated approximation for neutrinos in the full treatment, as long as $k$ is smaller than $k_{c}$. This is due to that sub-horizon perturbations of these low-$k$ modes are hardly coupled to metric perturbations (Phase(ib)). Noticeable differences occur only for $k>k_{c}$ since these modes respond to photon oscillations driven by metric perturbations after entering horizon (Phase(iib)), and a decrement in the metric perturbation due to the vanishing neutrino contribution decreases the oscillation amplitude. These high-$k$ modes, however, have vanishingly small power at the radiation-matter equality and hence the differences are insignificant for practical purposes.

In Appendix (B), we show in detail how the photon fluid and the baryon fluid, coupled by the relative drag, evolve in the full treatment. To summarize, passive evolution of $\psi$DM captures most essential physics in the evolution of full treatment. Hence we shall continue to take advantage of the understanding gained from passive evolution for the following discussions.

V Critical Mode, Matter Power Spectrum and Sub-Horizon Dynamics

We saw in the previous discussions that high-$k$ modes go through Phases (iia), (iib) and (ib), and low-$k$ modes go through Phase (iia), (ia) and (ib). In between the two, there exists a critical $k_{c}$ mode crudely defined as $k_{c}=ma_{m}=2H_{m}=(2H_{m}ma_{m})^{1/2}=(2Hma)^{1/2}$, a redshift-independent wave number, that characterizes a transition. This particular mode enters horizon at the moment when the mass oscillation starts, and so the mode can be viewed as having a wavelength equal to the Compton wavelength when $\psi$DM is crystalized into a real particle. For modes of $k<k_{c}$, they have wavelengths greater than the Compton length ever since the particle becomes real; for modes of $k>k_{c}$, there is a finite period after the particle becomes real where their wavelengths are smaller than the Compton length. One expects that structures smaller than the Compton length can hardly exist and hence must be suppressed. This is indeed what we saw in the drastically different dynamics for modes of $k<k_{c}$ and $k>k_{c}$.

Plotted in Fig. (3) are the $\psi$DM transfer function in reference to CDM, $|\Delta_{\psi}|^{2}/|\Delta_{CDM}|^{2}$, evaluated at $a_{eq}$ as functions of $k$ with particle mass, $m=10^{-22}$eV. The solid line denotes the passive evolution and the dashed line the full treatment. Comparison of the two curves shows minor ($<10\%$) differences for $k<k_{c}/2$, demonstrating the passive evolution is a good approximation for the entire radiation-dominant era for these long waves. For $k>k_{c}/2$, our discussion about the mode suppression still holds and details are given below. In addition, Fig. (3) exhibits substantially different oscillation phases for the two (passive versus full treatment) evolutionary paths. This is caused by the self-gravity in full treatment that lowers the matter-wave oscillation frequency. But the self-gravity comes to play late in the evolution, and hence little affects our dynamical picture of Phase (ib).

In Fig. (3), we also provide the power spectrum evaluated at $2.5a_{eq}$, serving as the initial condition for any calculation starting in the matter-dominated era.

Below, we shall explain how the short-wave suppression occurs more quantitatively. All modes first follow the rapid $a^{6}$ growth of phase (iia). Modes with $k<k_{c}$ continue this growth until the mass oscillation starts at $a=a_{m}$, and then they follow the CDM growth before entering horizon. Modes with $k>k_{c}$ exit the rapid growth of Phase (iia) early on when becoming sub-horizon and enter the slower $a^{2}\cos(k(\tau-\tau_{k}))$ growth of Phase (iib) before mass oscillation, and after mass oscillation, they enter Phase (ib) directly without going through the CDM phase, Phase (ia). Hence, these high-$k$ modes never have a chance to resemble CDM perturbations. The amplitudes differ roughly by $(k_{c}/k)^{4}$ for modes with wave numbers few times above $k_{c}$ in comparison with the critical mode at the onset of mass oscillation. Thus, the sharp spectral transition has already been established when mass oscillation starts and real particles of finite mass $m$ are crystalized.

Additional suppression also occurs in Phase (ib), which mostly affects modes near $k_{c}$. This suppression has nothing to do with the Compton length effect mentioned earlier, but something to do with the Jeans length. While the CDM mode still grows mildly as $\ln(a/a_{k})$, the $\psi$DM mode in Phase (ib) oscillates with a constant amplitude, as $(2Hma/k^{2})\sin[(k^{2}/2Hma)\ln(a/a_{k})]$, leading to a $k^{-2}$ suppression for shorter waves. This Jeans suppression in the radiation-dominant regime is an extension of the more familiar counterpart in the matter-dominant regime.

Despite the formula given in the last paragraph for $\psi$DM modes was derived for $k\ll k_{c}$, it still roughly holds for $k\sim k_{c}$. Take the critical mode with particle mass $m=10^{-22}$ eV as an example. The mass oscillation begins around the photon temperature $T_{\gamma}\sim 0.5$ keV, occurring $1/500$ times smaller in $a$ than that at the radiation-matter equality with $T_{\gamma}\sim 1$ eV, and hence the oscillation amplitude of the critical mode is $(\ln(a_{eq}/a_{m}))^{-1}$ smaller than that of the CDM mode of the same $k$. We thus have an oscillation amplitude suppression factor of about $6.2$, pretty consistent with Fig.(3). In general, the oscillation amplitude suppression factor is approximately $(\ln(500k/k_{c}))^{-1}(k_{c}/k)^{2}$ around $k_{c}$. This formula is still valid for a different particle mass $m$. For that purpose one simply replaces $500$ by $a_{eq}/a_{m}$.

Accompanying the suppression of the high-$k$ power is the oscillation pattern for shorter waves in the power spectrum. There are two types of oscillations separately exhibited by modes of $k\gg k_{c}$ and $k\sim k_{c}$. For longer waves, the oscillation pattern arises primarily from the matter-wave oscillation and for shorter waves from the photon fluid oscillation followed by the matter-wave oscillation. Although the oscillation pattern in the power spectrum for shorter waves is complicated, due to the mixture of two oscillations of different frequencies, the one for longer waves is predictable. The phase of the latter oscillation is given by $(k/k_{c})^{2}\ln(a_{eq}/a_{k})$ and hence the oscillation frequency in the power spectrum is $[(k/k_{c})^{2}(\ln(a_{eq}/a_{m})+\ln(k_{c}/k))]^{2}$. The transition between the two oscillation patterns for longer waves and for shorter waves is found roughly at $k=2.5k_{c}$ for $m\sim 10^{-22}$eV. (The power spectrum shown in Fig. (3) is hence not extended to cover such a short-wave regime.)

One can conveniently define $k_{1/2}$ as the wavenumber for which the modal power is suppressed by $1/2$ relative to CDM. Taking into account the phase and amplitude of the matter-wave oscillation, we find that $k_{1/2}/k_{c}\sim(\alpha/\ln(m/m_{28}))^{1/2}$ for any $m\gg m_{28}$ and $\alpha\approx 1.3$ for passive evolution and $\alpha\approx 1.7$ for active evolution, where $m_{28}\equiv 10^{-28}$eV. The mass dependence arises from that the onset of mass oscillations for a more massive particle occurs earlier $(a_{m}\propto m^{-1/2})$, and therefore it has a more ample time to execute matter-wave oscillations before the radiation-matter equality.

Interestingly the Jeans wave number $k_{J}$ in the matter dominated regime, for which modes with $k<k_{J}$ grow like CDM and modes with $k>k_{J}$ oscillate with constant amplitudes, is also $\sqrt{2Hma}$. Moreover in the matter-dominated regime, $H\propto a^{-1/2}$ and hence the Jeans length $k_{J}\propto a^{1/4}$, mildly evolving toward shorter wavelengths. This explains why $k_{J}$ is always larger than, but on the same order of, $k_{c}$. The oscillation pattern in the power spectrum near $k_{c}$ is therefore frozen in the matter-dominant era and grows self-similarly as $a^{2}$.

It is relevant to also compute the phase $S$ of $\hat{\psi}$, which, along with the power spectrum, can determine $\psi$ used for simulations in the matter-dominant regime as the initial condition. The power spectrum $\langle(\Re[\hat{\psi}])^{2}\rangle=\langle|\hat{\psi}|^{2}\rangle\cos^{2}(S)$ only provides a partial condition for this purpose. When $\Re[\hat{\psi}]/|\hat{\psi}|=\cos(S)$ and $\Im[\hat{\psi}]/|\hat{\psi}|=\sin(S)$ are available, one can construct $\exp(iS)$ for each $k$ mode. The amplitude $|\hat{\psi}|$, being a half-Gaussian random number, can be then fixed by constructing an exponential random number $|\hat{\psi}|^{2}\cos(S)^{2}$ and letting its average equal to the power spectrum. Figure (4) shows the ratio $\Re[\hat{\psi}]/|\hat{\psi}|$ as functions of $k$ at the matter-radiation equality ($a=a_{eq}$) and at $a=2.5a_{eq}$ using the full treatment. The first thing brings to notice is $\Im[\hat{\psi}]\gg\Re[\hat{\psi}]$ for long waves. This feature is apparent from the discussions of Phase (ib). Secondly, $\Im[\hat{\psi}]$ grows substantially faster than $\Re[\hat{\psi}]$ in between $a=a_{eq}$ and $a=2.5a_{eq}$ due to the acceleration of self-gravity, and the ratio $\Re[\hat{\psi}]/\Im[\hat{\psi}]\to 0$ from above for most $k<k_{c}$ modes. These modes are purely growing modes. For $k>k_{c}$, the wave function oscillates with a time-varying frequency, and in the range $k_{c}<k<1.3k_{c}$ plotted here, the oscillation frequency is nearly zero, indicative of that the wavenumber is close to the Jeans wavenumber $k_{J}$.

VI Adiabatic Perturbations for Superhorizon Modes

In the multi-fluid model, the condition for adiabatic perturbation is Eq. (19). Here $\langle w\rangle$ is the short-time averaged equation of state parameter $\langle w\rangle=\langle P/\epsilon\rangle$.

Before the mass oscillation in Phase (iia), there is no fast mass oscillation and the equation of state for $\psi$DM is

[TABLE]

and in this phase $\psi$DM behaves similar to the dark energy. In Appendix (A) we see $\delta_{\gamma}\sim-2\phi$ in this phase, approximately a constant. Hence the adiabatic condition yields $\delta_{\psi}=-(3/25)(m^{2}a^{2}/H^{2})\phi$, the same as we found by solving the approximate dynamical equation.

After mass oscillation in Phase (ia), we have used the adiabatic condition to pin down the value of the super-horizon energy perturbation, which involves an undetermined constant. In fact, one may question whether Eq. (19) in Phase (ia) really needs to take a short-time average. This question can be answered by examining Eq. (17) and the fact that $1+w_{\psi}\propto(\Psi^{\prime})^{2}\propto\sin^{2}(mt+\alpha)$. Since $\delta_{\gamma}/(1+w_{\gamma})$ has no mass oscillation but both $\delta_{\psi}$ and $w_{\psi}$ do, one should examine whether the respective oscillations of $\delta_{\psi}$ and $1+w_{\psi}$ cancel. We had before $\Re[\hat{\psi}]\sim-(3/4)\phi_{0}$, and $\Im[\hat{\psi}]_{p}=-(ma/2H)\phi_{0}$. Upon substituting them into Eq. (17), we find $\delta_{\psi}$ is also proportional to $\sin^{2}(mt+\alpha)$, and $\delta_{\psi}/(1+w_{\psi})$ indeed has no oscillation. Therefore the adiabatic condition holds instantaneously without the need to take a short-time averages over $\delta_{\psi}$ and $w_{\psi}$, even though the perturbation has fast oscillations.

Next, one may wonder why the long-wavelength adiabatic perturbations of $\psi$DM and CDM should track each other so well after the mass oscillation. A rough answer is that since they are both adiabatic and when the equations of state are the same, the two must be identical when super-horizon. As we have seen, the average equation of state parameter $<w_{\psi}>=0$ in Phase (ia), the same as that of CDM, which justifies the above statement. Further considerations show that once the mode enters Phase (ib), we also have the same fluid description for $\psi$DM as that for CDM in the long wave limit, since $\delta P_{\psi}=O((k/ma)^{2}\delta\epsilon_{\psi})\to 0$ when $k\to 0$. Therefore, for sufficiently long waves, $\psi$DM and CDM become almost identical after mass oscillation, despite the two are very different before mass oscillation.

VII Axion Model

Axion has a field potential $V=(fm)^{2}(1-\cos(\Psi/f))$, where there appears a new energy scale, the axion decay constant $f$, believed to be close to the GUT scale or the Planck scale Chiueh (2014); Hui et al. (2016). In the limit $\Psi/f\to 0$, the axion potential is reduced to the free-particle harmonic potential discussed previously. Normally without fine tuning, it is expected that the initial axion angle $\theta(\equiv\Psi/f)$ is of order unity at the very early epoch. (Here we have taken a simple axion model for which the field potential has been established long before the mass oscillation.)

Unlike the free particle model, the zero-order field obeys a nonlinear equation:

[TABLE]

Similar to Eq. (1), $\theta$ can oscillate, but initially executing anharmonic oscillation. Subject to the Hubble friction, $\theta$ declines in amplitude after a couple of oscillations where $\sin(\theta)\to\theta$, and Eq. (33) becomes Eq. (1). Thus, the axion model is a straightforward extension of the free-particle model, and can address the initial phase-lock problem of $\psi$DM explained in Sec. (I).

The perturbed field now obeys

[TABLE]

which recovers Eq. (2) when $\theta<<1$. The only difference of the axion model from the free-particle models occurs near the $ma=2H$ transition, and the one additional free parameter for the axion model is the initial $\theta$, which we denote as $\theta_{0}$.

In contrast to the free-particle model, the axion model has no analytical solutions for the zero-order field, and hence the analytical solution to the first-order field is not attainable even for passive evolution. We replace Eqs. (1) and (2) by Eq. (33) and Eq. (34) and numerically integrate the coupled equations of the axion model as we did for the full treatment of the free-particle model. The results are plotted in Fig. (5) for $\theta_{0}=\pi/16,\pi/2,35\pi/36$, where the $\theta_{0}=\pi/16$ case is essentially the free-particle model.

Perturbations of the axion model also roughly have the four asymptotic solutions corresponding to the four asymptotic phases of the free-particle model. This is actually not surprising since the major difference between the two models is just during the time near the onset of mass oscillation; for a difference in such a short time, the asymptotic phases are little affected. Generally speaking, the initial anharmonic oscillation of the axion model has longer periods than the harmonic oscillation of the free-particle model and would lead to some phase delay of mass oscillation. But solutions at $a\gg a_{m}$ turn out not to be sensitive to such a moderate phase shift. This is clearly shown in the power spectrum of Fig. (6) by a comparison of solutions of the free-particle model and the axion model of a moderate initial angle $\theta_{0}=\pi/2$. The difference is at most $10\%$ for $k\leq k_{1/2}$. Only in the extreme case when $\theta_{0}\to\pi$, where the onset of mass oscillation is significantly delayed causing the axion’s first cycle of matter-wave oscillation discussed in Phase (ib) to take a substantially longer period than that of the free particle, can it make an appreciable difference. The longer matter-wave oscillation cycle changes the phase evaluated at $a=a_{eq}$ and can yield a higher $k_{1/2}$, as a comparison of $k=k_{c}$ cases for three different $\theta_{0}$’s hints in Fig. (5). This feature makes the axion of extreme $\theta_{0}$ appear as if the particle had a higher mass assuming it is free particle that follows $k_{1/2}\sim O(k_{c})\propto m^{1/2}$. Surprisingly, the axion model of extreme initial angle can also produce a slightly higher power than the CDM model for the extreme $\theta_{0}\to\pi$ case, which has never been so for the free-particle model, and the amplitude excess occurs near $k=k_{c}/2$ for $\theta_{0}=35\pi/36$. For the fiducial particle mass $10^{-22}$ eV, the power excess is in the wavelength range $>1$Mpc, which encloses about $10^{10}-10^{11}M_{\odot}$ of dark matter, a typical mass of first galaxies, and may have an impact on the first-galaxy formation.

How do we understanding the two new features, higher $k_{1/2}$ and amplitude excess near $k_{c}/2$ of the extreme axion model? Compared with the free-particle model, the axion perturbation $\delta\theta$ has much smaller sources initially to drive the growth when $\theta$ is close to $\pi$, as its potential gradient is almost zero. This is why $\Delta_{axion}\ll\Delta_{\psi}$ in Phase (iia). In the beginning, this small source is balanced by both ”inertia” and ”friction”, similar to the free-particle case, and $\Delta_{axion}\propto a^{6}$. Shortly after, $\theta$ is rolling down the potential hill, the field force much exceeds the friction, and the perturbation grows more rapidly than $a^{6}$. This is almost an exponential growth, but it only lasts briefly for a fraction of the first mass oscillation cycle. In the first few mass oscillation cycles after the rapid growth, the axion perturbation grows roughly as $a^{3}$, as contrasted to the $a^{2}$ growth of the free-particle perturbation, and does the final catch-up with the free particle perturbation. The subsequent evolution of the axion perturbation after entering horizon pursues a slightly different evolutionary path from the free-particle perturbation, albeit both types of perturbations at this point satisfy almost the same Eq. (22) except for a small nonlinear correction in the axion case.

To assess this problem properly, we must keep the nonlinearity, albeit small, of the background $\theta$. The nonlinear mass oscillation causes nonlinear frequency shifts in both $\theta$ and $\delta\theta$. But two frequency shifts are slightly different, leading to additional detuning between the response and the driving source, similar to the role of $k$. (See the first paragraph of Phase (ia).) Since the matter-wave oscillation frequency is the detune frequency, the nonlinear effect changes the matter-wave oscillation frequency. The nonlinearity decreases with time, and hence the earlier the mode enters horizon, the more prominent this effect is. It explains why the very-long wavelength extreme axion perturbation, entering horizon very late, has no effect and is almost identical to the free-particle perturbation.

Moreover, the detune frequency yielded by the nonlinear oscillation of $\theta_{0}$ has two components, a quasi-static component and a rapidly oscillating component. It turns out the rapidly oscillating component drives the perturbed field into a parametrically unstable state and produces the power excess.

The above mechanisms qualitatively explain the most important two features of the extreme axion model. Given the implications that the extreme axion model can yield a higher particle mass effect in the spectral cutoff and can also produce a higher power than CDM in the wavelength range that yields first galaxies, this subject warrants a quantitative analysis, and we report it in a separate paper Zhang and Chiueh (2017).

VIII Conclusions

In this paper, why the matter power spectrum at the end of radiation era exhibits a CDM-like spectrum for long waves and a highly suppressed spectrum for short waves and why the transition of these two types of spectra is sharp are fully explored. Physically the sharp power suppression is due to the smearing of fluctuations below the Compton length. Dynamically it is due to short waves prematurely entering horizon before mass oscillations, so that quantum pressure becomes so effective as to slow down the otherwise very rapid growth which long waves entertain.

To summarize, the following issues are addressed and results obtained:

(1) We elucidate four asymptotic phases with passive evolution— before mass oscillations, after mass oscillations, super-horizon wavelengths and sub-horizon wavelengths. The transition occurs at the intersection of the four asymptotic phases, i.e., at the onset of mass oscillations for the wavelength just beginning to enter the horizon, i.e., $k=2H=ma$. It defines a critical wavenumber $k_{c}$ which coincides the transition of the two distinct behaviors in the $\psi$DM spectrum. This $k_{c}$ is the predecessor in the radiation-dominant era of the Jeans wavenumber $k_{J}$ in the matter-dominant era.

(2) We demonstrate passive evolution of $\psi$DM to be a good approximation of the full treatment, and we also address the details of the full treatment, including collisionless neutrinos and Thomson scattering coupled photons and baryons, through which numerical solutions are obtained.

(3) As a bonus of our analyses, we show that the adiabatic condition for super-horizon perturbations holds instantaneously even when the field potential produces fast mass oscillations.

(4) The phase and amplitude of the complex $\hat{\psi}$ are important quantities for any $\psi$DM simulation to start the initial condition faithfully in the matter dominant era. For this reason, we extend our numerical solutions into the matter-dominant era and compute the phase of $\hat{\psi}$ and the power spectrum beyond the radiation-matter equality.

(5) We also consider perturbations of the axion model, a plausible extension of the free-particle model of $\psi$DM. It is found that the evolution of axion perturbations is generally almost identical to that of free-particle perturbations, except for the extreme case where the initial value of the axion angle is near the field potential top, yielding a substantially longer frequency of matter wave oscillation and shifting the transition wavenumber to a higher value. The upward shift in the transition wavenumber is equivalent to having a higher particle mass for the free particle model. This extreme axion model also produces a higher power than CDM in a wavelength range about half of the critical wavenumber. For particle mass $\sim 10^{-22}$eV, the enclosed mass of these wavelengths is the typical mass of first galaxies, which may have impact to first galaxies at the cosmological redshift $z\sim 10$, discovered by recent observations Bouwens et al. (2011); Oesch et al. (2012); Zheng et al. (2012); Coe et al. (2013); Oesch et al. (2014); Zitrin et al. (2014); Bouwens et al. (2015); Bernard et al. (2016).

This work is supported in part by MOST of Taiwan under the grant MOST 103-2112-M-002-020-MY3.

Appendix A Passive Evolution

In this appendix, governing equations and analytic solutions for passive evolution are presented.

The action can be written as follows,

[TABLE]

where $g$ is the determinant of the metric $g_{\mu\nu}$, $R$ the scalar curvature, $L_{M}$ the Lagrangian density of matter and $G$ the gravitational constant. Here we choose the speed of light and the Planck constant $\hbar$ to be $1$.

The Lagrangian density $L_{M}$ consists of two components, the radiation fluid and the dark matter. The radiation can be described as a perfect fluid with the equation of state $P=\epsilon/3$, where $P$ is the pressure and $\epsilon$ is the energy density. On the other hand, two different dark matter models, the cold dark matter (CDM) model and the $\psi$ dark matter ($\psi$DM), will be introduced. CDM is a pressureless perfect fluid, and $\psi$DM is a spin-zero real scalar field with a Lagrangian density

[TABLE]

where $V$ is the scalar field potential. The corresponding stress-energy tensor is

[TABLE]

The metric is the flat Friedmann-Lemaitre-Robertson-Walker (FLRW) metric with the scalar metric perturbation, i.e.,

[TABLE]

where $\tau$ is the conformal time, $x^{i}$ with $i=1,2,3$ is the comoving coordinate, $\phi$ the metric perturbation, $a$ the scale factor, $\delta_{ij}$ the Kronecher delta function, and the Newtonian gauge has been chosen. Similarly, the energy-momentum stress tensor is also written as a uniform background with a small perturbation. Thus $\psi$DM field can be expressed as $\Psi(\tau)+\psi(x^{i},\tau)$ for calculating the corresponding stress-energy tensor.

Gathering the above and substituting into Eq. (A1) and applying Euler-Lagrange equation, one obtains the Einstein equations and the conservation of energy and momentum.

The zeroth order equations are just Friedmann equations,

[TABLE]

[TABLE]

with the following conservation laws,

[TABLE]

where $()^{{}^{\prime}}$ is the derivative with respect to $\tau$, $H\equiv a^{{}^{\prime}}/a$ the Hubble parameter, and the subscript denotes different components. Here we use $\gamma$ to denote the radiation fluid and $D$ the dark matter.

For the CDM model, the pressure $P_{D}$ is zero, and Eq. (A7) becomes,

[TABLE]

where we have changed the subscript ”$D$” to ”$C$” to label CDM. For the $\psi$DM model, we have the following relations

[TABLE]

[TABLE]

where the subscript ”$\psi$” is used to replace ”$D$”. These relations can be obtained by Eq. (A3). Substituting Eqs. (A9) and (A10) into Eq. (A7), it follows that

[TABLE]

The first order perturbed equations consist of the perturbed Einstein equations

[TABLE]

[TABLE]

the continuity equation

[TABLE]

and the Euler equation

[TABLE]

Here $\delta\epsilon_{\alpha}$, $\delta P_{\alpha}$ and $\theta_{\alpha}$ are the energy density perturbation, the pressure perturbation and the velocity potential of the $\alpha$ species, respectively. The quantity $\delta_{\alpha}$ is the fractional overdensity of the $\alpha$ species defined as $\delta_{\alpha}\equiv\delta\epsilon_{\alpha}/\epsilon_{\alpha}$. The spatial Fourier transformation has been applied in the comoving coordinate on Eqs. (A12), (A13), (A14) and (A15).

For the CDM model, $P_{D}=\delta P_{D}=0$, Eqs. (A14) and (A15) become

[TABLE]

[TABLE]

On the other hand for $\psi$DM, from Eq. (A3), the quantities $\delta\epsilon_{D}$, $\delta P_{D}$ and $\theta_{D}$ can be expressed as

[TABLE]

[TABLE]

[TABLE]

Substituting Eqs. (A18), (A19) and (A20) into Eq. (A14), it follows that,

[TABLE]

In addition, the Euler equation (Eq. (A15)) is reduced to Eq. (A11), as shown in Sec. (II).

In the radiation-dominant era, the dark matter energy density is much smaller than the radiation energy density, so are the perturbed quantities. Under this condition, dark matter can be treated as passive particles (fields) tracing the geometry of universe that is determined by radiation, and therefore the dark matter contribution in Friedmann equations (Eq. (A5) and (A6)) and perturbed Einstein equations (Eqs. (A12) and (A13)) can be ignored.

The form of scalar field potential $V(\psi)$ should be specified. For a non-self-interacting free field, we have $V(\psi)=m^{2}\psi^{2}/2$, where $m$ is the mass of the dark matter particle.

First, we should solve the zeroth order equations. From Eq. (A7) for the photon fluid, the photon energy density $\epsilon_{\gamma}\propto a^{-4}$ due to $3P_{\gamma}=\epsilon_{\gamma}$. Substituting the above relation into Friedmann’s equations, we find $a\propto\tau$ and $H\propto a^{-1}$.

When the CDM model is considered, it is straightforward to find the energy density $\epsilon_{C}\propto a^{-3}$. For the $\psi$DM model, substituting above scalar potential $V(\psi)$ into Eq. (A11), Eq. (1) follows. It is natural to define a dimensionless quantity $u_{m}=ma/2H$ which normalizes the time to the onset of mass oscillation. With above relations, Eq. (1) becomes

[TABLE]

Equation (A22) has following solution,

[TABLE]

where $C^{\Psi}_{r}$ and $C^{\Psi}_{s}$ are constants, and $\Psi_{r}$ and $\Psi_{s}$ are

[TABLE]

[TABLE]

where $\Gamma$ is the gamma function, $J$ and $Y$ are Bessel function of the first and second kind, respectively, which are orthogonal to each other. It is noted that $\Psi_{s}$ has a singularity as $u_{m}\rightarrow 0$ while $\Psi_{r}$ remains regular, and hence only $\Psi_{r}$ is the valid zero-order field.

When $u_{m}\gg 1$, the energy density $\propto u_{m}^{-3/2}\propto a^{-3}$. This means the energy density of $\psi$DM follows CDM after the mass oscillation. We further define $\Psi_{0}\equiv C^{\Psi}_{r}\Psi_{r}(0)=C^{\Psi}_{r}$. In the following, $C^{\Psi}_{r}$ will be replaced by $\Psi_{0}$.

Next, we turn to the first order perturbed equations. The metric perturbation $\phi$ should be obtained first. To solve $\phi$, it is natural to define a dimensionless quantity $u_{k}\equiv k/\sqrt{3}H$. This quantity approximately measures the epoch of the horizon entry. Using the chain rule, Eq. (A13) becomes

[TABLE]

where $g_{\gamma}\equiv 4H\theta_{\gamma}$. On the other hand, from Eqs. (A5), (A6), (A12), (A13), (A14) and (A15) with the radiation fluid’s equation of state $P_{\gamma}=\epsilon_{\gamma}/3$ and $\delta P_{\gamma}=\delta\epsilon_{\gamma}/3$, $g_{\gamma}$ obeys the following equation

[TABLE]

Combining Eq. (A26) and (A27), one finds $\phi$ satisfies the following second order equation

[TABLE]

Equation (A28) has a solution:

[TABLE]

where $C^{\phi}_{r}$ and $C^{\phi}_{s}$ are constants, and

[TABLE]

[TABLE]

From Eqs. (A14) and (A15), the radiation energy perturbation $\delta_{\gamma}$ becomes,

[TABLE]

with

[TABLE]

[TABLE]

Since $\phi_{s}$ and $\delta_{\gamma}^{s}$ have a singularity at $u_{k}=0$, only the regular modes, $\phi_{r}$ and $\delta_{\gamma}^{r}$, are retained. Similar to the $\psi$DM background field case, we can define $\phi_{0}\equiv C^{\phi}_{r}\phi_{r}(0)=C^{\phi}_{r}$. In the following, $C^{\phi}_{r}$ will be replaced by $\phi_{0}$.

It is noted that $\phi(u_{k})=\phi_{0}(1+O(u_{k}^{2}))$, dominated by a constant $\phi_{0}$. On the other hand, from Eq. (12) ignoring the dark matter contribution, the gauge covariant photon energy density perturbation grows as $\Delta_{\gamma}\propto u_{k}^{2}\propto a^{2}$ before entering horizon.

Moreover, from Eq. (A33), $\delta_{\gamma}(u_{k})=\phi_{0}(-2+O(u_{k}^{2}))$, when $u_{k}\ll 1$. This relation is used in the discussion of phase (ia) in Sec. (III).

The perturbed dark matter equations are solvable once the background solution and the metric perturbation are given. For CDM model, Eqs. (A16) and (A17) become

[TABLE]

[TABLE]

where $g_{C}\equiv 3H\theta_{C}$ and the metric perturbation $\phi=\phi_{0}\phi_{r}(u_{k})$. It is straightforward to solve Eq. (A36) for the solution of $g_{C}$, as

[TABLE]

Here the homogeneous solution of $g_{C}$ has been dropped since it has the singularity at $u_{k}=0$. Substituting Eq. (A37) into Eq. (A35) and integrating on both sides, one obtains

[TABLE]

where $Ci(u_{k})$ is the cosine integral function and $\gamma$ is the Euler-Mascheroni constant. The last term, $1/2$, is fixed by the adiabatic condition (Eq. (19)). Finally, the gauge covariant energy perturbation $\Delta_{C}$ is obtained to be

[TABLE]

For $\psi$DM model, Eq. (A21) becomes Eq. (2)for free particle. Using $u_{m}$ as the argument, Eq. (2) further becomes

[TABLE]

where $A_{k}=k^{2}/(2maH)$, a dimensionless constant, and $S$ is the driving source. We note that $u_{k}^{2}=4A_{k}u_{m}/3$. Equation (A40) has following solution,

[TABLE]

Here $C^{\psi}_{r}$ and $C^{\psi}_{s}$ are constants, $\psi_{r}$ and $\psi_{s}$ are the homogeneous solutions of Eq. (A40) with the following forms,

[TABLE]

[TABLE]

where $M$ is the Kummer’s function, and

[TABLE]

the particular integral of Eq.(A40) with initial conditions $\psi_{p}(0)=d\psi_{p}/du_{m}(0)=0$. Here $W(\psi_{r},\psi_{s})(u_{m})=\psi_{r}d\psi_{s}/du_{m}-\psi_{s}d\psi_{r}/du_{m}=-u_{m}^{-3/2}/2$, the Wronskian of $\psi_{r}$ and $\psi_{s}$.

Again, the constants, $C^{\psi}_{r}$ and $C^{\psi}_{s}$, can be determined by the adiabatic condition for the super-horizon mode. Particularly, these two constants must be zero if only regular mode of metric perturbation is considered. Therefore $\psi=\psi_{p}$, representing the adiabatic perturbation.

Once the perturbed field’s solution is obtained, the gauge covariant energy density perturbation follows from Eq. (9), and becomes

[TABLE]

In the following, the asymptotic behavior of the gauge covariant dark matter energy density perturbation is to be derived. The CDM model is considered first. From Eq. (A39), it is straightforward to find

[TABLE]

For $\psi$DM model, it is more complicated and has four different phases as discussed in Sec. (III), which are to be verified from Eqs. (A41) to (A44). In the following presentation, we will organize the four phases in the chronological order. The phases before mass oscillation (phase (ii) in Sec. (III)) are dealt with first and then we will show the results of the phases after mass oscillation (phase (i)).

(A) Before mass oscillation and super-horizon mode (phase (iia))

This phase can be dealt with using Eq. (A44) directly. The background field $\Psi$ has the following asymptotic form $\Psi=\Psi_{0}[1-{u_{m}^{2}/5}+O(u_{m}^{4})]$, and the metric perturbation $\phi$ can be expressed $\phi=\phi_{0}[1-u_{k}^{2}/10+O(u_{k}^{4})]=\phi_{0}[1-2A_{k}u_{m}/15+O(u_{m}^{2})]$. To find the asymptotic form of the particular solution $\psi_{p}$, we first have,

[TABLE]

[TABLE]

Combining with above relations, it follows that $\psi_{p}$ becomes

[TABLE]

Finally, the denominator of Eq. (A45) is

[TABLE]

and the numerator

[TABLE]

where the leading order terms cancel. We thus have

[TABLE]

The gauge covariant $\psi$DM energy density perturbation $\propto u_{k}^{2}u_{m}^{2}\propto a^{6}$ when both $u_{k}$ and $u_{m}\ll 1$.

(B) Before mass oscillation and sub-horizon mode (phase (iib))

This phase is equivalent to $u_{m}\ll 1$ and $u_{k}\gg 1$, which imply $A_{k}\gg 1$. Unlike the previous case, $u_{k}$ is the natural variable here instead of $u_{m}$. The variables $\psi_{r}$, $\psi_{s}$ and the driving source $S$ will be expressed as functions of $u_{k}$.

From the definition of Kummer’s function $M$, the homogeneous solutions $\psi_{r}$ and $\psi_{s}$ are

[TABLE]

[TABLE]

On the other hand, the background field $\Psi$ has the asymptotic form $\Psi=\Psi_{0}[1-9A_{k}^{-2}u_{k}^{4}/80+O(A_{k}^{-4})]$. Hence, the driving source $S$ becomes

[TABLE]

Substituting Eqs. (A53), (A54) and (A55) into Eq.(A44) and changing the variable from $u_{m}$ to $u_{k}$, the particular solution $\psi_{p}$ becomes

[TABLE]

Therefore, the gauge covariant energy density perturbation $\Delta_{\psi}$ is

[TABLE]

The gauge covariant energy density perturbation grows as $u_{k}^{2}\propto a^{2}$ modulating the oscillation of radiation perturbation.

(C) After mass oscillation and super-horizon mode (phase (ia))

In this phase, rather than dealing with Eq. (A44) directly we adopt an alternative approach to find the asymptotic solution. This is because the Green’s function integral of Eq. (A44) is quite cumbersome after mass oscillation. Before proceeding, the background fields $\Psi_{r}$ and $\Psi_{s}$ in Eqs. (A24) and (A25) are redefined by dividing a factor $2^{1/4}\Gamma(5/4)$, i.e., $\Psi_{r}=J_{1/4}(u_{m})/u_{m}^{1/4}$ and $\Psi_{s}=Y_{1/4}(u_{m})/u_{m}^{1/4}$, to avoid irrelevant numerical factors appearing in the derivation. With this definition, the particular integral $\psi_{p}$ will also be divided by the same factor. This change does not affect the values of either $\delta_{\psi}$ or $\Delta_{\psi}$.

First, $\psi_{p}$ can be expanded as follows,

[TABLE]

for small $A_{k}$. Similarly, the same expansion applies to the driving source $S$,

[TABLE]

The metric perturbation $\phi(u_{k})$ inside $S(u_{m},A_{k})$ can also be expanded as

[TABLE]

Substituting above relations into the perturbed field equation, it follows that,

[TABLE]

where $\widetilde{S}^{(n)}=S^{(n)}-\psi_{p}^{(n-1)}/u_{m}$ with $\psi_{p}^{(-1)}(u_{m})=0$ and $\psi_{p}^{(n)}(0)=d\psi_{p}^{(n)}/du_{m}(0)=0$ for any $n=0,1,2,3,\cdots$. It turns out that $\psi_{p}^{(n)}$ satisfies

[TABLE]

with $W(\Psi_{r},\Psi_{s})(u_{m})=2u_{m}^{-3/2}/\pi$, which resembles Eq. (A44) since $\psi\rightarrow\Psi$ as $k\rightarrow 0$ or $A_{k}\rightarrow 0$.

Now, $\widetilde{S}^{(0)}(u_{m})=-2\phi_{0}\Psi(u_{m})$. Substituting it into Eq. (A62), one obtains the two integrals in Eq. (A62) for $n=0$,

[TABLE]

and

[TABLE]

where ${}_{1}F_{2}$ is the generalized hypergeometric function. Due to the fact $u_{m}\gg 1$ in this phase, we can have the following asymptotic forms,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with $\Theta=u_{m}-3\pi/8$. After straightforward substitution, $\psi_{p}^{(0)}$ becomes

[TABLE]

It is difficult to get an analytic form for $\psi_{p}^{(n)}$ with $n=1,2,3,\cdots$. However, we can take advantage of $u_{m}\gg 1$, and separate the integral $\int_{0}^{u_{m}}=\int_{0}^{u_{m}^{0}}+\int_{u_{m}^{0}}^{u_{m}}$ where $u_{m}^{0}$ is a constant $\gg 1$. Decompose $\psi_{p}^{(n)}$ as follows:

[TABLE]

The first two integrals give constants of order $O((u_{m}^{0})^{n})$, but the last two integrals are at least order of $O(u_{m}^{n})\gg O((u_{m}^{0})^{n})$. Now, the last two integrations can be carried out via the above approximate formula (Eqs. (A60), (A65), (A66) and (A69)). The leading order forms of $\psi_{p}^{(n)}$ in the expansion of $u_{m}^{-1}$ can be obtained as follows:

[TABLE]

Substituting Eqs. (A69) and (A71) into Eq. (A58) and summing up all order $A_{k}$, we arrive at

[TABLE]

where $a_{l}$, $b_{l}$ and $c_{l}$ are constants which make the infinite summation in Eq. (A72) converges for $u_{k}\leq O(1)$. Here, we identify the coefficient of $\Psi_{r}$ to be $\Re[\hat{\psi}]$ and that of $\Psi_{s}$ to be $\Im[\hat{\psi}]$ discussed in phase (ia) of the main text.

The energy density perturbation $\delta_{\psi}$ can be express as

[TABLE]

Similarly, the coefficient $d_{l}$ plays the same role with $a_{l}$, $b_{l}$ and $c_{l}$. From Eq. (A39), one can verify $\Delta_{C}\rightarrow 0$ as $u_{k}\rightarrow 0$. This is consistent with the adiabatic condition discussed in the main text.

Finally, the asymptotic form of the gauge covariant energy density perturbation $\Delta_{\psi}$ is followed by adding the last term of the numerator in Eq. (A45) to Eq. (A73). This additional term exactly cancel the second term in Eq.(A73) up to the leading order. Hence the gauge covariant energy density perturbation becomes

[TABLE]

Again, the constant $e_{l}$ makes the infinite summation in Eq. (A74) converges when $u_{k}\leq O(1)$. It recovers the CDM result (Eq.(A46)) to the leading term. The high order terms can be ignored as long as $A_{k}\ll 1$ or $u_{k}\leq O(1)$. Since $u_{k}\ll 1$ and $u_{m}\gg 1$ in this phase, this implies $A_{k}\ll 1$. Hence, the leading term of Eq. (A74) is a good approximation, and the gauge covariant energy density perturbation of $\psi$DM follows the CDM case and grows as $u_{k}^{2}\propto a^{2}$.

(D) After mass oscillation and sub-horizon mode (phase (ib))

The order of $A_{k}$ can be arbitrary in this phase and so high order terms in Eq. (A74) cannot be ignored. We shall deal with Eq. (A44) directly. Noted that

[TABLE]

where $C_{1}$ and $C_{2}$ are constants depending only on $A_{k}$. This simplification is based on that the driving source $S(u_{m})$ diminishes as $u_{m}\rightarrow\infty$, and the integrals in Eq.(A75) can be approximated by constants as long as $u_{k}\geq O(1)$ and $u_{m}\gg 1$.

From the asymptotic relations for $\psi_{r}$ and $\psi_{s}$,

[TABLE]

[TABLE]

one can obtain $\psi_{p}$ from Eq. (A44) in terms of the above expressions of $\psi_{r}$ and $\psi_{s}$ as well as $C_{1}(A_{k})$ and $C_{2}(A_{k})$. A straightforward calculation from Eq. (A45) shows that the gauge covariant energy density perturbation can be expressed as

[TABLE]

where the mass oscillation of $\psi_{p}$ is eliminated by that of $\Psi_{r}$ and it is left with only the matter-wave oscillation. Here $B_{r}(A_{k})$ and $B_{s}(A_{k})$ are the linear combination of $C_{1}(A_{k})$ and $C_{2}(A_{k})$. The gauge-covariant energy density oscillates with a constant frequency $A_{k}$ in $\ln(u_{k})$ space.

Note that Equation (A78) is valid for any $A_{k}$. One can identify constants $B_{r}$ and $B_{s}$ in the long-wave limit, i.e., $A_{k}\ll 1$, where Eq. (A78) can be approximated as $\Delta_{\psi}\sim B_{r}(A_{k})+B_{s}(A_{k})A_{k}\ln(u_{k})$. The coefficients $B_{r}(A_{k})$ and $B_{s}(A_{k})$ can be fixed by solution matching, which demands that Eq. (A78) with $u_{k}\rightarrow 1$ from above should match Eq. (A74) with $u_{k}\rightarrow 1$ from below. Comparing these two equations, one can find

[TABLE]

This expression agrees with the result of Phase (ib) in the main text.

Appendix B Full treatment Evolution

Governing equations are presented here with considerations of neutrino decoupling and photon-baryon coupling. For the zeroth order equations, the Friedmann equations with all relevant species are

[TABLE]

[TABLE]

In addition, all species obey their own conservation laws (Eq. (A7)).

To enclose above equations, the equation of state is specified for the neutrino and baryon. Since neutrinos are relativistic particles, the equation of state is the same with photons, i.e., $P_{\nu}=\epsilon_{\nu}/3$. Baryons are non-relativistic and described by a pressureless fluid, i.e., $P_{b}=0$.

First-order perturbed Einstein equations have the same structures as Friedmann equations, i.e.,

[TABLE]

[TABLE]

We need the continuity (density) equation and Euler (momentum) equation for every species to close perturbed equations. For $\psi$DM, it obeys Eqs. (A14) and (A15). On the other hand, neutrinos have decoupled from other species and become collisionless particles very early in the radiation-dominant era. Free-streaming damping will make its energy density perturbation diminished abruptlyMa, and Bertschinger (1995) upon the mode entering horizon. However, before entering horizon, neutrino perturbations have the same behavior as photon perturbations under the adiabatic condition. Hence the energy density and velocity potential perturbations of neutrinos can be approximated as,

[TABLE]

[TABLE]

where $\tau_{1}$ is the conformal time of the first oscillation null of $\Delta_{ph}(=\delta_{ph}-4H\theta_{ph})$ and is $k$-dependentMa, and Bertschinger (1995).

The photon and baryon density and momentum equations read:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the subscripts ”ph” and ”b” stand for photon and baryon respectively, $Q\equiv\theta_{ph}-\theta_{b}$, $\beta=(4\epsilon_{ph}/3\epsilon_{b})\alpha$ and $\alpha=l_{T}^{-1}$ with $l_{T}$ being the Thomson scattering mean free path. The last quantity is defined to be $l_{T}^{-1}\equiv an_{e}\sigma_{T}$, where $n_{e}$ is the electron number density and $\sigma_{T}$ the Thomson cross-section. The photon and baryon fluids are coupled through the relative velocity potential $Q$ via Thomson scattering.

Note that $l_{T}^{-1}$ decreases with time as $a^{-2}$ for fully ionized baryons, and at some point $l_{T}^{-1}$ can be smaller than $k$. After then, photons can no longer be described by a fluid and the above photon fluid equations fail to hold. Hence, for the analysis below to be valid we must demand $kl_{T}\ll 1$. Moreover, we shall be focusing on modes that have already entered horizon, $k/H\gg 1$. It turns out that $k_{c}$ of $m=10^{-22}$ eV satisfies the first criteria throughout the radiation dominant era and satisfies the second in the late epoch of the radiation era.

In the sub-horizon regime, the photon momentum equation is dominated by the photon pressure and the coupling term, and hence the gravity can be ignored. This situation is similar to the short wave limit of the self-gravitating nonrelativistic fluid where perturbations do not feel the gravity and become sound waves. We therefore drop $\phi$ in Eqs. (B7), (B8), (B9) and (B10).

Subtracting Eq. (B8) from Eq. (B10), we can derive an equation for $Q$:

[TABLE]

Substituting $(1/4)\delta_{ph}$ from Eq. (B11) into Eq. (B8) and solving for $\theta_{ph}$ and $\theta_{ph}^{{}^{\prime}}$, we arrive at

[TABLE]

[TABLE]

Now, since we are in the regime $k\gg H$, any term of order $O(H^{2}/k^{2})$ or $O(H^{{}^{\prime}}/k^{2})$ can be ignored. However, we need to keep terms of order $O(H/k)$ as will become clear later. The reason for the coupling terms, such as $\alpha Q$ or $\beta Q$, to be finite leading-order terms is that we have large $\alpha$ and $\beta$ despite having a small relative velocity potential $Q$. Therefore, $Q^{{}^{\prime}}$ is small compared to $\alpha Q$ or $\beta Q$, as $Q^{{}^{\prime}}=O((k/\sqrt{3})Q)$ due to the photon oscillation to be shown below and furthermore $kl_{T}\ll 1$.

After the above considerations, we find

[TABLE]

[TABLE]

Terms in the brace brackets are terms of order $O(H/k)$. If we ignore those high-order terms in Eqs. (B14) and (B15), it follows

[TABLE]

The quantity $Q^{{}^{\prime}}+(\alpha+\beta)Q$ oscillates with the photon sound frequency modified by the presence of baryons, and this quantity is nothing more than $-(1/4)\Delta_{ph}$ from Eq. (B11), where $\Delta_{ph}$ is the covariant energy density perturbation of photon. Thus the leading order terms give rise to a fast photon sound oscillation, and so the slow evolution of drag damping entails the next order correction to Eq. (B16), which all involve $Q^{{}^{\prime}}$. To the leading order, $Q^{{}^{\prime}}=ik(\beta/(3(\alpha+\beta)))^{1/2}Q$.

The approach by which we will derive the next order equation is similar to Phase (i) of Sec. (III), where we let $Q^{{}^{\prime}}+(\alpha+\beta)Q=A(\tau)\exp[ik\int(\beta/(3(\alpha+\beta)))d\tau]$ and $A(\tau)$ is a slowly varying function of $\tau$. Substituting this expression into Eqs. (B14) and (B15), we find that

[TABLE]

where $\eta\equiv\ln a$. This equation has an analytical solution.

It is clear from that Eq. (B17) that for very long waves, $k^{2}\to 0$, the amplitude $A$ of $\Delta_{ph}$ will decrease when $\alpha/\beta(=(3/4)\epsilon_{b}/\epsilon_{ph})$ becomes non-negligible occurring shortly after radiation-matter equality. On the other hand for shorter waves where $k^{2}>H(\alpha+\beta)$, the amplitude $A$ will further decrease and occur earlier in time. We find that $k_{c}^{2}/H\beta\sim 1$ shortly before the radiation-matter equality and the drag damping begins to take effect for the critical mode near the end of radiation-dominant era.

Figure (7) shows the evolution of $\Delta_{ph}$ for $k=k_{c}$ and $k=0.1k_{c}$, where both solutions of Eq. (B17) and of the full treatment are given. Good agreement is found. Baryon density perturbations can also be calculated from Eqs. (B9) and (B14) with an understanding that $\Delta_{b}=\delta_{b}$ and $\theta_{ph}=\theta_{b}$ to the leading order. It yields $\Delta_{b}=(3/4)\Delta_{ph}$. We also plot baryons solutions of the above approximate treatment and the full treatment. Again, good agreement is found, till after radiation-matter equality, as afterwards the gravity of matter can no longer be ignored.